Numpy Correlate 不提供偏移量
Numpy Correlate is not providing an offset
我正在尝试使用 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,并且机器人显示了同样的东西(尽管我被引导相信它们本质上是相同的功能)。
为什么互相关没有按我预期的方式运行,如何让它以有用的方式运行?
您遇到的问题可能是因为您的光谱不zero-centered;无论您绘制的是什么单位,它们的 RMS 值看起来都在 100 左右。这是一个问题的原因是 convolution/cross-correlation 函数必须 用零 填充您的光谱,以便在 "same" 模式下计算完整响应。因此,即使您的信号最相似且偏移量约为 50 个样本,但当两个信号未完全对齐时,您只是对它们重叠的乘积进行积分,并丢弃所有偏移值,因为它们已乘以零。这是有问题的,因为您的光谱不是 zero-mean,并且它们的相关性在重叠时几乎呈线性增加。
请注意,您的 cross-correlation 结果看起来像一个三角形脉冲,这正是您从两个方形脉冲 (c.f. Convolution of a Rectangular "Pulse" With Itself。那是因为你的光谱,一旦被填充,看起来就像一个从零到 100 左右的轻微噪声值脉冲的阶梯函数——实际上是矩形脉冲与高斯噪声的卷积。你可以尝试与 mode='full'
进行卷积以查看您正在关联的两个光谱的整个响应,或者,请注意,对于 mode='valid'
您应该只获得 一个值 return,因为你的两个光谱长度完全相同,所以只有 一个偏移量 (零!),你可以将它们完全排列起来。
为了回避这个问题,您可以尝试减去光谱的 RMS 值,使它们 zero-centered,或者在两侧的 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')
,我们必须 "pad" g边的一半长度为零(默认情况下;您可以填充其他值。)如果我们 不 填充 g(如 np.correlate(g,f,mode='valid')
),我们只会在return中得到一个值(与零偏移量的相关),因为f和g 是相同的长度,并且没有一个信号相对于另一个信号移位的空间。
当你在填充后计算 g 和 f 的相关性时,你会发现它在 时达到峰值non-zero 部分信号完全对齐,即当原始 f 和 f 之间没有偏移 时=39=]g。这是因为信号的 RMS 值远高于零——f 和 g 的重叠大小依赖性更强在这个高 RMS 水平上重叠的元素数量比在每个函数周围相对较小的波动上。我们可以通过从每个系列中减去 RMS 水平来消除对相关性的巨大贡献。在下图中,右侧的灰线显示 cross-correlation 之前的两个系列 zero-centering,蓝绿色线显示 cross-correlation 之后的两个系列。灰线与您的第一次尝试一样,是两个 non-zero 信号重叠的三角形。正如我们所期望的那样,蓝绿色线更好地反映了两个信号波动之间的相关性。
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()
我正在尝试使用 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,并且机器人显示了同样的东西(尽管我被引导相信它们本质上是相同的功能)。
为什么互相关没有按我预期的方式运行,如何让它以有用的方式运行?
您遇到的问题可能是因为您的光谱不zero-centered;无论您绘制的是什么单位,它们的 RMS 值看起来都在 100 左右。这是一个问题的原因是 convolution/cross-correlation 函数必须 用零 填充您的光谱,以便在 "same" 模式下计算完整响应。因此,即使您的信号最相似且偏移量约为 50 个样本,但当两个信号未完全对齐时,您只是对它们重叠的乘积进行积分,并丢弃所有偏移值,因为它们已乘以零。这是有问题的,因为您的光谱不是 zero-mean,并且它们的相关性在重叠时几乎呈线性增加。
请注意,您的 cross-correlation 结果看起来像一个三角形脉冲,这正是您从两个方形脉冲 (c.f. Convolution of a Rectangular "Pulse" With Itself。那是因为你的光谱,一旦被填充,看起来就像一个从零到 100 左右的轻微噪声值脉冲的阶梯函数——实际上是矩形脉冲与高斯噪声的卷积。你可以尝试与 mode='full'
进行卷积以查看您正在关联的两个光谱的整个响应,或者,请注意,对于 mode='valid'
您应该只获得 一个值 return,因为你的两个光谱长度完全相同,所以只有 一个偏移量 (零!),你可以将它们完全排列起来。
为了回避这个问题,您可以尝试减去光谱的 RMS 值,使它们 zero-centered,或者在两侧的 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')
,我们必须 "pad" g边的一半长度为零(默认情况下;您可以填充其他值。)如果我们 不 填充 g(如 np.correlate(g,f,mode='valid')
),我们只会在return中得到一个值(与零偏移量的相关),因为f和g 是相同的长度,并且没有一个信号相对于另一个信号移位的空间。
当你在填充后计算 g 和 f 的相关性时,你会发现它在 时达到峰值non-zero 部分信号完全对齐,即当原始 f 和 f 之间没有偏移 时=39=]g。这是因为信号的 RMS 值远高于零——f 和 g 的重叠大小依赖性更强在这个高 RMS 水平上重叠的元素数量比在每个函数周围相对较小的波动上。我们可以通过从每个系列中减去 RMS 水平来消除对相关性的巨大贡献。在下图中,右侧的灰线显示 cross-correlation 之前的两个系列 zero-centering,蓝绿色线显示 cross-correlation 之后的两个系列。灰线与您的第一次尝试一样,是两个 non-zero 信号重叠的三角形。正如我们所期望的那样,蓝绿色线更好地反映了两个信号波动之间的相关性。
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()