MATLAB FIR陷波滤波器代码报错求助：无法滤除音频特定频率

问题总结 #

编写MATLAB代码为SunshineSquare.wav设计FIR陷波滤波器滤除特定干扰频率，但初始代码和修改后代码均存在错误导致运行失败。

初始代码错误分析 #

[x, fs] = audioread('SunshineSquare.wav'); 
 
X = fft(x); 
f = linspace(0, fs/2, length(X)/2);
[~, f0] = max(abs(X(1:length(X)/2))); 
 
notchFilter = fdesign.notch('N,Fc', 100, f0, fs); 
b = firpm(notchFilter); 

 
y = filter(b, 1, x); 


sound(x, fs);

• 干扰频率计算错误：max返回的f0是频谱的索引值，不是实际频率，直接传给fdesign.notch的Fc参数（需要实际频率或归一化频率）会导致滤波器参数完全错误。
• 滤波器设计函数使用错误：firpm无法直接接收fdesign.notch对象作为输入，正确做法是用design函数基于fdesign对象生成滤波器。
• 输出错误：最后播放的是原始信号x，而非滤波后的y，无法验证滤波效果。

修改后代码错误分析 #

[x, fs] = audioread('SunshineSquare.wav'); 
x = y; fs = Fs;
X = fft(x); 
L = floor(length(X)/2);
f = linspace(0, fs/2, L);
[~, f0idx] = max(abs(X(1:L))); 
f0 = (f0idx-1)/L;
Q = 2.5;
notchFilter1 = fdesign.notch('N,F0,Q', 6, f0, Q); 
notchFilter2 = fdesign.notch('N,F0,Q', 6, f0*2, Q);
Hd1 = design(notchFilter1, 'SystemObject', true);
Hd2 = design(notchFilter2, 'SystemObject', true);
y = Hd2(Hd1(x));
sound(y, fs)
tiledlayout('flow');
nexttile
plot(f, abs(X(1:L))); title('original spectrum');
nexttile
Y = fft(y);
plot(f, abs(Y(1:L))); title('filtered spectrum');

• 未定义变量赋值：第二行x = y; fs = Fs;中，y和Fs均未定义，直接触发运行错误。
• 频率归一化错误：f0 = (f0idx-1)/L得到的是[0,1)范围的相对值，但fdesign.notch的F0参数若使用实际频率，需要结合采样率转换；若用归一化频率，应对应到[0,1]（代表fs/2），此处计算未关联采样率，导致滤波器目标频率错误。
• 滤波器阶数与设计方法缺失：仅指定阶数N=6但未明确FIR设计方法，可能导致滤波器性能不达标。

修正后的完整代码 #

% 读取音频文件
[x, fs] = audioread('SunshineSquare.wav'); 
% 处理立体声（如果是双声道，转为单声道）
if size(x,2) > 1
    x = mean(x, 2);
end

% 计算原始信号频谱，定位干扰频率
X = fft(x);
L = floor(length(X)/2);
f = linspace(0, fs/2, L);
% 获取频谱幅值
X_amp = abs(X(1:L));
% 找到幅值最大的频率点（干扰频率）
[~, f0idx] = max(X_amp);
% 转换为实际频率值
f0 = f(f0idx);
% 倍频干扰频率
f0_2 = 2*f0;

% 设计FIR陷波滤波器：指定阶数、中心频率、品质因数，使用firpm方法
Q = 2.5;
% 第一个陷波滤波器（针对f0）
notchSpec1 = fdesign.notch('N,F0,Q', 100, f0, Q, fs);
Hd1 = design(notchSpec1, 'firpm', 'SystemObject', true);
% 第二个陷波滤波器（针对2*f0）
notchSpec2 = fdesign.notch('N,F0,Q', 100, f0_2, Q, fs);
Hd2 = design(notchSpec2, 'firpm', 'SystemObject', true);

% 执行滤波
y = Hd2(Hd1(x));

% 播放滤波后的音频
sound(y, fs);

% 绘制频谱对比图
tiledlayout('flow');
nexttile
plot(f, X_amp); 
title('原始信号频谱');
xlabel('频率 (Hz)'); ylabel('幅值');
nexttile
Y = fft(y);
Y_amp = abs(Y(1:L));
plot(f, Y_amp); 
title('滤波后信号频谱');
xlabel('频率 (Hz)'); ylabel('幅值');

关键修正点说明 #

• 增加双声道处理：若音频为立体声，转为单声道避免滤波异常。
• 正确计算干扰频率：通过f(f0idx)直接获取实际频率值，确保滤波器目标频率准确。
• 明确滤波器设计方法：指定'firpm'方法设计FIR陷波滤波器，同时提高阶数到100以获得更陡峭的陷波特性。
• 删除未定义变量赋值：移除导致直接报错的x = y; fs = Fs;语句。
• 完善频谱图标注：添加坐标轴标签，提升可读性。

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