python - Numpy Correlate 不提供偏移量

Question

我正在尝试使用 Python 查看天文光谱，并且我正在使用 numpy.correlate 尝试找到径向速度偏移。我正在将我拥有的每个光谱与一个模板光谱进行比较。我遇到的问题是，无论我使用哪个光谱， numpy.correlate 都表明相关函数的最大值出现零像素偏移，即光谱已经对齐，这显然不是真的.下面是一些相关的代码:

corr = np.correlate(temp_data, imag_data, mode='same')
ax1.plot(delta_data, corr, c='g')
ax1.plot(delta_data, 100*temp_data, c='b')
ax1.plot(delta_data, 100*imag_data, c='r')

此代码的输出如下所示:

What I Have

请注意，尽管模板(蓝色)和观察到的(红色)光谱清楚地显示了偏移，但互相关函数在零像素的偏移处达到峰值。我希望看到的会有点像(虽然不完全像；这只是我能产生的最接近的表示):

What I Want

在这里，我在模板数据中引入了一个 50 像素的人工偏移量，现在它们或多或少地对齐了。我想要的是，对于这样的情况，峰值出现在 50 个像素的偏移处而不是零处(我不在乎底部的光谱是否排列；这仅用于视觉表示) .但是，尽管在网上进行了几个小时的工作和研究，但我找不到甚至描述此问题的人，更不用说解决方案了。我曾尝试使用 ScyPy 的关联和 MatLib 的 xcorr，并且 bot 显示了同样的事情(尽管我相信它们本质上是相同的功能)。

为什么互相关没有按照我期望的方式运行，我该如何让它以有用的方式运行？

Answer 1

您遇到的问题可能是因为您的光谱不是以零为中心的；无论您绘制哪种单位，它们的 RMS 值看起来都约为 100。这是一个问题的原因是因为卷积/互相关函数必须用零填充您的频谱，以便在“相同”模式下计算完整响应。因此，即使您的信号最相似，偏移量约为 50 个样本，但当两个信号未完全对齐时，您仅对它们重叠的乘积进行积分，并丢弃所有偏移值，因为它们乘以零。这是有问题的，因为您的光谱不是零均值，并且它们的相关性在重叠时几乎呈线性增加。

请注意，您的互相关结果看起来像一个三角脉冲，这可能是您从两个方波脉冲(参见 Convolution of a Rectangular "Pulse" With Itself)的互相关中所期望的结果。那是因为您的频谱，一旦被填充，看起来就像一个从零到一个的阶跃函数100 左右略带噪声值的脉冲——实际上是矩形脉冲与高斯噪声的卷积。您可以尝试使用 mode='full' 进行卷积以查看您正在关联的两个光谱的整个响应，或者，请注意，使用 mode='valid' 时，您应该只得到一个值作为返回，因为您的两个光谱长度完全相同，因此您可以将它们完全对齐的偏移量(零!)。

为了避免这个问题，您可以尝试减去光谱的 RMS 值，使它们以零为中心，或者用它们在任一侧的 RMS 值中的长度填充两个光谱。

编辑 :
为了回答您在评论中提出的问题，我想我会附上一个图形来使我试图更清楚地描述这一点。

假设我们有两个值向量，与您的光谱并不完全不同，每个向量都与零有很大的偏移。

# Generate two noisy, but correlated series
t = np.linspace(0,250,250)
f = 10*np.exp(-((t-90)**2)/8) + np.random.randn(250) + 40
g = 10*np.exp(-((t-180)**2)/8) + np.random.randn(250) + 40

f 在 t=90 附近有一个尖峰，g 在 t=180 附近有一个尖峰。因此，我们希望 g 和 f 的相关性在 90 个时间步长(或频率区间，或您相关的函数的任何参数)附近有一个尖峰。

但是为了获得与我们的输入形状相同的输出，如 np.correlate(g,f,mode='same') ，我们必须在任一侧“填充” g 的一半长度为零(默认情况下；您可以填充其他值。)如果我们不要填充 g (如 np.correlate(g,f,mode='valid') )，我们只会得到 一个 值作为返回(与零偏移的相关性)，因为 f 和 g 的长度相同，并且没有空间相对移动其中一个信号到另一个。

当您计算填充后的 g 和 f 的相关性时，您会发现当信号的非零部分完全对齐时，即原始 f 和 g 之间没有偏移时，它会达到峰值。这是因为信号的 RMS 值远高于零——f 和 g 的重叠大小更依赖于在这个高 RMS 级别重叠的元素数量，而不是每个函数具有的相对较小的波动周围。我们可以通过从每个系列中减去 RMS 水平来消除对相关性的巨大贡献。在下图中，右侧的灰线表示两个系列归零前的互相关，青色线表示归零后的互相关。灰线与您的第一次尝试一样，是两个非零信号重叠的三角形。正如我们所希望的，青色线更好地反射(reflect)了两个信号波动之间的相关性。

xcorr = np.correlate(g,f,'same')
xcorr_rms = np.correlate(g-40,f-40,'same')
fig, axes = plt.subplots(5,2,figsize=(18,18),gridspec_kw={'width_ratios':[5,2]})
for n, axis in enumerate(axes):
    offset = (0,75,125,215,250)[n]
    fp = np.pad(f,[offset,250-offset],mode='constant',constant_values=0.)
    gp = np.pad(g,[125,125],mode='constant',constant_values=0.)

    axis[0].plot(fp,color='purple',lw=1.65)
    axis[0].plot(gp,color='orange',lw=lw)
    axis[0].axvspan(max(125,offset),min(375,offset+250),color='blue',alpha=0.06)
    axis[0].axvspan(0,max(125,offset),color='brown',alpha=0.03)
    axis[0].axvspan(min(375,offset+250),500,color='brown',alpha=0.03)
    if n==0:
        axis[0].legend(['f','g'])
    axis[0].set_title('offset={}'.format(offset-125))


    axis[1].plot(xcorr/(40*40),color='gray')
    axis[1].plot(xcorr_rms,color='teal')
    axis[1].axvline(offset,-100,350,color='maroon',lw=5,alpha=0.5)
    if n == 0:
        axis[1].legend(["$g \star f$","$g' \star f'$","offset"],loc='upper left')

plt.show()