Python 与直方图和高斯卷积

Question

我有一个显示为直方图的模拟信号。我想使用具有特定宽度的高斯卷积来模拟真实测量信号，因为在真实实验中，检测器在测量通道中具有一定的不确定性。

我尝试使用 np.convolve 和 scipy.signal.convolve 进行卷积，但似乎无法正确进行过滤。不仅预期的形状被关闭，这将是直方图和 x 轴的轻微涂抹版本，例如能量等级也关闭了。

我尝试将宽度为 20 keV 的高斯定义为：

gauss = np.random.normal(0, 20000, len(coincidence['esum']))
hist_gauss = plt.hist(gauss, bins=100)[0]

其中 len(coincidence['esum']) 是我的 coincidence 数据帧 column.This 列的长度，我使用：

counts = plt.hist(coincidence['esum'], bins=100)[0]

除了这种生成合适高斯分布的方法外，我还尝试了 scipy.signal.gaussian(50, 30000)，不幸的是，它生成了 抛物线形 曲线并且没有表现出特征尾部。

我尝试使用 coincidence['esum'] 和 counts 两种高斯方法进行卷积。请注意，当根据 Finding the convolution of two histograms 与标准示例进行简单卷积时，它可以正常工作。

有人知道如何在 python 中做这样的卷积吗？我将用于直方图的 coincidende['esum'] 列导出到一个 pastebin，以防有人感兴趣并想用特定数据 https://pastebin.com/WFiSBFa6

重新创建它

Answer 1

如您所知，对具有相同 bin 大小的两个直方图进行卷积将给出将一个样本的每个元素与另一个样本的每个元素相加的结果的直方图。

我看不清楚你在做什么。您似乎没有做的一件重要事情是确保直方图的箱子具有相同的宽度，并且您必须注意第二个箱子边缘的位置。

在代码中我们有

def hist_of_addition(A, B, bins=10, plot=False):
    A_heights, A_edges = np.histogram(A, bins=bins)
    # make sure the histogram is equally spaced
    assert(np.allclose(np.diff(A_edges), A_edges[1] - A_edges[0]))
    # make sure to use the same interval
    step = A_edges[1] - A_edges[0]
    
    # specify parameters to make sure the histogram of B will
    # have the same bin size as the histogram of A
    nBbin = int(np.ceil((np.max(B) - np.min(B))/step))
    left = np.min(B)
    B_heights, B_edges = np.histogram(B, range=(left, left + step * nBbin), bins=nBbin)
    
    # check that the bins for the second histogram matches the first
    assert(np.allclose(np.diff(B_edges), step))
    
    C_heights = np.convolve(A_heights, B_heights)
    C_edges = B_edges[0] + A_edges[0] + np.arange(0, len(C_heights) + 1) * step
    
    if plot:
        plt.figure(figsize=(12, 4))
        plt.subplot(131)
        plt.bar(A_edges[:-1], A_heights, step)
        plt.title('A')
        plt.subplot(132)
        plt.bar(B_edges[:-1], B_heights, step)
        plt.title('B')
        plt.subplot(133)
        plt.bar(C_edges[:-1], C_heights, step)
        plt.title('A+B')
    return C_edges, C_heights

然后

A = -np.cos(np.random.rand(10**6))
B = np.random.normal(1.5, 0.025, 10**5)
hist_of_addition(A, B, bins=100, plot=True);

给予

Python 与直方图和高斯卷积

Python convolution with histogram and Gaussian

python

numpy

scipy

signal-processing

convolution