FFT加窗处理出现异常多小峰问题求助

你的问题核心出在窗函数手动实现的索引错误，以及降采样未做抗混叠滤波两个关键点，导致加窗后频谱出现额外小峰：

1. 手动Hanning窗的索引错误 #

你用n=1:N计算窗函数，不符合Hanning窗的标准定义。标准Hanning窗针对0到N-1的采样点设计：
w(n) = 0.5*(1 - cos(2πn/(N-1)))，其中n∈[0, N-1]

当用n=1:N时，最后一个点代入会得到非零值，导致窗函数结尾无法归零，信号加窗后两端不能平滑过渡，引入额外频谱泄漏，表现为峰值间的小峰。

解决方法：直接用MATLAB内置的hann()函数生成窗，避免手动计算错误：

window = hann(N); % 或 hann(N, 'periodic')，按需选择周期/对称模式

如果坚持手动实现，必须用0-based索引：

n = 0:N-1; n = n';
hanning = 0.5*(1 - cos(2*n*pi/(N-1)));

2. 降采样未做抗混叠滤波 #

你直接对原信号抽点降采样，没有先做抗混叠滤波。原采样率44100Hz，降采样后为2756.25Hz，奈奎斯特频率1378.125Hz。若原弹跳球信号存在高于该频率的成分，会直接混叠到低频段，形成额外小峰。

解决方法：降采样前先对原信号应用低通滤波器，截止频率设为降采样后奈奎斯特频率的90%左右：

% 设计抗混叠低通滤波器
order = 8;
fSamp_raw = 44100;
downsample_ratio = 16;
fSamp = fSamp_raw / downsample_ratio;
cutoff_norm = 0.9*(fSamp/2)/(fSamp_raw/2); % 归一化截止频率
[b, a] = butter(order, cutoff_norm);
d_filtered = filtfilt(b, a, A(:,1));
y_filtered = filtfilt(b, a, A(:,2));

% 再执行降采样
d = d_filtered(1:downsample_ratio:2^16);
y = y_filtered(1:downsample_ratio:2^16);

3. 额外优化点 #

• 对FFT结果做归一化，方便对比加窗与不加窗的幅度：

  Pdd = (abs(fft(d))/N).^2;
  Pyy = (abs(fft(y))/N).^2;
  Pzw = (abs(fft(zw))/N).^2;
  

• 加窗后乘以能量补偿因子，修正幅度偏差：

  window_factor = sum(window.^2)/N;
  Pzw = Pzw / window_factor;
  

修正后的完整代码示例 #

% 加载数据（假设A已提前加载）
d_raw = A(:,1);   % 驱动信号
y_raw = A(:,2);   % 弹跳球信号

% 降采样前做抗混叠滤波
fSamp_raw = 44100;
downsample_ratio = 16;
fSamp = fSamp_raw / downsample_ratio;
nyquist_new = fSamp / 2;

% 设计巴特沃斯低通滤波器
order = 8;
cutoff_norm = 0.9 * nyquist_new / (fSamp_raw/2);
[b, a] = butter(order, cutoff_norm);
d_filtered = filtfilt(b, a, d_raw);
y_filtered = filtfilt(b, a, y_raw);

% 执行降采样
d = d_filtered(1:downsample_ratio:2^16); 
y = y_filtered(1:downsample_ratio:2^16);

N = length(y);
% 使用内置Hann窗
window = hann(N);
zw = y .* window;

% 计算归一化功率谱
Pdd = (abs(fft(d))/N).^2; 
Pyy = (abs(fft(y))/N).^2; 
Pzw = (abs(fft(zw))/N).^2;

% 补偿窗的能量损失
window_factor = sum(window.^2)/N;
Pzw = Pzw / window_factor;

% 提取正频率部分
stop = floor(N/2); 
Pdd = Pdd(1:stop);
Pyy = Pyy(1:stop);
Pzw = Pzw(1:stop);

% 生成频率轴
fdd = (1:stop) * (fSamp/2) / stop; 

% 绘图
semilogy(fdd,Pdd,fdd,Pyy,fdd,Pzw)
title('Unwindowed/Windowed')
xlabel('frequency (Hz)')
legend('Driving (Unwindowed)','Bouncing Ball (Unwindowed)','Bouncing Ball (Windowed)')

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