MATLAB实战：用频率抽样法设计FIR低通滤波器（附完整代码）

MATLAB实战用频率抽样法设计FIR低通滤波器附完整代码最近在做一个音频信号处理的小项目需要滤除一些高频噪声自然而然地就想到了FIR滤波器。和IIR滤波器比起来FIR没有反馈回路天生就是稳定的线性相位特性对保持信号波形也特别友好。不过设计FIR滤波器的方法有好几种窗函数法简单直接但过渡带和阻带衰减的权衡总让人头疼等波纹法Parks-McClellan算法性能最优可参数调起来又有点复杂。而频率抽样法在我看来是介于两者之间一个非常直观且灵活的选择——它允许你直接在频域上“描绘”出你想要的滤波器响应然后通过逆傅里叶变换得到时域系数整个过程就像在频域里“画图”一样清晰。尤其当你对滤波器的频率响应有非常具体、非标准的要求时这个方法就显得格外得心应手。这篇文章就是为你——无论是正在完成课程设计的学生还是需要快速实现一个定制化滤波器的工程师——准备的一份实战指南。我们将完全聚焦于MATLAB环境手把手地使用fir2函数从原理到代码从参数设置到效果调优完整地走一遍频率抽样法设计FIR低通滤波器的流程。我会分享一些调整过渡点的实用技巧并对比不同设计对最终滤波效果的影响所有代码都是可运行、可修改的希望能帮你把书本上的理论快速转化为解决实际问题的能力。1. 频率抽样法在频域“绘制”你的滤波器在深入代码之前我们有必要花点时间理解频率抽样法到底在做什么。它的核心思想异常直接既然数字滤波器的频率响应是周期性的那么我们能否直接在频域上在一个周期内通常是0到π或0到Fs/2指定一组离散频率点上的理想幅值然后通过数学方法反推出滤波器的时域冲激响应也就是滤波器系数呢频率抽样法给出的答案是肯定的。这种方法的美妙之处在于其直观性。你不需要像窗函数法那样先选窗、再算阶数也不像等波纹法那样去迭代逼近一个最优解。你就像个设计师在频率坐标纸上用点标出你期望的滤波器在各个频率点应该有的增益。比如对于低通滤波器在通带频率范围内你希望增益接近10 dB在阻带范围内你希望增益接近0负很多dB。这些你标出的点就是“频率抽样点”。注意这里有一个关键点我们指定的幅值通常是实数代表该频率点的幅度。但为了最终得到的滤波器系数是实数这是物理可实现的因果系统所要求的我们指定的频率响应必须满足共轭对称条件。好在对于实系数的滤波器其频率响应本身就具有共轭对称性。在使用fir2函数时我们只需要指定从0到1对应0到Fs/2的频率点及其幅值函数会自动帮我们处理对称部分。那么指定了这些离散点的幅值后如何得到连续的频率响应和时域系数呢这里就用到了离散傅里叶逆变换IDFT。我们指定的离散频率点幅值可以看作是对理想连续频率响应进行等间隔采样的结果。对这套采样值进行IDFT就能得到有限长的时域序列即我们的FIR滤波器系数。这个过程可以概括为以下几个步骤确定滤波器阶数N阶数决定了滤波器的长度系数个数为N1也直接影响频率分辨率。N越大你能在频域设置的抽样点就越密理论上越能逼近你想要的响应。在归一化频率轴上设置抽样点将0到采样频率一半Fs/2的频率范围映射到0到1的区间。在这个区间上等间隔或非均匀地选取M个点M通常与N有关。指定抽样点的理想幅值根据滤波器的类型低通、高通、带通等为每个抽样点赋予一个幅度值A。例如低通滤波器在通带内A1在阻带内A0。构造频域抽样序列并做IDFT将上一步得到的幅度序列结合必要的相位信息通常为了获得线性相位会赋予一个特定的线性相位构成一个完整的频域序列。对这个序列做N点的IDFT取实部就得到了FIR滤波器的N1个实系数。验证与调整计算实际滤波器的频率响应与理想响应对比。通常需要反复调整阶数N或在通带和阻带之间设置过渡点来优化性能比如减小通带波纹、增加阻带衰减。其中过渡点的设置是频率抽样法的一个精髓也是我们后面实战中要重点玩弄的“旋钮”。直接在通带和阻带之间进行0到1的跳变会导致时域系数序列产生剧烈的吉布斯振荡反映在频域就是通带和阻带内出现较大的起伏波纹。在跳变处插入一个或多个具有中间幅值比如0.4, 0.6的过渡点可以有效地平滑这种突变显著改善滤波器的幅频响应特性。当然这会以稍微加宽过渡带为代价。2. 实战准备理解MATLAB的fir2函数在MATLAB中频率抽样法的实现被封装成了一个非常便捷的函数fir2。它的存在让我们无需手动去计算IDFT和构造对称序列只需关心核心的设计意图。我们先来彻底搞懂它的调用方式。fir2函数的基本语法是b fir2(N, F, A)b: 输出的一维向量即设计得到的FIR滤波器系数单位冲激响应长度为N1。N: 滤波器的阶数。阶数越高滤波器越长通常性能越好过渡带更陡阻带衰减更大但计算量也越大。F: 一个从0到1递增的向量指定了归一化的频率抽样点。F1对应于奈奎斯特频率即采样频率Fs的一半。这个向量必须从0开始到1结束。A: 一个与F长度相同的向量指定了在F中各频率点处所期望的滤波器幅度响应。例如你想设计一个截止频率在0.4π归一化频率0.4的低通滤波器你可以这样定义F [0, 0.4, 0.41, 1]; % 频率点0, 0.4π, 0.41π, π A [1, 1, 0, 0]; % 对应幅度通带内为1从0.41π开始阻带为0这里在0.4和0.41之间设置了一个非常窄的过渡带实际上就是模拟了一个近似的理想矩形响应。fir2函数还有一个增强版本允许你指定窗函数来进一步优化系数b fir2(N, F, A, window)window: 一个长度为N1的窗函数向量。默认情况下fir2使用汉明窗Hamming。使用不同的窗函数可以调整滤波器的旁瓣特性。例如使用凯泽窗Kaiser可以更灵活地控制主瓣宽度和旁瓣衰减的权衡。为了更清晰地展示不同参数组合下的设计思路我们可以参考下表设计目标阶数N策略频率向量F设置要点幅度向量A设置要点可选窗函数window建议标准低/高通根据过渡带宽和阻带衰减需求估算初始值再微调。明确标出通带截止频率、过渡带起点/终点、阻带起点。通带内为1阻带内为0。过渡点可设中间值。默认汉明窗均衡性好。追求高衰减可用布莱克曼窗。需灵活权衡用凯泽窗。多带通/带阻由于响应复杂通常需要更高的阶数。需要定义多个通带和阻带区间频率点会更多。在不同频带区间内A在1和0之间切换。汉明或汉宁窗较为常用以减少带间波纹。任意形状响应取决于响应曲线的细节复杂度细节越多要求N越大。在响应变化剧烈的区域频率点需要设置得更密集。直接描绘出你想要的幅度曲线值。通常先尝试默认窗若纹波过大再尝试其他窗进行平滑。提示fir2函数内部会自动将你指定的(F, A)对进行线性插值形成一个连续的期望频率响应然后进行频率抽样和加窗设计。因此你不需要设置非常密集的频率点只需在响应发生变化的关键频率点如截止频率、过渡带边界进行指定即可。3. 从零开始一个低通滤波器设计实例现在让我们设定一个具体的滤波器指标并从头开始用fir2实现它。假设我们的需求如下采样频率 Fs: 4000 Hz通带边界频率 Fpass: 800 Hz阻带边界频率 Fstop: 1000 Hz通带最大波纹: 0.5 dB (近似认为通带内增益在 -0.5dB ~ 0dB 之间)阻带最小衰减: 40 dB (即阻带内增益 -40dB)首先我们需要将实际频率Hz转换为fir2函数所需的归一化频率0到1之间。归一化频率 实际频率 / (Fs/2)。通带边界归一化频率:Fpass_norm 800 / (4000/2) 0.4阻带边界归一化频率:Fstop_norm 1000 / (4000/2) 0.5一个最直接的想法是在0.4处从通带跳到阻带。我们尝试设计一个阶数N20的滤波器。clc; clear; close all; % 设计参数 Fs 4000; % 采样频率 (Hz) Fpass 800; % 通带边界 (Hz) Fstop 1000; % 阻带边界 (Hz) N 20; % 滤波器阶数 % 计算归一化频率 F_norm [0, Fpass/(Fs/2), Fstop/(Fs/2), 1]; % 频率向量: 0, 0.4, 0.5, 1 A [1, 1, 0, 0]; % 幅度向量: 通带1阻带0 % 使用fir2设计滤波器 b_direct fir2(N, F_norm, A); % 分析频率响应 [H_direct, f_direct] freqz(b_direct, 1, 1024, Fs); mag_direct abs(H_direct); Hdb_direct 20*log10(mag_direct); % 绘图 figure(Position, [100, 100, 800, 400]); plot(f_direct, Hdb_direct, b-, LineWidth, 1.5); grid on; xlabel(频率 (Hz)); ylabel(幅度 (dB)); title(直接跳变设计 (N20) 的幅频响应); xline(Fpass, r--, 通带边界, LabelVerticalAlignment, middle); xline(Fstop, g--, 阻带边界, LabelVerticalAlignment, middle); yline(-0.5, k:, 通带波纹 0.5dB); yline(-40, k:, 阻带衰减 40dB); legend(滤波器响应, Location, best); axis([0, Fs/2, -80, 5]);运行这段代码你会看到滤波器的响应曲线。很大概率上这个简单的设计无法满足我们的要求通带波纹可能超过0.5dB阻带衰减也可能远未达到40dB。问题就出在从0.4到0.5的幅度从1直接跳变到0引发了严重的吉布斯现象。4. 性能优化的关键过渡点魔法为了改善性能我们必须引入过渡点。过渡点是在通带和阻带之间的频率点上赋予一个介于1和0之间的幅度值。这相当于告诉滤波器“不要在这里突然变化平滑地过渡过去。” 这能有效抑制通带和阻带内的波纹振荡。4.1 插入单个过渡点我们可以在通带边界0.4和阻带边界0.5之间插入一个过渡频率点比如0.45并给它一个合适的幅度例如0.6这是一个经验值可以调整。% 插入一个过渡点 F_norm_1trans [0, 0.4, 0.45, 0.5, 1]; % 频率点 A_1trans [1, 1, 0.6, 0, 0]; % 对应幅度 b_1trans fir2(N, F_norm_1trans, A_1trans); [H_1trans, f] freqz(b_1trans, 1, 1024, Fs); Hdb_1trans 20*log10(abs(H_1trans));4.2 插入两个过渡点更进一步我们可以插入两个过渡点提供更平滑的过渡。例如在0.42和0.48处设置幅度0.8和0.2。% 插入两个过渡点 F_norm_2trans [0, 0.4, 0.42, 0.48, 0.5, 1]; A_2trans [1, 1, 0.8, 0.2, 0, 0]; b_2trans fir2(N, F_norm_2trans, A_2trans); [H_2trans, f] freqz(b_2trans, 1, 1024, Fs); Hdb_2trans 20*log10(abs(H_2trans));现在让我们将三种设计无过渡点、1个过渡点、2个过渡点的结果绘制在一起进行对比。figure(Position, [100, 100, 900, 500]); subplot(2,1,1); plot(f, Hdb_direct, k-, LineWidth, 1.2); hold on; plot(f, Hdb_1trans, r--, LineWidth, 1.5); plot(f, Hdb_2trans, b:, LineWidth, 1.8); hold off; grid on; xlabel(频率 (Hz)); ylabel(幅度 (dB)); title(不同过渡点设计对幅频响应的影响 (N20)); xline(Fpass, k--); xline(Fstop, k--); yline(-0.5, k:); yline(-40, k:); legend(无过渡点, 1个过渡点 (0.450.6), 2个过渡点 (0.420.8, 0.480.2), Location, southwest); axis([0, Fs/2, -80, 5]); % 局部放大观察过渡带和阻带 subplot(2,1,2); plot(f, Hdb_direct, k-, LineWidth, 1.2); hold on; plot(f, Hdb_1trans, r--, LineWidth, 1.5); plot(f, Hdb_2trans, b:, LineWidth, 1.8); hold off; grid on; xlabel(频率 (Hz)); ylabel(幅度 (dB)); title(局部放大图 (重点关注过渡带和阻带)); xline(Fpass, k--); xline(Fstop, k--); yline(-40, k:); legend(无过渡点, 1个过渡点, 2个过渡点, Location, southwest); xlim([600, 1200]); % 放大到通带和阻带附近 ylim([-70, 2]);通过对比图你可以清晰地观察到无过渡点过渡带最陡峭几乎是垂直下降但通带和阻带内波纹振荡非常剧烈通带波纹远超0.5dB阻带衰减在很多频点也达不到40dB要求。这是吉布斯现象的典型表现。1个过渡点过渡带明显变宽了但通带和阻带内的波纹被显著压制。通带更加平坦阻带衰减整体加深更容易满足40dB的要求。2个过渡点过渡带进一步加宽但阻带衰减性能通常是最好的通带也更为平滑。它用更宽的过渡带换来了更优的带内特性。提示过渡点的位置和幅度值是需要反复调试的“艺术”。没有绝对的最优解只有针对你特定指标通带波纹、阻带衰减、过渡带宽的权衡。一个常见的策略是先根据经验公式或简单估算设置初始值然后通过循环微调并观察响应曲线找到满足所有约束的最小阶数N和过渡点配置。4.3 结合窗函数微调除了调整过渡点我们还可以在fir2中更换窗函数。默认的汉明窗在旁瓣抑制和主瓣宽度之间取得了很好的平衡。如果你需要更深的阻带衰减可以尝试布莱克曼窗blackman(N1)它的旁瓣更小但主瓣更宽意味着过渡带会更缓。如果你需要精确控制波纹凯泽窗kaiser(N1, beta)提供了参数beta来调节这种权衡。% 使用布莱克曼窗设计两个过渡点案例 window_blackman blackman(N1); b_2trans_blackman fir2(N, F_norm_2trans, A_2trans, window_blackman); [H_2trans_b, f] freqz(b_2trans_blackman, 1, 1024, Fs); Hdb_2trans_b 20*log10(abs(H_2trans_b)); % 对比汉明窗默认和布莱克曼窗 figure; plot(f, Hdb_2trans, b-, LineWidth, 1.5); hold on; plot(f, Hdb_2trans_b, r--, LineWidth, 1.5); hold off; grid on; xlabel(频率 (Hz)); ylabel(幅度 (dB)); title(不同窗函数对比 (2个过渡点N20)); legend(默认汉明窗, 布莱克曼窗, Location, best); axis([0, Fs/2, -100, 5]); % 扩大纵轴范围看阻带衰减你会发现使用布莱克曼窗后阻带衰减在大部分频段可能更深但过渡带的边缘可能会变得更“软”且通带波纹特性也可能发生变化。这再次体现了滤波器设计中的权衡。5. 完整设计流程与代码整合综合以上所有技巧一个稳健的频率抽样法FIR滤波器设计流程可以归纳如下明确指标确定Fs, Fpass, Fstop, 通带波纹(Apass), 阻带衰减(Astop)。估算初始阶数N可以根据过渡带宽Δf (Fstop - Fpass)/Fs进行粗略估算。一个经验起点是N ≈ 4 / Δf。例如本例中 Δf 200/40000.05 N≈80。我们之前用N20是为了演示过渡点效果实际为满足40dB衰减可能需要更高的阶数。设置归一化频率向量F和幅度向量A总是以F [0, ... , 1]和A [..., ...]的格式。在通带和阻带边界之间插入1到3个过渡点。过渡点幅度可以尝试等间隔值如0.66, 0.33或根据优化算法确定。首次设计使用b fir2(N, F, A)进行设计。性能验证使用freqz计算频率响应检查通带最大波纹和阻带最小衰减是否达标。计算通带波纹max(20*log10(abs(H_pass))) - min(20*log10(abs(H_pass)))计算阻带衰减-max(20*log10(abs(H_stop)))(因为阻带内响应小于1对数为负值取负后得到正衰减值)迭代优化如果阻带衰减不足尝试a) 增加过渡点b) 优化过渡点幅度c) 增加滤波器阶数N。如果通带波纹过大尝试a) 调整过渡点位置和幅度b) 更换窗函数如尝试凯泽窗并调整beta值c) 增加阶数N。如果过渡带过宽而带内性能已达标可以尝试a) 减少过渡点b) 微调过渡点位置使其更靠近边界c) 在满足性能前提下尝试降低阶数N。最终实现与测试得到满意的系数向量b后可以使用filter(b, 1, x)对信号x进行滤波。下面是一个整合了自动性能检查的、更完整的示例代码框架它尝试通过微调过渡点幅度来逼近指标clc; clear; close all; % 1. 设计指标 Fs 4000; Fpass 800; Fstop 1000; Apass_dB 0.5; % 通带最大波纹 (dB) Astop_dB 40; % 阻带最小衰减 (dB) % 2. 初始参数设置 (可调整) N 60; % 从经验公式估算的更高阶数 F_norm [0, Fpass/(Fs/2), 0.45, 0.5, 1]; % 一个过渡点 % 过渡点幅度初始猜测值 trans_amp 0.5; % 3. 优化循环示例 (这里简化仅微调过渡点幅度) best_b []; best_metrics [inf, 0]; % [通带波纹, 阻带衰减] for amp 0.3:0.05:0.7 A [1, 1, amp, 0, 0]; b fir2(N, F_norm, A); [H, f] freqz(b, 1, 4096, Fs); H_db 20*log10(abs(H)); % 提取通带和阻带响应 idx_pass (f 0) (f Fpass); idx_stop (f Fstop) (f Fs/2); Hpass_db H_db(idx_pass); Hstop_db H_db(idx_stop); % 计算性能指标 ripple_pass max(Hpass_db) - min(Hpass_db); % 通带波纹(dB) atten_stop -max(Hstop_db); % 阻带衰减(dB) % 判断是否优于当前最佳且满足衰减要求 if atten_stop Astop_dB ripple_pass best_metrics(1) best_metrics [ripple_pass, atten_stop]; best_b b; best_amp amp; fprintf(发现更优解: 过渡点幅度%.2f, 通带波纹%.2f dB, 阻带衰减%.2f dB\n, ... amp, ripple_pass, atten_stop); end end % 4. 展示最终设计结果 if ~isempty(best_b) figure(Position, [100, 100, 800, 600]); [H_final, f_final] freqz(best_b, 1, 4096, Fs); Hdb_final 20*log10(abs(H_final)); plot(f_final, Hdb_final, b-, LineWidth, 1.8); grid on; hold on; xline(Fpass, r--, F_{pass}800Hz, LineWidth, 1.2); xline(Fstop, g--, F_{stop}1000Hz, LineWidth, 1.2); yline(-Astop_dB, k:, sprintf(阻带衰减要求: %d dB, Astop_dB), LineWidth, 1.2); % 绘制通带波纹范围 yline(0, k-.); yline(-Apass_dB, k-.); patch([0, Fpass, Fpass, 0], [-Apass_dB, -Apass_dB, 0, 0], r, FaceAlpha, 0.1, EdgeColor, none); hold off; xlabel(频率 (Hz)); ylabel(幅度 (dB)); title(sprintf(最终设计: N%d, 过渡点幅度%.2f\\n通带波纹≈%.2f dB, 阻带衰减≈%.2f dB, ... N, best_amp, best_metrics(1), best_metrics(2))); legend(滤波器响应, Location, southwest); axis([0, Fs/2, -80, 2]); % 显示滤波器系数 disp(滤波器系数 (前10个):); disp(best_b(1:10)); else disp(未找到满足阻带衰减要求的设计请尝试增加阶数N或增加过渡点。); end这个框架提供了一个自动搜索过渡点幅度的简单示例。在实际更复杂的设计中你可能需要将阶数N、过渡点数量和位置、窗函数类型等都纳入优化变量甚至使用更高级的优化算法如fminsearch来寻找满足多重指标约束的最优解。不过对于大多数工程应用通过手动调整和观察响应曲线已经能够快速得到一个可用的优质滤波器。