重采样FFT中的频谱泄漏问题及解决方案探究

问题背景 #

我有一段时长10秒、频率50Hz的纯正弦信号，采样率为40.96kSa/s（共409600个样本）。施加Hanning窗后，numpy FFT的幅度谱（忽略负镜像分量）在49.9Hz、50.0Hz、50.1Hz处出现3个谱峰，符合预期。

随后对信号进行重采样（抽取每第50个样本），得到8192个样本（等效采样率819.2Sa/s），此时FFT幅度谱在距50Hz约n*5Hz（n为整数）的位置出现频谱泄漏。观察到提升重采样频率、增加信号时长或补零、改用Flattop窗可降低泄漏幅度。

核心问题 #

1. 该频谱泄漏的成因是什么？
2. 为何将重采样信号的FFT样本数改为8193时泄漏完全消失？该方法是否具备实用性？
3. 能否在保持FFT样本数为2^n（n为整数，满足C++实现要求）的前提下避免此类泄漏？

可复现最小验证代码 #

import numpy as np
import plotly.express as px


# ___ 输入参数 ___
m = 1  # 信号幅度 (-)
H = 1  # 谐波次数 (-)
phi = 0  # 相移 (°)
frequency_H1 = 50  # 基波频率 (Hz)
T_max = 10  # 信号时长 (s)
sampling_rate = 40960  # 原始信号采样率 (Sa/s)
fft_size_res = int(8192)  # 重采样信号的FFT点数 (-)
fft_y_axis_limit = 0.01  # FFT幅度图的y轴上限 (-)
# ___ 衍生参数 ___
omega = 2*np.pi*frequency_H1
t = np.arange(0, T_max, 1/sampling_rate)


# 1. 原始信号x(t)（10秒，40960Hz采样）_______________________________________________________________________
# 1.1 生成信号 _________________________________________________________________________________________
x = m * np.cos(H * omega * t + np.radians(phi))
#px.line(x=t, y=x, labels={"x": "时间 (s)", "y": "幅度"}, title="原始信号x").show()

# 1.2 加窗 _______________________________________________________________________________________________
x_windowed = x*np.hanning(len(x))
#px.line(x=t, y=x_windowed, labels={"x": "时间 (s)", "y": "幅度"}, title="加窗后原始信号").show()

# 1.3 FFT分析 _____________________________________________________________________________________________________
x_windowed_fft = np.fft.fft(x_windowed)
# 1.3.1 FFT幅度谱 _________________________________________________________________________________________
x_windowed_fft_abs = np.abs(x_windowed_fft)
x_windowed_fft_abs_normalized = x_windowed_fft_abs/len(x)*2   # 因Hanning窗使信号幅度降低约一半，故乘以2归一化
x_windowed_fft_frequencies = np.fft.fftfreq(len(x), d=1/sampling_rate)
(px.scatter(x=x_windowed_fft_frequencies, y=x_windowed_fft_abs_normalized,
           labels={"x": "频率 (Hz)", "y": "幅度"}, title="原始信号加窗后FFT归一化幅度谱").
           update_xaxes(range=[30, 70]).update_yaxes(range=[0, fft_y_axis_limit]).show())


# 2. 重采样信号x_res(t)（10秒，819.2Hz采样）______________________________________________________________
# 2.1 重采样/生成信号 _____________________________________________________________________________
resample_indices = np.linspace(0, len(x)-1, fft_size_res, dtype=int)
x_res = x[resample_indices]
t_res = t[resample_indices]
sampling_rate_resampling = sampling_rate/(len(x)/len(x_res))  # 重采样后的等效采样率
#px.line(x=t_res, y=x_res, labels={"x": "时间 (s)", "y": "幅度"}, title="重采样信号x_res").show()

# 2.2 加窗 _______________________________________________________________________________________________
x_res_windowed = x_res*np.hanning(len(x_res))
#px.line(x=t_res, y=x_res_windowed, labels={"x": "时间 (s)", "y": "幅度"}, title="加窗后重采样信号").show()

# 2.3 FFT分析 _____________________________________________________________________________________________________
x_res_windowed_fft = np.fft.fft(x_res_windowed)
# 2.3.1 FFT幅度谱 _________________________________________________________________________________________
x_res_windowed_fft_abs = np.abs(x_res_windowed_fft)
x_res_windowed_fft_abs_normalized = x_res_windowed_fft_abs/len(x_res)*2   # 因Hanning窗使信号幅度降低约一半，故乘以2归一化
x_res_windowed_fft_frequencies = np.fft.fftfreq(len(x_res), d=1/sampling_rate_resampling)
(px.scatter(x=x_res_windowed_fft_frequencies, y=x_res_windowed_fft_abs_normalized,
           labels={"x": "频率 (Hz)", "y": "幅度"}, title="重采样信号加窗后FFT归一化幅度谱").
           update_xaxes(range=[30, 70]).update_yaxes(range=[0, fft_y_axis_limit]).show())

问题解答 #

1. 频谱泄漏的成因 #

重采样后的等效采样率为819.2Hz，此时50Hz信号在10秒时长内的采样周期数为 819.2*10/50=163.84——这不是整数，意味着重采样后的信号并非整周期截断。加Hanning窗后，窗函数的截断效应与非整周期信号叠加，触发了频谱泄漏。

同时，FFT点数设为8192时，频率分辨率为 819.2/8192=0.1Hz，泄漏旁瓣的频率（n*5Hz）正好是分辨率的整数倍，导致旁瓣能量集中落在FFT栅格点上，被显著观测到。

2. FFT样本数改为8193时泄漏消失的原因及实用性 #

当FFT点数改为8193时，频率分辨率变为 819.2/8193≈0.0999878Hz，此时泄漏旁瓣的频率不再与FFT栅格点对齐，旁瓣能量被分散到多个FFT bin中，幅度降低到无法被观测的水平，看起来像是泄漏消失。

这种方法实用性极低：一是仅适配当前特定参数组合，信号频率、采样率或时长变化后需重新调整FFT点数，无通用性；二是非2^n的FFT点数在C++实现中效率远低于基2FFT，会大幅增加实时处理的性能开销。

3. 保持FFT样本数为2^n时避免泄漏的方法 #

可以通过以下几种方式实现：

• 调整重采样率为信号频率的整数倍：比如将等效采样率改为800Hz（50*16），此时10秒内信号周期数为整数，整周期截断加窗后可避免泄漏。
• **整周期截断后补零到2^{n点数**：若无法调整采样率，先截取重采样信号的整周期片段（比如164个周期，时长3.28秒），再补零到最近的2}n点数（如4096）进行FFT。
• 使用低旁瓣窗函数：改用Flattop窗等旁瓣幅度极低的窗型，即使存在泄漏，其幅度也会被抑制到可接受水平，同时保持FFT点数为2^n。
• 重采样时采用插值方法：替换直接抽取的重采样方式，使用线性插值或sinc插值，减少重采样过程中的失真，降低后续FFT的泄漏程度。

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