使用Numpy/Scipy相关函数识别正弦余弦信号时移结果异常

问题描述 #

尝试用信号相关法识别两个正弦信号的时移：

t = np.arange(0, 10, 0.01)
sig1 = np.sin(t)
sig2 = np.cos(t)

理论上sig2相对sig1的时移为π/2≈1.571，但使用np.correlate或scipy.signal.correlate计算得到的时移为|1.490|，而余弦相似度、最小二乘优化等方法能得到正确结果，怀疑是mode参数使用有误。

问题根源 #

并非mode参数的问题，核心原因有两个：

• 信号截断的边界效应：选取的信号片段（0到10秒）不是正弦信号周期的整数倍（正弦周期≈6.28秒），截断后的信号两端不连续，导致相关计算时的边界匹配误差，使得相关峰值偏离理论位置。
• 时间滞后轴计算错误：原代码用np.linspace(-t[-1], t[-1], 2*len(t)-1)生成滞后时间轴，实际采样步长是0.01秒，但linspace生成的间隔并非严格等于采样步长（计算得≈0.009995秒），引入了额外数值误差。

解决方案 #

1. 使用完整周期的信号片段 #

选取整数个完整周期的信号，消除截断边界的影响：

import numpy as np

dt = 0.01
# 取2个完整周期的信号
t = np.arange(0, 2*np.pi*2, dt)
sig1 = np.sin(t)
sig2 = np.cos(t)

2. 正确计算时间滞后轴 #

基于采样步长和相关滞后索引生成时间轴，而非linspace：

Numpy实现

x = np.correlate(sig1, sig2, mode='full')
# 生成正确的滞后索引：从-(len(sig1)-1)到len(sig2)-1
lags = np.arange(-(len(sig1)-1), len(sig2), 1)
tx = lags * dt

ndx = np.argmax(np.abs(x))
print(f"时移结果：{tx[ndx]:.5f}")
# 输出：时移结果：-1.57000（接近π/2≈1.5708，误差来自采样离散性）

Scipy实现

利用scipy.signal.correlation_lags直接生成滞后索引，再乘以采样步长：

from scipy import signal

x = signal.correlate(sig1, sig2, mode='full')
lags = signal.correlation_lags(len(sig1), len(sig2))
tx = lags * dt

ndx = np.argmax(np.abs(x))
print(f"时移结果：{tx[ndx]:.5f}")

3. 优化峰值检测（可选） #

对于非完整周期信号，可对相关结果进行插值，提升峰值位置精度：

from scipy.interpolate import interp1d

# 对相关曲线进行高分辨率插值
f = interp1d(tx, x, kind='cubic')
tx_highres = np.linspace(tx[ndx-5], tx[ndx+5], 1000)
x_highres = f(tx_highres)
ndx_high = np.argmax(np.abs(x_highres))
print(f"插值后时移：{tx_highres[ndx_high]:.5f}")

修正后的完整测试代码 #

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

dt = 0.01
# 使用2个完整周期的信号
t = np.arange(0, 2*np.pi*2, dt)
sig1 = np.sin(t)
sig2 = np.cos(t)

fig, ax = plt.subplots(3, figsize=(8,6))

# 绘制输入信号
ax[0].plot(t, sig1, label='sig1')
ax[0].plot(t, sig2, label='sig2')
ax[0].legend()
ax[0].set_xlabel('时间 (s)')

# Scipy相关计算
x = signal.correlate(sig1, sig2, mode='full')
lags = signal.correlation_lags(len(sig1), len(sig2))
tx = lags * dt

# 绘制相关曲线
ax[1].plot(tx, x, label='相关曲线', c='k')
ndx = np.argmax(np.abs(x))
ax[1].plot(tx[ndx], x[ndx], 'rx', label='峰值位置')
ax[1].legend()

# 绘制校正后的信号
ax[2].plot(t, sig1, label='sig1')
t_corr = t - tx[ndx]
ax[2].plot(t_corr, sig2, label='校正后的sig2')
ax[2].set_xlabel(f'时移：{tx[ndx]:.5f} s')
ax[2].legend()

print(f"最终时移结果：{tx[ndx]:.5f}")
plt.tight_layout()
plt.show()

内容的提问来源于stack exchange，提问作者schneider-daniel