很多人刚开始用Matlab做数字信号处理时最容易碰到的坎就是“低通滤波器到底用哪个函数”。butter、cheby1、cheby2、ellip、fir1、firpm这一排函数摆在面前每个都能设计出低通滤波器但出来的结果差异比想象中大得多。这篇博文我就用一组统一指标把这些函数逐个拉出来实测对比把阶数、幅频响应、相位特性、群延迟这些关键性能讲透顺便整理我在实际调试中踩过的坑。不管你是课程设计要交滤波器实验报告还是做嵌入式项目需要选一个实时滤波方案或者纯粹想把Matlab的滤波器设计工具箱用明白这篇文章都值得你花几分钟看完。1. 先把需求翻译成指标一切对比的前提1.1 通带纹波与阻带衰减被忽视的“硬指标”很多人设计低通滤波器时只知道说“我要一个截止频率100Hz的低通”到了Matlab里就写一句[b,a] butter(4, 0.2)觉得大功告成。但实际工程里光有截止频率远远不够。一个专业的滤波器设计需求通常包含四个关键参数通带边界频率、阻带起始频率、通带最大纹波、阻带最小衰减。通带纹波指的是通带内幅度响应允许的最大起伏单位为dB阻带衰减则要求阻带内的信号被压低到多少dB以下。比如心电信号预处理常常要求通带纹波小于1dB、阻带衰减大于60dB这个60dB就意味着阻带信号能量被抑制到千分之一以下。为什么不直接说一个“截止频率”就完事因为模拟理想低通在频域上是矩形的但物理可实现的所有滤波器都有一个从通带到阻带的过渡过程。过渡带的宽度、通带是否平坦、阻带压得够不够深这些才是真正决定滤波器性能的维度。把需求翻译成指标才有办法横向对比不同设计函数。1.2 过渡带宽度与滤波器阶数计算资源的第一道门槛过渡带就是从通带边缘到阻带起始之间的那段频率区间。过渡带越窄滤波器的选择性越好但代价是阶数必须提高。阶数越高意味着计算量越大对于实时系统来说就是延迟和CPU负担。这里有个容易混淆的点IIR滤波器的阶数和FIR滤波器的阶数物理意义完全不同。IIR的4阶意味着系统函数分母是4阶多项式对应4个极点而FIR的4阶严格说是4个系数实际阶数N-1可能连一个像样的过渡带都做不出来。所以后面做对比的时候不能光看“谁阶数低谁厉害”要看同等的指标约束下谁用的资源更少。1.3 相位响应与群延迟波形失真的隐形杀手幅度响应只告诉你“哪些频率被放大了哪些被衰减了”相位响应则告诉你“不同频率的信号经过滤波器后延迟了多久”。如果相位响应不是线性的信号经过滤波器后波形就会畸变——这在时域里表现为过冲、拖尾、波形不对称。群延迟就是相位响应对频率的导数单位是采样点或秒。如果群延迟在通带内是个常数说明各频率成分延迟一致波形形状能保持如果群延迟随频率剧烈变化波形就会乱。IIR滤波器天然是非线性相位的越靠近截止频率群延迟变化越剧烈FIR滤波器如果系数对称可以实现严格的线性相位群延迟是恒定值但代价是你得付出高得多的阶数。1.4 本次对比的测试基准为了公平起见我给所有设计函数设置了同一组指标采样率1000Hz通带截止频率100Hz阻带起始频率150Hz通带最大纹波1dB阻带最小衰减60dB。这个指标组合过渡带只有50Hz阻带要求60dB在工程里算是比较典型的中等强度需求既能体现各种滤波器的特性差异又不至于极端到脱离实际。后续所有对比都围绕这组指标展开。2. Matlab里6个常用滤波器设计函数逐个拆解2.1 butter与cheby1最常用的两支IIRbutter函数设计的是巴特沃斯滤波器它的最大特点是通带内幅度响应最平坦没有纹波从零频一直平缓下降到截止频率。这种特性让它在音频、传感器信号等对通带平坦度要求高的场景很受欢迎。代价是过渡带相对较宽达到同样阻带衰减需要的阶数往往比切比雪夫和椭圆滤波器高。% 设计巴特沃斯低通滤波器 Fs 1000; % 采样率 Fpass 100; % 通带截止频率 Fstop 150; % 阻带起始频率 Apass 1; % 通带纹波(dB) Astop 60; % 阻带衰减(dB) % 方式一根据阶数设计 [b, a] butter(8, 100/(Fs/2), low); % 方式二根据指标自动估算最低阶数 [N, Wn] buttord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b, a] butter(N, Wn, low);cheby1设计的是切比雪夫I型滤波器特点是通带内是等纹波的允许在1dB范围内上下波动但过渡带比同阶的巴特沃斯更陡峭。如果你对通带内的小幅度波动不敏感但希望用更低的阶数达到同样的阻带衰减切比雪夫I型是很好的选择。% 切比雪夫I型 [N, Wn] cheb1ord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b, a] cheby1(N, Apass, Wn, low);这里有个细节cheby1设计时通带纹波出现在通带内阻带是单调下降的。所以如果你的信号本身有微弱的有用成分恰好落在通带边缘纹波会导致这些频率被放大或衰减的程度不一致选型时要留意这一点。2.2 cheby2与ellip纹波位置的取舍cheby2是切比雪夫II型特点是通带内平坦无纹波纹波出现在阻带。这正好和I型相反。它的适用场景很明确你特别在意通带内的信号保真度不希望有用信号有任何幅度抖动但阻带里的泄漏噪声多一点少一点无所谓。不过切比雪夫II型的过渡带通常比I型更宽一些在同样指标下阶数会稍高。% 切比雪夫II型 [N, Wn] cheb2ord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b, a] cheby2(N, Astop, Wn, low);ellip函数设计椭圆滤波器也叫考尔滤波器它在通带和阻带内部都有纹波但换来的是同等指标下所有滤波器中最小的阶数。椭圆滤波器的过渡带最窄选择性最高但相位响应也是所有类型里最扭曲的。% 椭圆滤波器 [N, Wn] ellipord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b, a] ellip(N, Apass, Astop, Wn, low);实际项目中椭圆滤波器常常用来做“指标苛刻但算力有限”的场合。比如某款嵌入式设备需要70dB阻带衰减处理器资源又很紧张椭圆滤波器可能是唯一能在给定阶数内满足要求的方案。前提是你能接受通带纹波和较差的相位特性。2.3 fir1与firpm从窗函数到最优等波纹fir1是基于窗函数法的FIR设计函数原理是先计算理想低通的无限长脉冲响应再用窗函数截断。因为窗函数的存在最终的频率响应是理想响应和窗函数频谱的卷积所以会出现吉布斯现象——过渡带两侧产生纹波。fir1的优点是简单直观、结果稳定配合kaiser窗可以在阻带衰减和过渡带宽度之间做平衡。% fir1配合kaiser窗 fc (Fpass Fstop) / 2 / (Fs/2); % 截止频率取通带和阻带中点 % 根据阻带衰减估算kaiser窗的beta A Astop; if A 50 beta 0.1102*(A - 8.7); elseif A 21 beta 0.5842*(A - 21)^0.4 0.07886*(A - 21); else beta 0; end N 80; % 阶数需要试凑或估算 h fir1(N, fc, kaiser(N1, beta));firpm实现的是Parks-McClellan最优等波纹算法它把误差均匀地分布在通带和阻带内而不是像窗函数法那样误差集中在截止频率附近。这意味着在相同阶数下firpm能达到更小的最大误差或者说达到同样指标需要的阶数比fir1更低。% firpm设计与指标评估 f [0 Fpass Fstop Fs/2] / (Fs/2); % 频带边界 a [1 1 0 0]; % 期望幅度 dev [(10^(Apass/20)-1)/(10^(Apass/20)1), 10^(-Astop/20)]; % 通带和阻带纹波 [n, fo, ao, w] firpmord(f, a, dev, Fs); h_firpm firpm(n, fo, ao, w);用firpmord估算所需阶数时要注意它返回的n是最小阶数估计值实际设计后必须用freqz验证是否真的满足全部指标。我遇到过估算阶数偏小的情况设计出来的滤波器阻带衰减只有55dB差一点点不达标把阶数加两位数就好了。2.4 designfilt把参数统一起来的入口从某个版本开始Matlab提供了designfilt这个统一接口可以把上述所有设计方法封装成一套输入格式。好处是参数名称直观不容易搞混Wn的归一化方式而且可以直接生成一个digitalFilter对象后续分析都用对象方法。% 用designfilt设计巴特沃斯低通 d designfilt(lowpassiir, ... FilterOrder, 8, ... HalfPowerFrequency, 100, ... SampleRate, 1000); % 用designfilt设计等波纹FIR d2 designfilt(lowpassfir, ... PassbandFrequency, 100, ... StopbandFrequency, 150, ... PassbandRipple, 1, ... StopbandAttenuation, 60, ... SampleRate, 1000, ... DesignMethod, equiripple);不过designfilt也不是万能神器它内部其实还是调用上面那些底层函数而且有时候参数组合过严会直接报错。我实际用下来的体会是初期学习和调试用designfilt很友好但到了要精细控制设计细节、或者要把IIR转成零极点形式做数值优化时还是回到butter、cheby1这些底层函数更顺手。3. 同一组指标下的实测对比阶数、响应与相位差异3.1 最小阶数对比IIR和FIR不是同一个量级我按照前面的指标Fs1000HzFpass100HzFstop150HzApass1dBAstop60dB分别用各种函数设计了滤波器记录了达到指标所需的最小阶数。这里的“最小阶数”是用各函数自带的阶数估算函数如buttord、cheb1ord自动算出来的结果如下设计函数滤波器类型满足指标的最小阶数过渡带宽度butterIIR巴特沃斯约18阶较宽cheby1IIR切比雪夫I型约9阶较窄cheby2IIR切比雪夫II型约9阶中等ellipIIR椭圆约5阶最窄fir1kaiser窗FIR窗函数法约70阶中等firpmFIR最优等波纹约55阶较窄这组数据充分说明了IIR和FIR的本质差异IIR靠反馈结构用较少的阶数实现了陡峭的过渡带而FIR没有反馈只能靠增加阶数来逼近理想响应。5阶的椭圆滤波器能达到60dB阻带衰减这在FIR世界里几乎是不可能的任务。3.2 幅频响应、相位响应和群延迟的具体表现做完阶数对比我画了所有滤波器的幅频响应、相位响应和群延迟曲线几个关键现象值得展开说说。从幅频响应看巴特沃斯滤波器的通带最平坦从零频到接近截止频率几乎是一条直线但过渡带明显偏宽切比雪夫I型通带内有等幅度的波纹过渡带比巴特沃斯陡椭圆滤波器的通带和阻带都有波纹但过渡带最窄幅频响应曲线看起来最“矩形”。相位响应方面三个IIR滤波器的相位曲线在低频段都接近线性到截止频率附近开始明显弯曲椭圆滤波器弯曲得最厉害。FIR滤波器因为系数对称相位响应在整个通带内都是严格线性的这也是FIR在时域波形保真上最大的优势。群延迟曲线更能说明问题IIR滤波器的群延迟在通带内不是常数越靠近截止频率延迟越大最大群延迟往往出现在截止频率附近。椭圆滤波器的群延迟波动尤其剧烈通带边缘的延迟可能比低频段高出好几倍。这意味着如果一个信号包含100Hz和50Hz两个频率成分经过椭圆滤波器后两个成分的延迟时间会差很多叠加后波形自然就变形了。FIR滤波器的群延迟则是完美的常数等于(N-1)/2个采样周期所有频率成分同等延迟波形形状完全不变。3.3 一段脚本完成全部横向对比为了让你自己也能复现这个对比我给出完整脚本。这段脚本会生成对比图分别展示幅频响应和群延迟曲线一眼就能看出各滤波器的差异。% 低通滤波器设计函数性能对比 Fs 1000; Fpass 100; Fstop 150; Apass 1; Astop 60; % IIR滤波器设计使用最小阶数 [Nb, Wnb] buttord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b1, a1] butter(Nb, Wnb, low); [Nc1, Wnc1] cheb1ord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b2, a2] cheby1(Nc1, Apass, Wnc1, low); [Nc2, Wnc2] cheb2ord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b3, a3] cheby2(Nc2, Astop, Wnc2, low); [Ne, Wne] ellipord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b4, a4] ellip(Ne, Apass, Astop, Wne, low); % FIR滤波器设计 fc (Fpass Fstop) / 2 / (Fs/2); Nf 72; beta kaiserbeta(Astop); h1 fir1(Nf, fc, kaiser(Nf1, beta)); dev [(10^(Apass/20)-1)/(10^(Apass/20)1), 10^(-Astop/20)]; f [0 Fpass Fstop Fs/2] / (Fs/2); a [1 1 0 0]; [nf, fo, ao, w] firpmord(f, a, dev, Fs); h2 firpm(nf, fo, ao, w); % 幅频响应对比 figure; freqz(b1, a1, 1024, Fs); hold on; freqz(b2, a2, 1024, Fs); freqz(b3, a3, 1024, Fs); freqz(b4, a4, 1024, Fs); freqz(h1, 1, 1024, Fs); freqz(h2, 1, 1024, Fs); legend(butter, cheby1, cheby2, ellip, fir1-kaiser, firpm); grid on; % 群延迟对比 figure; grpdelay(b1, a1, 1024, Fs); hold on; grpdelay(b2, a2, 1024, Fs); grpdelay(b3, a3, 1024, Fs); grpdelay(b4, a4, 1024, Fs); grpdelay(h1, 1, 1024, Fs); grpdelay(h2, 1, 1024, Fs); legend(butter, cheby1, cheby2, ellip, fir1-kaiser, firpm); grid on;3.4 怎么读这些图别只盯着阻带衰减很多初学者看滤波器对比只关心阻带衰减够不够深、过渡带窄不窄然后草率下结论“椭圆滤波器最好”。但如果你的应用是时域波形分析比如心电信号、振动信号的特征提取相位失真比阻带衰减更致命。我不止一次见过有人用椭圆滤波器做完预处理后发现波形的峰值位置偏移了、形态变得很奇怪找半天原因才发现是相位响应的锅。所以在读对比图时先把幅频响应看个大概确认所有方案都满足指标要求然后重点看群延迟曲线在通带内的平坦程度。如果群延迟在通带内有剧烈波动说明这个滤波器对时域波形的影响很大。对于离线数据处理可以用filtfilt做零相位滤波来补救对于实时系统只能老老实实选FIR或者接受IIR带来的波形失真。4. IIR与FIR阵营的本质差异相位、稳定性和实时性怎么权衡4.1 反馈结构带来的高效与隐患IIR滤波器的系统函数里有分母多项式意味着输出不仅取决于当前和过去的输入还取决于过去的输出。这个反馈结构让IIR能用很少的阶数实现很陡的过渡带计算效率极高。同样是60dB阻带衰减、50Hz过渡带IIR只需要5~18阶FIR却要50~80阶在实时系统里这个差距直接决定了能不能跑得动。但反馈结构也带来两个隐患。第一个是稳定性问题。IIR的极点必须全部落在单位圆内否则滤波器会发散。设计阶段Matlab会保证极点位置合理但你用自己的代码处理系数时要小心比如直接对传输函数系数做浮点运算后可能改变极点位置。第二个是系数灵敏度问题。高阶IIR直接型结构对系数量化非常敏感尤其阶数超过10以后很小的量化误差就可能导致频率响应严重变形甚至不稳定。这也是为什么工程上高阶IIR都要用二阶节SOS级联实现而不是直接使用传输函数形式存系数。4.2 线性相位到底有多重要什么时候必须用FIR线性相位意味着所有频率成分经过滤波器的延迟完全相同波形不会发生相位失真。在以下场景里线性相位是刚需生物医学信号处理比如心电图、脑电图分析。医生看的是波形形态如果P波、QRS波群的相对位置被相位失真改变可能导致误判。通信系统中的匹配滤波器为了保证符号判决的最佳信噪比通常也要求线性相位。数据采集系统中的抗混叠滤波如果后面要重建原始波形相位失真同样不可接受。FIR滤波器只要系数满足对称或反对称条件就天然具有严格线性相位特性。代价是阶数高、延迟大、计算量大。FIR的延迟等于(N-1)/2个采样点阶数70的滤波器意味着输出比输入延迟了35个采样点在1000Hz采样率下就是35毫秒。对于实时控制系统35毫秒的延迟可能是致命的。4.3 实时系统里的选型逻辑如果是实时嵌入式系统算力有限、中断周期固定我一般按下面的逻辑选型先看是否允许相位失真。不允许就直接FIR然后在FIR里选firpm尽可能降低阶数。允许相位失真就优先考虑IIR然后在IIR里优先选椭圆滤波器因为它阶数最低计算量最小。如果阻带衰减要求特别高比如超过80dB椭圆滤波器容易出现数值问题这时候切比雪夫I型或巴特沃斯配合更高阶数反而更稳。还要考虑处理器的定点能力。IIR的反馈结构在定点实现时更容易受舍入误差影响可能出现极限环振荡。FIR没有反馈定点实现几乎不会出现稳定性问题。如果目标平台是低端MCU而且只能用定点数FIR反而更容易调通。我曾经在一个项目里用某款低端MCU做实时信号处理最开始选了6阶椭圆滤波器各项指标都完美但定点移植后输出总在零输入时出现小幅振荡排查了很久才锁定是IIR反馈路径的舍入误差问题。后来把设计改成等波纹FIR虽然阶数高了不少但定点实现极其稳定问题当场消失。4.4 filtfilt离线处理中的“作弊器”如果你做的是离线数据分析而不是实时滤波IIR相位失真带来的问题有一个非常实用的解决办法filtfilt双向零相位滤波。它的思路是把信号正向滤波一次再把滤波结果反向滤波一次两次滤波的相位失真方向相反正好相互抵消输出信号的相位响应为零。% 使用filtfilt进行零相位滤波 [b, a] butter(8, 100/(Fs/2), low); y filtfilt(b, a, x);filtfilt相当于对信号做了两次滤波所以等效的阻带衰减翻倍达到60dB指标实际会得到120dB的效果等效阶数也是原来阶数的两倍。这意味着你不需要把阶数设得特别高就能获得很陡的过渡带和零相位特性。但要注意filtfilt只适用于离线场景实时系统没法“反着滤波”。而且它在信号边界处会产生较大的瞬态误差处理短数据时要留意开头和结尾的失真。实际使用时我会先对信号边缘做延拓比如把信号开头和结尾的几个数据点按照原始序列的特性做镜像或线性外推滤波后再裁掉能明显减少边界畸变。5. 实测中躲不开的坑数值稳定性、边界效应与调试技巧5.1 IIR阶数一高就数值发散用butter直接设计18阶滤波器后如果你直接用[b,a]形式存系数并调用filter在高精度输入下有时会发现输出在某个点突然变成NaN或者剧烈振荡。这不是滤波器本身的问题而是传输函数直接型结构在数值上“病态”了。高阶多项式的根对系数极其敏感浮点数表示误差会随阶数指数级放大。解决方法是把滤波器转换成二阶节级联形式。Matlab里用tf2sos或zp2sos即可% 高阶IIR滤波器的正确使用姿势 [N, Wn] buttord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [z, p, k] butter(N, Wn, low); [sos, g] zp2sos(z, p, k); % 用二阶节形式滤波 y sosfilt(sos, x) * g; % 或者用filtfilt配合SOS y2 filtfilt(sos, g, x);经验是IIR阶数大于5就尽量用SOS形式不要图省事直接用[b,a]。这是一个在Matlab里看不出来、但到了嵌入式定点实现时差别巨大的细节。5.2 滤波起始段波形被拉偏用filter滤波后输出信号的最前面几十个点经常出现明显的瞬态失真波形被“拉”一下才恢复正常。原因是滤波器初始状态为零但输入信号在第一个采样点未必是零或者说信号的实际初始状态和滤波器默认的初始状态不一致滤波器需要一段时间“建立”稳态响应。三种常用处理方式一是对信号开头补一段“预热”数据滤波后裁掉二是直接丢弃输出信号前N个点N取滤波器阶数的几倍三是改用filtfilt它内部会自动处理初始条件边界效应小得多。我在处理实际采集的传感器数据时会在滤波前先截掉原始信号前几十毫秒的数据因为那部分往往还是上电瞬态本身就不包含有效信息一箭双雕。5.3 归一化频率搞错Matlab里的Wn参数大多数情况是归一化频率以奈奎斯特频率为基准取值范围0到1。比如采样率1000Hz奈奎斯特频率500Hz100Hz对应的归一化频率就是100/5000.2。很多人直接把100当Wn传进去结果滤波器截止频率变成了100倍的奈奎斯特频率设计出来的完全是另一个东西。freqz也一样如果不指定采样率横轴默认是归一化频率单位是π rad/sample。用freqz(b,a,1024,Fs)指定采样率后横轴才会变成Hz。调试时先确认坐标轴单位再判断曲线是否符合预期能省很多冤枉时间。5.4 群延迟图看出来“振荡”是正常的吗第一次看IIR滤波器的群延迟曲线时很多人会被那些在截止频率附件猛冲的尖峰吓到以为是计算错误。其实这就是IIR相位响应的真实表现。相位变化率在靠近极点位置时本来就很大群延迟出现尖峰是数学上的必然结果。判断一个滤波器相位特性的好坏要看通带内群延迟的平坦度而不是看有没有尖峰。5.5 问题速查表现象常见原因解决方案滤波后输出出现NaN或发散高阶IIR直接型数值不稳定改用zp2sos二阶节形式信号开头波形被拉偏滤波器初始状态与信号不匹配丢弃前若干点、信号延拓、改用filtfilt截止频率完全不对归一化频率计算错误确认Wn f/(Fs/2)绘图时指定采样率幅频响应达标但时域波形很怪相位失真严重改用FIR或filtfilt零相位滤波实时系统跑不动IIR阶数低但定点不稳 / FIR阶数高计算量大评估椭圆IIR、或优化FIR用firpm降阶两个滤波器幅频看起来一样单看幅频忽略相位特性增加群延迟对比观察通带内平坦度阻带衰减差一点点没到指标阶数估算公式的保守余量不足手动加2~4阶再验证不要硬凑调试滤波器时我最推荐的工具是fvtool。它可以把多个滤波器放在同一个窗口里对比显示幅频、相频、群延迟、零极点图而且能量化显示通带纹波和阻带衰减的实际值省得肉眼从曲线上估读。语法很简单fvtool(b1, a1, b2, a2, h1, 1, Fs, Fs);还可以配合measure(d)方法自动测量滤波器的各项指标是否满足设计规格输出结果一目了然。用fvtool把幅频响应和群延迟两个图并排对比比看任何文字描述都直观。我个人的习惯是设计完滤波器后永远不直接跳到信号处理先跑一遍fvtool看图再用一段标准的正弦叠加信号测试滤波效果确认幅频、相频都符合预期后再接入真实数据。这个习惯帮我避开过很多“指标全达标但波形全乱了”的隐蔽问题。不同函数设计低通滤波器这件事说到底没有绝对的“最好”只有针对你的信号类型、算力预算和相位要求的“最合适”。希望这篇对比能帮你少走点弯路。