使用FFT计算两信号零滞后互相关结果错误的原因排查

你的实现结果和预期值存在巨大差异，核心是三个层面的错误：

具体错误点 #

• 跳过逆FFT步骤，未做对应归一化
  基于FFT计算互相关的标准流程是：对两信号的FFT结果取共轭乘积后，必须经过逆FFT变换回时域，才能得到不同滞后位置的互相关值。numpy.fft的正变换无归一化系数，逆变换自带1/N（N为FFT计算长度）的归一化，你直接对频域数组做内积，得到的是N倍的零滞后循环互相关值，结果自然比正确值大3个数量级（测试用例中N=1000），还会残留浮点计算带来的无意义虚部。
  实际上你得到的-353199.837除以1000之后数值量级就和正确值匹配，符号异常可以检查代码中相位偏移、共轭相乘的顺序是否写反。
• 混淆循环互相关与线性互相关的计算要求
  直接对原长信号做FFT相乘得到的是循环互相关结果，而scipy.signal.correlate默认计算的是full模式的线性互相关。如果要和scipy的FFT计算结果完全对齐，做FFT前需要将两个信号补零到长度至少为len(y1)+len(y2)-1=1999，避免循环移位带来的混叠误差。
• 误读scipy的计算结果
  你观察到的475左右的数值是互相关的峰值，不是零滞后位置的结果。测试用例中y2比y1滞后pi/4相位，对应1000个采样点下滞后约125个点，零滞后位置的互相关理论值约为353.5，和时域直接计算np.sum(y1*y2)的结果一致。

正确实现参考 #

1. 时域直接计算零滞后互相关（最贴合定义，用于原理验证） #

import numpy as np
# 关闭端点避免0和2pi重复采样引入误差
x1 = np.linspace(0, 2*np.pi, 1000, endpoint=False)
y1 = np.sin(x1)
y2 = np.sin(x1 - np.pi/4)
zero_lag_direct = np.sum(y1 * y2)
print(zero_lag_direct) # 输出约353.5，和理论计算值匹配

2. FFT方式计算全滞后线性互相关，取零滞后值 #

N = len(y1)
L = 2 * N - 1 # 线性互相关要求的最小FFT长度
# 对信号补零避免循环混叠
y1_pad = np.pad(y1, (0, L-N))
y2_pad = np.pad(y2, (0, L-N))
# 频域计算
y1_fft = np.fft.fft(y1_pad)
y2_fft = np.fft.fft(y2_pad)
ccf_fft = np.fft.ifft(y1_fft * np.conj(y2_fft)).real
# 移位对齐滞后轴，和scipy返回的lags顺序对应
ccf_fft = np.fft.fftshift(ccf_fft)
zero_lag_fft = ccf_fft[L//2]
print(zero_lag_fft) # 输出和时域计算结果一致

3. scipy接口结果验证 #

from scipy import signal as sg
corr = sg.correlate(y1, y2, method='fft')
lags = sg.correlation_lags(N, N)
zero_lag_scipy = corr[lags == 0][0]
print(zero_lag_scipy) # 和上述两种方法结果完全匹配

内容的提问来源于stack exchange，提问作者Pasindu Samaranayake