我试着用Python来观察天文光谱,用numpy.correlate来寻找径向速度偏移。我在比较每个光谱,我有一个模板光谱。我遇到的问题是,无论使用哪种谱,NUPYP.关联状态表明,相关函数的最大值以零像素的移位发生,即光谱已经排队,这显然不是真的。以下是一些相关代码:

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像素的偏移处,而不是0像素(我不在乎底部的光谱是否排成一行;这仅仅是为了视觉表现)。然而,尽管在网上工作和研究了几个小时,我还是找不到一个能描述这个问题的人,更别说找到解决方案了。我尝试使用ScyPy的correlate和MatLib的xcorr,bot也显示了同样的东西(尽管我相信它们本质上是相同的函数)。
为什么互相关没有按我期望的方式进行,如何让它以一种有用的方式进行操作?

最佳答案

您遇到的问题可能是因为您的光谱不是以零为中心的;它们的RMS值看起来是100左右,以您正在绘制的任何单位表示。这是一个问题,因为卷积/互相关函数必须用零填充你的光谱,以便在“相同”模式下计算完整响应。因此,即使您的信号与偏移量约为50个样本的信号最为相似,但当两个信号没有完全对齐时,您将只对其重叠部分的乘积进行积分,并丢弃所有偏移量值,因为它们乘以零。这是有问题的,因为你的光谱不是零均值的,而且它们的相关性在重叠中几乎线性增加。
请注意,您的互相关结果看起来像一个三角形脉冲,这是您可能从两个正方形脉冲的互相关(c.f.Convolution of a Rectangular "Pulse" With Itself)中得到的结果。这是因为你的光谱,一旦被填充,看起来就像是一个阶跃函数,从零到一个脉冲的噪声值约为100——实际上是矩形脉冲与高斯噪声的卷积。您可以尝试用mode='full'卷积来查看您正在关联的两个光谱的整个响应,或者,请注意,用mode='valid'卷积,您应该只得到一个返回值,因为您的两个光谱的长度完全相同,所以只有一个偏移量(零!)你可以把它们排成一行。
为了避免这个问题,您可以尝试减去光谱的均方根值,使其为零居中,或者用两边的均方根值填充两个光谱的长度。
编辑:
在回答你在评论中提出的问题时,我想我会附上一个图表来说明我试图描述的要点。
假设我们有两个向量的值,和你们的光谱不完全不同,每个向量都有很大的偏移量。

# 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

python - Numpy Correlate没有提供补偿-LMLPHP
f的峰值在t=90左右,g的峰值在t=180左右。因此,我们期望g和f的相关性在90个时间步(或频率箱,或任何相关函数的参数)的滞后上有一个峰值
但是为了得到一个与我们的输入形状相同的输出,如在np.correlate(g,f,mode='same')中,我们必须用一半长度的0来“填充”g(默认情况下,你可以用其他值填充)。如果我们不填充g(如在np.correlate(g,f,mode='valid')中),我们将只得到一个值(零偏移的相关性),因为f和g是相同的长度,并且没有一个信号相对于另一个信号的移动空间。
当你计算填充后g和f的相关性时,你会发现它在信号的非零部分完全对齐时达到峰值,也就是说,当原始f和g之间没有偏移量时,这是因为信号的均方根值远远高于零——f和g的重叠部分的大小更依赖于在这个高均方根水平上重叠的元素的数量,而不是每个函数周围相对较小的波动。我们可以通过从每个序列中减去均方根水平来消除这种对相关性的巨大贡献。在下图中,右边的灰色线表示两个序列在归零之前的互相关,而teal线表示在归零之后的互相关。像你第一次尝试一样,灰线是三角形的,两个非零信号重叠。teal线更好地反映了两个信号波动之间的相关性,正如我们所希望的那样。
python - Numpy Correlate没有提供补偿-LMLPHP
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()

10-06 05:20
查看更多