先说明一点这个标题乍看有个小矛盾低通和带阻一个是“只留低频”一个是“掐掉中间一段”严格说不是同一种滤波器。但你把它俩放在一起并不奇怪——实际项目里传感器信号经常要先用低通剔掉高频噪声再对特定频段比如50Hz工频干扰做带阻陷波反过来设计带阻时也常常从低通原型推导。这篇就以MATLAB为主线把双线性变换法下的巴特沃斯低通滤波器和带阻滤波器从指标设置、参数计算到代码实现、实测验证完整走一遍。这篇内容主要面向正在做《数字信号处理》课程设计、信号采集预处理、或者刚接触MATLAB滤波器设计的朋友。你会看到两个工程中最常用的组合巴特沃斯滤波器的平坦幅频特性加上双线性变换的稳定映射。全文不绕弯子直接给可复现的脚本同时把容易踩的坑一并说清楚。1. 先说清楚低通和带阻为什么总被放在一起1.1 两种滤波器分别干什么活低通滤波器的作用用一句话概括就是“频率越高压得越狠”。设一个截止频率低于它的信号基本无损通过高于它的信号按固定斜率快速衰减。这个特性在信号采集里太常用了传感器出来的数据高频部分往往不是有用信息而是环境噪声、运放毛刺、开关尖峰。先做一级低通把高频垃圾清掉后面再分析慢变趋势或者做特征提取结果会干净得多。带阻滤波器则是反过来它放行低频和高频只针对某一个连续频段做大幅度衰减。最典型的场景就是50Hz工频干扰。实验室里只要设备接地不好采集出来的信号波形上大概率会叠一层50Hz的波浪纹。这时候用带阻或者陷波器就能精准压制这个频段而不去碰有用信号。注意带阻并不是只压某一个单频点它一个频带都压这跟陷波器有点区别。这两个滤波器表面上一放一挡实际上设计思路同源——很多高阶带阻滤波器本质上是先设计一个低通原型再做频率变换得到的。这就是为什么标题里“低通带阻”并列出现并不违和它其实是一条技术路线上的两种输出形态。1.2 巴特沃斯与双线性变换的“老搭档”关系先聊聊为什么偏偏是巴特沃斯。数字滤波器家族里巴特沃斯最出名的特点是通带内部“绝对平坦”——从0Hz一直到截止频率附近幅度响应几乎是一条直线不像切比雪夫滤波器那样在通带有波纹。用示波器看滤波后的波形巴特沃斯不会给你带来可察觉的通带起伏这在很多对幅值精度敏感的测量场景里非常重要。代价是它的过渡带相对平缓想要更陡的衰减斜率就得提高阶数同样指标下切比雪夫或椭圆滤波器的阶数通常更低。至于双线性变换它是把模拟滤波器“翻译”成数字滤波器的一把标准钥匙。模拟滤波器设计理论很成熟已经有几十年的积累各种归一化原型、频率变换公式都是现成的。但数字滤波器工作在离散域不能直接把模拟电路里的s域传递函数换成z域完事。双线性变换的做法是把整个模拟频率轴逐点映射到数字频率轴上[ s \frac{2}{T} \cdot \frac{1 - z^{-1}}{1 z^{-1}} ]这个变换最大的好处是只要模拟原型稳定变换出来的数字滤波器就一定稳定。这一点重要到什么程度呢脉冲响应不变法做类似映射时如果原型的阻带衰减不够数字化后高频分量会混叠回有效频带滤波器实际性能大打折扣双线性变换则通过频率轴的压缩彻底避免了混叠——代价是频率位置会被“挤歪”这就是常说的频率畸变。好在畸变有明确的数学关系设计时预先做一次反向补偿预畸变就能精确落在目标频率上。整个过程像把一张世界地图投影到平面上面积和形状会略变但你知道投影规则就能提前校准回来。MATLAB里的butter函数之所以能直接在数字域工作底层走的就是这条路先构造满足指标的模拟巴特沃斯原型再用双线性变换数字化。理解这个底层逻辑后面看代码就不会一头雾水。2. 动手前先把三个关键参数想清楚很多同学一上来就写butter拿一组随便猜的频率就滤波最后波形不对也不知道错在哪。做数字滤波器开工之前有三个参数必须先定否则后面全是白忙。2.1 采样率与归一化频率第一个是采样率fs。它不只是决定时间轴还决定了数字音频、数字信号处理系统里一切频率的基准单位。根据奈奎斯特定理能无混叠处理的最高信号频率是fs/2也就是折叠频率。MATLAB的滤波器函数一律不用你输入“多少Hz”而是输入归一化频率“除以Nyquist频率之后的比值”。比如fs1000Hz时折叠频率是500Hz那要设计一个100Hz截止的低通滤波器传给butter的截止频率就是Wn 100 / (fs/2); % 也就是100/500 0.2这个Wn最大不会超过11就代表折叠频率本身。很多新手拿5000Hz的采样率却往butter里写一个100以为那是100Hz——实际上那代表的是100倍Nyquist频率滤波器设计直接报错或者行为完全错误。这里必须强调所有MATLAB滤波器设计函数里的频率要么是归一化数字频率要么是弧度/秒很少有直接填Hz的。你需要在最前面定义一个fs变量所有频率都从它换算出来。第二个参数是有用信号的最高频率成分。比如振动信号关心0到80Hz的频段噪声主要在200Hz以上那低通截止频率选100Hz左右、阻带边界选150~250Hz是合理的。如果这个判断本身错了滤波器设计得再漂亮也没有意义。2.2 设计指标怎么写才不会被系统坑滤波器不是一根针做不到“小于100Hz全留、大于100Hz全删”。真实滤波器在通带和阻带之间有一条过渡带通带内允许有一定衰减阻带内也要求衰减足够大。这些要求就构成了设计指标通带最大衰减Ap通带里幅度允许掉的dB数通常取1dB或3dB。阻带最小衰减Ast阻带里至少要压掉多少dB工程上40dB起步要求高一点可以到60dB。通带边界频率和阻带边界频率这两个频率夹出来的区间就是过渡带。过渡带越窄滤波器阶数越高。阶数高了不是问题但系数数值稳定性、运算量都会跟着上来。所以在定指标时除非真实需求就那么苛刻否则尽量给过渡带留一点余量。用MATLAB的buttord函数算出最少需要的阶数和截止频率是高阶还是低阶一目了然。2.3 一个低通例子先算好预期我拿一个很常见的场景来说采集系统采样率1000Hz有用信号在50Hz附近200Hz以上的高频噪声需要压掉。目标指标设为通带边界100Hz、阻带边界250Hz、通带衰减1dB、阻带衰减40dB。对应的归一化频率是fs 1000; fp 100; % 通带边界 Hz fst 250; % 阻带边界 Hz Ap 1; % 通带最大衰减 dB Ast 40; % 阻带最小衰减 dB Wp fp / (fs/2); Ws fst / (fs/2);照这个参数算出来的结果在我的环境下是6阶巴特沃斯归一化截止频率大约0.24。阶数不高运算负担很小滤波后100Hz以内几乎无损250Hz以上衰减超过40dB符合直觉预期。如果你把阻带边界改成150Hz过渡带立刻变窄阶数可能直接跳到13阶左右——这就是我前面说的“别把过渡带交界处压太近”。3. 低通滤波器实现与双线性变换的手动验证3.1 用buttord加butter一把梭参数算好之后设计低通滤波器的代码短得惊人% 低通巴特沃斯滤波器设计 [n, Wn] buttord(Wp, Ws, Ap, Ast); [b, a] butter(n, Wn, low); fprintf(设计阶数 n %d\n, n); fprintf(归一化截止频率 Wn %.4f\n, Wn);buttord根据阻带衰减要求倒推阶数butter再基于阶数和截止频率生成分子分母系数。b是分子系数a是分母系数后面直接用filter或者freqz就能派上用场。如果想看一眼幅频特性是否达标[h, f] freqz(b, a, 4096, fs); plot(f, 20*log10(abs(h))); grid on; xlabel(频率 (Hz)); ylabel(幅度 (dB)); xlim([0 fs/2]);你会看到典型的巴特沃斯曲线通带内几乎平直过了截止频率后以大约每十倍频120dB的速度滚降6阶对应-120dB/dec在250Hz处确实压到了-40dB以下。3.2 通过buttap加bilinear走一遍底层路线用一行butter就能出结果为什么还要学底层路线因为面试、答辩或者课程报告里老师很可能追问一句“双线性变换体现在哪”你总不能在代码里找到一行bilinear也没有。所以下面这段代码值得仔细看它和butter做的事完全等价只是把内部步骤拆开了。先对两个边界频率做预畸变从数字频率反向换算成模拟角频率% 双线性变换手动版 wp_analog 2 * fs * tan(pi * fp / fs); ws_analog 2 * fs * tan(pi * fst / fs); % 在模拟域计算阶数 [n2, Wn2] buttord(wp_analog, ws_analog, Ap, Ast, s); % 生成归一化巴特沃斯模拟原型 [z, p, k] buttap(n2); [ba, aa] zp2tf(z, p, k); % 把原型截止频率从1 rad/s缩放到Wn2 [ba, aa] lp2lp(ba, aa, Wn2); % 双线性变换数字化 [bd, ad] bilinear(ba, aa, fs);逐行解释一下。tan(pi * fp / fs)就是双线性变换的预畸变核心先把目标频率“提前拉偏”让数字化之后刚好回到100Hz原处。buttap(n2)生成的是截止频率为1 rad/s的归一化巴特沃斯原型系数很标准但截止频率不对所以用lp2lp把它缩放到Wn2。最后bilinear把模拟传递函数映射成数字传递函数得到bd、ad。3.3 两种方法的结果一致性验证把butter得到的b、a和手动路线得到的bd、ad放在一起画幅频响应两条曲线基本重合。用MATLAB的比较方式[h1, f] freqz(b, a, 4096, fs); [h2, ~] freqz(bd, ad, 4096, fs); plot(f, 20*log10(abs(h1)), b-, LineWidth, 1.5); hold on; plot(f, 20*log10(abs(h2)), r--, LineWidth, 1.5); legend(butter 直接设计, 手动双线性变换);我这里两条曲线最大误差在-80dB以下肉眼几乎看不出来。这说明一个问题MATLAB的butter对你封装了很多细节但核心思想就是“模拟原型双线性变换”。理解这条底层路线后面碰着需要手动调整频率配置、或者要自己写zet域的推导题时你就不慌了。4. 带阻滤波器实战用巴特沃斯压掉50Hz工频干扰4.1 buttord和butter的带阻模式用法与坑带阻的设计思路跟低通一样但参数变成了两组频率。还是用采样率1000Hz来举例我需要在45Hz到55Hz之间压出至少40dB的衰减同时20Hz和80Hz附近的有用信号要保持住。于是设计指标写成通带边界是20Hz和80Hz阻带边界是40Hz和60Hz。% 带阻巴特沃斯滤波器设计 Wp_pass [20 80] / (fs/2); % 通带边界 Ws_stop [40 60] / (fs/2); % 阻带边界 [ns, Wns] buttord(Wp_pass, Ws_stop, 1, 40); [bs, as] butter(ns, Wns, stop);这里有个特别容易翻车的点对带阻滤波器Wp和Ws的顺序不等于“小频率在前、大频率在后”这么简单而是必须体现“通带在外、阻带在内”的关系。也就是说Wp的两个值必须比Ws的两个值更靠外0.04和0.16在两边0.08和0.12在中间。如果你条件反射把buttord的第一参数写成阻带范围第二参数写成通带范围MATLAB会直接报“频率范围不一致”错误或者干脆返回一个完全反了的滤波器——把有用的20Hz和80Hz信号削掉而50Hz纹丝不动。我第一次做带阻就吃过这个亏后来总结了一句口诀通带宽、阻带窄阻带永远被通带包在中间。butter函数里同样要指定滤波类型为stop千万别漏。4.2 对信号做滤波前后的对比带阻设计完随便造一段含50Hz干扰的信号实验一下t (0:fs-1)/fs; x 0.8 * sin(2*pi*20*t) 0.6 * sin(2*pi*50*t) 0.5 * sin(2*pi*120*t); y filter(bs, as, x);这组信号里20Hz和120Hz是有用成分50Hz是模拟工频干扰。滤波后看时域波形20Hz和120Hz的幅值变化不大50Hz那个明显的波动明显被压下去了。再用freqz看频域40到60Hz之间是深谷峰值处衰减超过40dB而20Hz和80Hz两侧只有不到1dB的损耗完全符合设计指标。这里顺便解释为什么工频干扰场合用带阻而不是普通低通。50Hz这种频率往往在有用信号的频带之内低通一刀切会把有用信号也干掉带阻则只针对窄带做处理两边都能保全。如果现场测量发现50Hz有轻微的频率漂移带阻的带宽还留了缓冲空间如果非要压单个频点也可以考虑iirnotch但它对频率漂移更敏感工程上我倾向先用带阻。4.3 由低通原型到带阻的变换原理与手动示例带阻与低通之间有一条经典的设计捷径先做一个低通原型再用频率变换把它“折叠”成带阻。对应MATLAB里就是lp2bs函数。给一个手动示意版本% 4阶低通原型 n4 4; [z4, p4, k4] buttap(n4); [b4, a4] zp2tf(z4, p4, k4); % 定义中心频率和带宽这里先不预畸变仅示意变换过程 W0 2 * pi * 50; % 中心频率 50Hz BW 2 * pi * 10; % 带宽 10Hz % 低通原型 - 模拟带阻 [bbs, abs] lp2bs(b4, a4, W0, BW); % 双线性变换数字化 [bd_bs, ad_bs] bilinear(bbs, abs, fs);这套代码能跑通频率响应也能看到阻带但我要提醒一句由于没有做严格的频率预畸变数字化后阻带中心会和50Hz存在偏差。设计带阻的严谨做法是先对阻带上下边界fl、fh做预畸变再用预畸变后的值换算中心频率和带宽wl 2 * fs * tan(pi * 45 / fs); wh 2 * fs * tan(pi * 55 / fs); W0 sqrt(wl * wh); BW wh - wl;用这套参数替换上面的W0和BW再走lp2bs加bilinear阻带就能精确落在45~55Hz。平时图省事直接用buttord加butter当然更简单这个手动流程的主要价值在于帮你理解“带阻到底是怎么从低通变出来的”考试和答辩都爱问这个。5. 实测验证与波形分析5.1 用合成信号验证低通的降噪效果光看幅频曲线还不够真实信号走一遍更直观。造一个模拟传感器信号50Hz是有用信号300Hz和500Hz是高频噪声再加一点随机白噪声。用前面设计好的6阶巴特沃斯低通滤波t (0:fs-1)/fs; x 0.5 * sin(2*pi*50*t) 0.3 * sin(2*pi*300*t) ... 0.2 * sin(2*pi*500*t) 0.05 * randn(size(t)); y filter(b, a, x);滤波前后画在同一张图上明显看到高频毛刺被清理掉50Hz的正弦波形重新变得光滑。再从频域看频谱X abs(fft(x)); Y abs(fft(y)); f_fft (0:length(t)-1) * fs / length(t);150Hz以上几乎没有残留300Hz和500Hz两个尖峰被压到接近底噪。这个实验结果和设计指标完全对得上说明这套滤波器在工程场景里是靠谱的。5.2 分析滤波器的通带平坦度与相位响应巴特沃斯滤波器的优点在通带平坦度上体现得很明显。画幅频曲线时100Hz前的增益波动不超过0.5dB在测试测量场景里可以忽略。代价是相位响应并非线性的——所以用filter做滤波之后输出波形相对输入会有不同程度的相位延迟不是每个频率延迟一样。如果只关心幅值那无所谓如果后面要做波形对齐、时延估计这个相位失真就会成为问题。解决办法是用零相位滤波filtfilt正向滤波一次再把序列倒过来反向滤波一次相位抵消波形位置基本不变y_zero filtfilt(b, a, x);但注意filtfilt有两个代价一是它是非因果的不能用于实时在线处理二是在序列的起始和结束位置容易出现预振铃或边界过渡失真。数据量很短时要特别小心最好掐头去尾再分析。5.3 实际项目中的数据场景我在一个采集声发射信号的小项目里用过类似配置。传感器输出的原始信号里目标频段集中在100kHz以下但采集系统本身引入了大量射频噪声和高频振铃直接分析根本看不出特征。我当时就是先用一个巴特沃斯低通把高频糊掉再用带阻把电路板带来的某个特定干扰频率压掉。那套流程跟上面写的几乎没有区别唯一的差别是采样率更高、阶数稍微调大了一点。所以不要觉得示例里的1000Hz采样率是拍脑袋换成任何采样率只要先把Nyquist频率算清楚思路完全不变。6. 我踩过的坑与给初学者的速成经验6.1 参数顺序、归一化频率与画图习惯第一个坑在前面反复提了带阻的buttord第一参数是通带第二参数才是阻带顺序写反轻则报错重则给你一个反着干活儿的滤波器。第二个坑是频率单位butter、buttord、freqz对频率的理解各不相同butter要归一化频率freqz可以被喂一个采样率参数但传进去的数值必须是Hz。如果混合用画出来的横轴经常就是0到1的归一化刻度而不是Hz看着看着就懵了。我自己的习惯是代码开头就用fs 1000定义一个变量所有频率全部写成Hz数值需要归一化时临时除一下既好读又不会错。另外一个小习惯值得养成设计完滤波器不要急着对信号做处理先把幅频响应画出来看截止频率、阻带衰减是否落在预期位置。这一步最快发现问题比处理完数据再回头查快得多。6.2 高阶滤波器的数值稳定性问题巴特沃斯阶数高了以后直接用b、a做filter容易出问题。阶数超过十几阶直接型结构对系数舍入误差非常敏感滤波器响应可能出现顶端毛刺甚至不稳定。碰到这种情况不要硬扛用二阶节级联形式sos zp2sos(z, p, k); y sosfilt(sos, x);zp2sos把高阶系统拆成一堆二阶子系统每个子系统的系数都控制在很小的动态范围里数值稳定性大幅提升。这也是为什么很多MATLAB老手设计高阶滤波器时很少直接看b、a而是去看sos。当然如果项目指标允许优先把过渡带放宽一点让阶数控制在10以内问题从源头就消失了。6.3 什么时候不该用巴特沃斯巴特沃斯不是万能的。如果指标里要求过渡带特别陡、阶数被压到15以上说明你选错了滤波器家族。同样一组指标切比雪夫I型可能只需要8阶椭圆滤波器甚至只要5阶。巴特沃斯的优点在于通带平坦但平坦是要用过渡带宽度去换的。工程中如果对通带纹波要求不高、但对实时性、对算力敏感果断换切比雪夫或椭圆。反过来如果要求幅值精度极高、哪怕0.1dB的波纹都不可接受那就老老实实用巴特沃斯通过放宽过渡带来解决问题。最后再分享一个我常用的实操技巧在写正式代码之前先用MATLAB自带的滤波器设计工具filterDesigner或者fdatool把指标拖一拖、看一看。工具界面里左边填通带、阻带、衰减、滤波器类型右边实时显示幅频曲线、零极点图。先用它确认“这个指标确实能实现、阶数能接受”再回命令行脚本复用数值能少走很多弯路。毕竟双线性变换和巴特沃斯明白了原理之后剩下的功夫就是在参数上多试几次。