Numpy Correlate没有提供补偿
问题内容:
我正在尝试使用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')
此代码的输出如下所示:
请注意,尽管模板(蓝色)和观察到的(红色)光谱清楚地显示了偏移,但互相关函数在零像素偏移处仍达到峰值。我希望看到的有点像(尽管不完全一样;这只是我可以产生的最接近的表示形式):
在这里,我在模板数据中引入了50个像素的人工偏移,现在它们或多或少地排列在一起。我想要的是,在这种情况下,峰出现在50个像素的偏移处而不是零处(我不在乎底部的光谱是否排成一行;这只是为了视觉表示)
。但是,尽管经过数小时的工作和在线研究,我找不到一个描述这个问题的人,更不用说解决方案了。我尝试使用ScyPy的相关函数和MatLib的xcorr,并且bot显示了相同的内容(尽管我被认为是本质上相同的功能)。
为什么互相关不按我期望的方式运行,如何使它以有用的方式运行?
问题答案:
您遇到的问题可能是因为您的光谱不是零中心的。无论您以哪种单位绘制,其RMS值均约为100。之所以这样,是因为卷积/互相关函数必须 在频谱上填充零
,以便在“相同”模式下计算出完整的响应。因此,即使您的信号与50个样本之间的偏移量最相似,但当两个信号没有完全对齐时,您将仅对其重叠量进行乘积积分,并丢弃所有偏移量值,因为它们乘以零。这是有问题的,因为您的频谱不是零均值的,并且它们的相关性在重叠中几乎呈线性增加。
请注意,您的互相关结果看起来像一个三角形脉冲,这是您从两个方形脉冲的互相关中所期望的结果(请 参见
矩形“脉冲”与其自身的卷积。这是因为一旦填充后的光谱看起来像一个从零到大约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')
,我们以“垫” 摹
与在零一半的长度两侧(默认情况下,你可以与其他值垫)。如果我们 不要t 填充 g (如所示),因为 f 和 g
的长度相同,并且没有空间将一个信号相对于另一个移位,所以np.correlate(g,f,mode='valid')
我们只会得到 一个
值作为回报(零偏移的相关性)。
在填充之后计算 g 和 f 的相关性时,您会发现当信号 的非零 部分完全对齐时,即原始 f 和 g 之间 没有偏移
时,它达到峰值。这是因为信号的RMS值远高于零 -f 和 g 重叠的大小 __
__与在每个RMS周围具有相对较小的波动相比,在此较高RMS级别上重叠的元素数量要更重要。
我们可以通过从每个序列中减去RMS水平来消除对相关性的巨大贡献。在下图中,右侧的灰线显示了零位居中之前两个序列的互相关,而蓝绿色线显示了之后的互相关。像您第一次尝试一样,灰线是三角形的,两个非零信号重叠。如我们所愿,蓝绿色线更好地反映了两个信号波动之间的相关性。
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()