前阵子帮一个做嵌入式采集的哥们调试传感器信号他拿着MATLAB问我为什么我设计的滤波器滤波后波形还是乱的我看了一眼代码发现他把模拟滤波器的系数直接套到filter函数里去处理数字信号这当然不对。后来我把模拟滤波器设计、数字IIR设计、数字FIR设计从头到尾跟他过了一遍也正好把“基于MATLAB的模拟滤波器和数字滤波器设计”这个经典课题整理了出来。这篇文章围绕低通、高通以及带通、带阻这些最常见需求讲清楚指标怎么定、原型怎么选、MATLAB代码怎么写、结果怎么验证。如果你是做课程设计或者刚接触信号处理可以直接照着流程走。1. 设计前的指标拆解与整体思路1.1 先分清模拟滤波器和数字滤波器设计路线完全不同模拟滤波器处理的是连续时间信号传递函数是s域的H(s)物理上通常用电阻、电容、电感或者运放搭建数字滤波器处理的是离散时间信号传递函数是z域的H(z)实际运行在DSP、FPGA或者PC软件里。MATLAB里设计这两类滤波器的函数名字虽然很像但参数含义、设计流程差别非常大。我见过太多新手把一个s域滤波器的系数直接丢给filter或者filtfilt去处理采样后的数据结果出来的波形不是畸变就是根本不起作用。原因很简单模拟滤波器的系数描述的是连续系统数字滤波器的系数描述的是差分方程两者不能用同一个函数直接对应。模拟设计用freqs看频响数字设计用freqz看频响调用的工具都不一样。所以设计之前第一件事是想清楚你的信号是已经采回来的离散序列还是要做硬件前端的抗混叠滤波如果信号已经变成了数字采样点那就老老实实设计数字滤波器如果你是搭电路之前先仿真模拟滤波器的可行性那再走s域路线。很多时候一个项目里两者都会用到但它们是串联关系不是替换关系。1.2 指标设置不要只给一个截止频率很多教程里“设计一个截止频率为1kHz的低通滤波器”这种说法其实是不完整的。滤波器设计需要一组完整的指标至少包括四个通带截止频率比如0到4kHz需要保留阻带截止频率比如6kHz以上需要被抑制通带最大纹波比如不超过1dB阻带最小衰减比如不低于40dB。这四个指标一起决定了滤波器阶数和过渡带宽度。单纯给一个截止频率设计出来的滤波器可能过渡带特别宽通带纹波超标或者阻带衰减根本达不到系统要求。举个例子假设一个信号采样率Fs 8000Hz需要保留1000Hz以下的低频成分滤掉1500Hz以上的高频成分要求通带纹波不超过1dB阻带衰减不低于40dB。这套指标看起来简单但后面所有MATLAB设计函数都靠它来算阶数、算归一化频率。指标越严阶数越高计算量越大实际工程里需要反复权衡而不是一味追求陡峭的过渡带。2. 模拟滤波器设计从归一化原型到实际截止频率2.1 原型选择巴特沃斯、切比雪夫、椭圆怎么选模拟滤波器设计里最常用的四类原型是巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆滤波器。它们的频率响应各有特点选型不能只看名字好不好听。巴特沃斯滤波器最突出的优点是通带内幅频特性非常平坦没有纹波但代价是过渡带比较宽。也就是说从通带到阻带需要更宽的频率变换范围才能达到同样的衰减指标。切比雪夫I型滤波器允许通带内有等纹波阻带是单调下降的相同阶数下过渡带比巴特沃斯陡适合对通带纹波要求不是特别苛刻的场景。切比雪夫II型反过来通带平坦阻带内有纹波。椭圆滤波器在通带和阻带里都有纹波但是相同阶数下过渡带最陡能以最低的阶数满足最严的指标代价就是相频特性和参数敏感度都比较复杂。MATLAB里对应的函数很直观butter、cheby1、cheby2、ellip。具体选哪个我的经验是课程设计和多数工程验证优先用巴特沃斯因为参数简单、结果好解释如果阶数受限制、过渡带要求比较严再切到切比雪夫I型或椭圆。下面这张表可以帮你快速判断原型类型通带特性阻带特性过渡带陡峭程度典型场景巴特沃斯最平坦单调下降较宽对通带平坦度要求高切比雪夫I型等纹波单调下降中等允许少量通带纹波切比雪夫II型平坦等纹波中等对阻带纹波不敏感椭圆等纹波等纹波最陡对阶数有严格限制2.2 模拟低通滤波器的MATLAB实现步骤模拟滤波器的标准设计流程是先确定低通原型再用频率变换得到高通、带通或带阻。MATLAB里其实已经把这些变换封装好了直接调用即可。比如我要设计一个模拟低通滤波器通带边界1000Hz阻带边界1500Hz通带纹波1dB阻带衰减40dB。注意模拟滤波器的频率单位不是Hz而是角频率rad/s。MATLAB的模拟设计函数里Wp和Ws默认用的就是角频率。所以需要先做一次换算Fs_analog 2 * pi * 1000; % 通带边界角频率 Fp 2 * pi * 1000; % 1000Hz - rad/s Fs_stop 2 * pi * 1500; % 1500Hz - rad/s Rp 1; % 通带纹波 1dB Rs 40; % 阻带衰减 40dB % 计算满足指标的最低阶数 [N, Wn] buttord(Fp, Fs_stop, Rp, Rs, s); [B, A] butter(N, Wn, s); % 绘制幅频响应 [H, W] freqs(B, A, 4096); figure; plot(W / (2 * pi), 20 * log10(abs(H))); xlabel(频率 (Hz)); ylabel(幅度 (dB)); grid on;这里buttord返回的N是最低阶数Wn是自动修正后的3dB截止频率。butter生成了分子系数B和分母系数A它们是s域传递函数的系数。用freqs可以在模拟频率轴上计算频响注意横坐标要除以2*pi才能显示成Hz。如果要做高通或者带通同样可以给butter传额外参数。例如设计一个高通模拟滤波器截止频率500Hz[B, A] butter(N, 2 * pi * 500, high, s);实际上MATLAB支持直接指定滤波器的类型low、high、bandpass、stop都可以bandpass时Wn[W1 W2]是一个二元向量。这里要特别提醒一点模拟滤波器的系数数量级常常很大因为角频率动辄几千甚至上万导致系数矩阵病态。画频响没问题但如果要转化成数字滤波器最好还是专门走数字设计流程不要手动把模拟系统离散化后面会讲为什么。3. 数字IIR滤波器设计双线性变换的核心逻辑3.1 为什么IIR要用双线性变换而不是冲激响应不变法数字IIR滤波器最经典的设计方法是从模拟滤波器原型变换过来。常见的变换方法有两种冲激响应不变法和双线性变换法。冲激响应不变法的思路是直接对模拟滤波器的冲激响应采样得到数字滤波器的脉冲响应。它的优点是数字滤波器的时域特性和模拟原型保持得比较好但问题是会发生频谱混叠。如果模拟原型的阻带衰减不够陡混叠的成分会直接叠加到数字滤波器的通带里导致设计结果不达标。更重要的是这种方法对高通和带阻滤波器是完全不实用的因为高通在奈奎斯特频率附近有很高的频率分量混叠会非常严重。双线性变换法则完全不同。它把s平面通过正切函数映射到z平面核心思路是用一个数值积分近似替代微分关系。这种变换不会产生频谱混叠代价是频率轴发生了非线性压缩也就是说数字滤波器实际的转折频率和模拟原型的转折频率不再相等。好在非线性关系是确定的可以通过预畸变公式把设计目标频率反过来修正。实际使用MATLAB时butter、cheby1这些函数在指定设计数字滤波器时内部会自动完成双线性变换和频率预畸变不需要自己手动做。但理解这个原理很重要否则一旦需要自定义模拟原型或者要用bilinear函数手动变换就可能搞不清为什么频率总差一截。3.2 IIR低通与高通的完整MATLAB代码直接设计数字IIR低通滤波器沿用前面那组真实工程指标采样率Fs 8000Hz通带边界1000Hz阻带边界1500Hz通带纹波1dB阻带衰减40dB。Fs 8000; % 采样率 Wp 1000 / (Fs / 2); % 归一化通带边界单位是奈奎斯特频率 Ws 1500 / (Fs / 2); % 归一化阻带边界 Rp 1; Rs 40; % 计算最低阶数 [N, Wn] buttord(Wp, Ws, Rp, Rs); % 设计数字低通IIR滤波器 [B, A] butter(N, Wn, low); % 查看幅频和相频响应 freqz(B, A, 1024, Fs);这里的关键点是归一化频率。数字滤波器设计中所有频率都要除以Fs/2也就是奈奎斯特频率。所以1000Hz在8kHz采样率下对应0.25但MATLAB里的归一化范围是0到11对应的就是4000Hz。很多人第一次用的时候直接把1000传给butterMATLAB会把这个数当作比1还大的归一化频率设计结果完全不对。高通设计几乎一模一样把low改成high就行[B, A] butter(N, Wn, high);带通或者带阻也类似把Wp和Ws换成两元素的向量。比如要保留1000到2000Hz滤掉其他频率可以这样Wp [1000 / (Fs/2), 2000 / (Fs/2)]; Ws [800 / (Fs/2), 2400 / (Fs/2)]; [N, Wn] buttord(Wp, Ws, Rp, Rs); [B, A] butter(N, Wn, bandpass);3.3 阶数选择与稳定性排查IIR滤波器阶数高的时候容易出问题。理论上双线性变换会保证模拟域稳定的系统变换到数字域仍然稳定即极点都落在单位圆内但实际数值计算时高阶系数的精度损失可能导致极点跑出单位圆。最直观的排查办法是用zplane画零极点图或者直接计算极点的模[z, p, k] tf2zp(B, A); if max(abs(p)) 1 warning(极点不在单位圆内滤波器不稳定); end如果发现不稳定通常有两种处理方式。一是降低阶数重新放宽指标二是使用更稳定的滤波器结构比如在FPGA上实现时用二阶节级联而不是直接使用高阶传递函数。MATLAB里的tf2sos可以把高阶系统分解成若干二阶节级联实际硬件实现时非常常用sos tf2sos(B, A);我自己调试时发现同样一组指标下IIR滤波器的阶数往往比FIR低很多。比如上面这个低通巴特沃斯8000Hz采样率、1000Hz通带边界可能只需要四阶就能满足衰减要求换成FIR用窗函数法可能需要几十阶。这是IIR的优势但别高兴太早非线性的相位特性有时候会让它不适用于某些场景这个后面专门说。4. 数字FIR滤波器设计窗函数法和频率采样法4.1 FIR为什么能实现线性相位FIR滤波器的脉冲响应是有限长的它是直接对某个期望频率响应做逆傅里叶变换然后截断成有限长度的系数。由于没有反馈结构FIR滤波器天生是稳定的而且只要系数满足对称性就能实现严格的线性相位。线性相位意味着所有频率成分经过滤波器后产生的群延迟是恒定的波形不会因为相位非线性而畸变。这一点在音频处理、通信信号整形、心电图滤波等场景里特别重要。IIR滤波器虽然阶数低、计算效率高但相位特性往往是弯曲的各频率成分到达时间不一致输出波形看起来就会“变味”。所以很多对波形保真度有要求的项目宁可多付出一些计算量也要用FIR。窗函数法的设计思路很直观理想低通滤波器在频域是一个矩形门做逆傅里叶变换得到的时域脉冲响应是一个无限长的sinc函数。显然不能直接使用无限长序列所以用一个窗函数去截断它。加窗会让频响出现过渡带和旁瓣不同窗函数对过渡带宽度和旁瓣衰减的权衡不同。矩形窗最窄主瓣但旁瓣很高汉宁窗、汉明窗、布莱克曼窗则逐步用更宽的过渡带换取更低的旁瓣。4.2 用fir1实现FIR低通和高通滤波MATLAB里的fir1封装了窗函数法最简单的用法是Fs 8000; fc 1000; N 32; % 阶数实际抽头数是 N1 b fir1(N, fc / (Fs/2), low, hamming(N1)); freqz(b, 1, 1024, Fs);这里的N是滤波器阶数也就是延迟的采样点数。高频截止频率除以Fs/2归一化和IIR设计一样。hamming(N1)是指定窗函数长度必须是N1这一点很容易被忽略。如果省略窗函数参数fir1默认使用汉明窗。高通设计也很简单b fir1(N, fc / (Fs/2), high, hamming(N1));带通、带阻同理b fir1(N, [1000/(Fs/2), 2000/(Fs/2)], bandpass, hamming(N1)); b fir1(N, [1000/(Fs/2), 2000/(Fs/2)], stop, hamming(N1));这里要注意一个问题如果N是奇数fir1对于高通和带阻滤波器可能会自动把阶数加1因为高通FIR要实现线性相位其脉冲响应必须满足特定的对称性条件。所以设计完成后最好用length(b)确认一下实际抽头数。4.3 用kaiserord根据指标估算阶数fir1需要自己给阶数但指标到底需要多少阶拍脑袋给一个数字经常造成两种情况阶数给低了过渡带不够陡阻带衰减不足阶数给高了计算量浪费延迟变大。更专业的做法是用kaiserord根据通带边界、阻带边界和允许的纹波自动估算最低阶数。还是用那组指标通带边界1000Hz阻带边界1500Hz采样率8000Hz通带纹波1dB阻带衰减40dB。通带纹波1dB对应线性幅度的允许偏差大约是0.05阻带衰减40dB对应偏差是0.01。代码可以这样写Fs 8000; fcuts [1000 1500]; devs [0.05 0.01]; % 通带纹波1dB - 0.05阻带40dB - 0.01 [N, Wn, beta, ftype] kaiserord(fcuts, [1 0], devs, Fs); b fir1(N, Wn, ftype, kaiser(N1, beta)); freqz(b, 1, 1024, Fs);这里kaiserord返回的N是估算出的阶数Wn是归一化边界频率beta是凯泽窗的形状参数ftype是滤波器类型。用凯泽窗的好处是它可以在主瓣宽度和旁瓣衰减之间连续调节比固定窗更灵活尤其当阻带衰减要求不是标准窗能精确满足时。窗函数法之外还有频率采样法fir2可以设计任意幅频响应的FIR滤波器。不过日常低通、高通、带通设计用fir1配合kaiserord已经足够。5. 模拟与数字滤波结果对比实战5.1 如何验证滤波器到底有没有用设计完滤波器不能只看频响曲线好看还要放到真实信号里试一下。我用一个简单例子来说明信号是一个500Hz的正弦波叠加了50Hz工频干扰和2000Hz的高频噪声。我期望保留500Hz成分滤掉50Hz和2000Hz。生成测试信号并分别用IIR和FIR滤波Fs 8000; t 0:1/Fs:1; x 1.5 * sin(2 * pi * 500 * t) 0.8 * sin(2 * pi * 50 * t) 0.5 * sin(2 * pi * 2000 * t); % 设计带通数字滤波器保留500Hz附近 Wp [400 / (Fs/2), 600 / (Fs/2)]; Ws [200 / (Fs/2), 1000 / (Fs/2)]; [N, Wn] buttord(Wp, Ws, 1, 40); [B, A] butter(N, Wn, bandpass); % 使用filter滤波 y_iir filter(B, A, x); % 设计FIR带通先用kaiserord估算 fcuts [200 400 600 1000]; devs [0.01 0.05 0.05 0.01]; [Nf, Wn, beta, ftype] kaiserord(fcuts, [0 1 0], devs, Fs); bf fir1(Nf, Wn, ftype, kaiser(Nf1, beta)); y_fir filter(bf, 1, x); % 对比时域波形和频谱 subplot(2,1,1); plot(t, y_iir); subplot(2,1,2); plot(t, y_fir);这里要注意模拟滤波器不能像这样直接对数字信号用filter因为filter输出的是数字滤波器的差分方程结果。如果想验证模拟滤波器应该用freqs观察频响或者用连续时间信号通过lsim仿真。但实际工程里如果输入信号已经采样成离散序列就应该直接用数字设计路线否则就会出现开头我朋友那种错误。5.2 群延迟、相位和计算量的实际差异用freqz可以看到相位响应IIR滤波器在高频和低频附近的相移变化往往很大群延迟也不固定。用grpdelay可以直观地看到滤波后各频率成分的到达时间差。grpdelay(B, A, 1024, Fs); grpdelay(bf, 1, 1024, Fs);FIR线性相位滤波器的群延迟是一个常数等于N/2个采样周期。这意味着所有频率成分被延迟了相同时间波形形状不会发生相位扭曲。IIR滤波器的群延迟是会波动的在设计带通或高通滤波器时尤其明显。我自己做音频均衡器时发现级联多个IIR滤波器后整体相位特性会变得很复杂输出波形和原始波形差距很大这时候换成等阶数的FIR虽然延迟更大但处理效果反而“干净”得多。从计算量来看IIR因为阶数低单位采样点要做的乘加运算也少在嵌入式实时处理时更省资源。但要在FPGA或者DSP上实现线性相位滤波选择FIR会更方便因为系数对称性可以大幅减少乘法器数量。下面这张表是我平时做选型决策时的参考对比维度IIR滤波器FIR滤波器阶数需求低计算量小高计算量大相位特性非线性可实现线性相位稳定性需要检查极点天然稳定设计方法从模拟原型变换窗函数法/频率采样典型应用资源受限的控制系统通信、音频、生物电信号6. 常见问题与排查技巧实录6.1 截止频率单位混乱结果完全不对这是最高频的坑。模拟滤波器设计要用角频率单位是rad/s如果直接传Hz设计出来的截止频率会差一个2*pi倍。数字滤波器设计要用归一化频率单位是奈奎斯特频率倍数所有Hz都要除以Fs/2。很多人拿数字设计函数butter去设计模拟滤波器或者反过来结果就是滤波器带宽完全不符合预期。我建议每次写代码前先做一个简单的换算表比如在注释里写明目标频率和归一化值目标频率(Hz)采样率(Fs)归一化频率100080000.25400080001.050080000.1256.2 滤波起始段有很大瞬态波形开头明显跳变这是因为filter函数默认初始状态为零相当于滤波器在第一个采样点突然被激励输出自然会产生一段暂态响应。对于无限冲激响应的IIR滤波器这段暂态可能持续好几个时间常数。解决办法有两种一是用filtfilt做零相位滤波它对数据进行正向和反向两次滤波消除了相位偏移但延迟会加倍且要求滤波器对反向数据同样稳定二是用filtic估算合适的初始状态。% 零相位滤波 y_zero filtfilt(B, A, x); % 计算初始状态 zi filtic(B, A, y_initial, x_initial); y_filtered filter(B, A, x, zi);FIR滤波器因为长度有限起始瞬态的持续时间等于滤波器长度如果阶数很高平滑效果会更好但响应延迟也更明显。对于实时系统filtfilt是不能用的因为它需要整段数据只能在离线处理时使用。6.3 滤波后幅度衰减严重和期望不一致滤波器在通带内的增益理论上应该是0dB但实际设计时如果阶数不足、过渡带太宽或者通带边界取得太靠近截止点通带边缘幅度就会有明显衰减。检查方法很简单用freqz看通带内的幅度是不是接近1。如果只是为了找到一个滤波器让整体效果更干净可以在设计后手动校准增益把系数乘以一个修正因子。另一个常见错误是设计高通或带通时把通带和阻带边界顺序写反了。比如带通要求Wp[400 600]结果写成Wp[600 400]MATLAB会报错或者设计出奇怪的滤波器。遇到这种情况先检查边界频率大小关系。6.4 常见问题速查表现象可能原因解决方法滤波器完全没效果模拟系数用在数字信号上重新设计数字滤波器截止频率偏移归一化频率算错用fc/(Fs/2)波形开头混乱filter零初始状态用filtfilt或filtic阻带衰减不够阶数太低提高N或用椭圆/凯泽窗出现NaN或无穷大极点不稳定降低阶数或查看极点相位失真严重IIR非线性相位改用FIR带通边界颠倒频率向量顺序错检查Wp、Ws升序最后再分享一个小技巧设计过程中多用fvtool这个图形化工具它可以直接查看幅频、相频、群延迟和零极点图。我一般会让代码先输出到变量然后运行fvtool(B,A)所有波形一屏看全比自己写绘图代码省事得多。另外如果只是临时验证一个滤波器效果也可以打开MATLAB自带的交互式设计工具filterDesigner把参数填进去再导出系数然后放到自己的项目里。但要注意交互式工具生成的代码风格和手写的有些差异最好还是理解原理后再决定要不要直接用。