如何对非严格周期随机微变信号应用FFT提取目标频率

处理严格周期信号时，直接调用scipy.fft做快速傅里叶变换即可精准提取频率，效果稳定。但面对存在轻微随机周期波动的非严格周期信号时，直接对整段长信号做FFT会出现大量杂乱频谱峰值，无法定位到靠近理论值1/period的目标频率范围。
原有实现代码如下：

df = pd.read_csv('data.txt', header=None)
x = np.asarray(df.iloc[:, 0])
y = np.asarray(df.iloc[:, 1])

yf = fft(y)
xf = fftfreq(len(y))

plt.plot(xf, abs(yf))

这类带周期抖动的信号示例可参考对应公开项目的说明文档。 #

可直接落地的解决方法 #

首先要修正原有代码的一个基础问题：fftfreq调用时没有传入实际采样间隔，默认采样间隔为1，输出的是归一化频率，和实际物理频率不对应，这本身就会导致峰值定位偏差，需要把x轴的平均时间差作为d参数传入。
以下几种方法都可以解决周期抖动导致的频谱杂峰问题，按实现难度从低到高排序：

• 加窗FFT+局部功率加权
  周期轻微抖动会让原本集中的单根谱线扩散成一个窄峰簇，不要直接找频谱单点最大值。先根据预估的理论周期算出中心频率，划定小范围的搜索区间，对区间内的频谱功率做加权平均就能得到稳定的等效基频；同时给时域信号加汉宁窗，可以抑制旁瓣泄露，大幅减少无关杂峰。
  核心实现代码：

  from scipy.signal import windows
  # 加汉宁窗抑制旁瓣
  win = windows.hann(len(y))
  y_win = y * win
  yf = fft(y_win)
  # 传入实际采样间隔，得到物理频率轴
  dt = np.mean(np.diff(x))
  xf = fftfreq(len(y), d=dt)
  # 替换为你预估的理论基频
  f0_guess = 1 / period_guess
  # 根据实际抖动幅度调整搜索范围，通常±10%~±20%足够覆盖轻微抖动
  mask = (xf >= f0_guess*0.9) & (xf <= f0_guess*1.1)
  # 功率加权计算实际基频
  power = np.abs(yf[mask])**2
  f0_real = np.sum(xf[mask] * power) / np.sum(power)
  

• 自相关法提取平均周期
  自相关函数对周期抖动的抗干扰能力远优于直接FFT：计算去均值后信号的自相关序列，跳过0延迟点后找到第一个峰值对应的时间延迟，就是信号的平均周期，取倒数即为平均基频，不会出现FFT那种谱线扩散导致的大量杂峰。
  核心实现代码：

  def calc_autocorr(sig):
      sig = sig - np.mean(sig)
      ac = np.correlate(sig, sig, mode='full')
      return ac[len(ac)//2:]
  
  acf = calc_autocorr(y)
  dt = np.mean(np.diff(x))
  min_lag = int(0.5 * period_guess / dt)
  max_lag = int(1.5 * period_guess / dt)
  # 在理论周期附近的区间找自相关峰值
  peak_idx = np.argmax(acf[min_lag:max_lag]) + min_lag
  period_real = peak_idx * dt
  f0_real = 1 / period_real
  

• 短时傅里叶变换跟踪瞬时频率
  如果周期波动是随时间缓慢变化的，可以把长信号切分成多个短时间窗，保证每个时间窗内信号可近似为严格周期信号，对每个窗单独做FFT提取峰值，就能得到随时间变化的瞬时频率序列，既可以算平均频率，也能量化周期的波动幅度。注意时间窗长度要短于周期波动的时间尺度，同时包含至少3~5个完整周期，保证频率分辨率。
• 窄带滤波+过零检测
  如果信号信噪比足够，可以先设计一个覆盖目标频率范围的窄带带通滤波器，用零相位滤波（调用scipy.signal.filtfilt实现，避免相位偏移）滤除带外噪声，之后统计信号的过零点间隔，就能得到每个周期的实际长度，适合需要同时获取周期波动统计特征的场景。

注意：如果信号的周期抖动幅度过大（超过理论周期的20%），属于准周期信号范畴，建议先做趋势项去除、去直流等预处理，再选择上述方法计算。

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