1. 从一道面试题说起为什么数字滤波器非要绕一道模拟的弯几年前我去某通信公司面试对方出了道题给一个采样率16kHz的系统要做截止频率2kHz的低通滤波要求通带纹波小于1dB阻带衰减大于40dB设计一个IIR滤波器并写出C语言实现。我当时第一反应是直接上Matlab的butter函数几秒钟就能算出系数。但面试官追了一句系数怎么来的双线性变换这一步到底在干什么——这才问到根子上了。很多做嵌入式、FPGA或者音频处理的朋友用FIR比较多遇到IIR时多半是拿现成工具箱点一下系数抄出来用。但对双线性变换这五个字往往停留在大概是把模拟滤波器变成数字滤波器的一种方法这种模糊认知。等到指标跑偏了、相位炸了、或者滤波器不稳定时根本不知道去哪排查。今天我就把这个过程彻底拆开从数学原理到C语言落地再到FPGA实现踩过的坑一次讲透。先说清楚双线性变换的作用它是一个从连续时间系统模拟滤波器映射到离散时间系统数字滤波器的桥梁。为什么需要这道桥因为模拟滤波器的设计理论极其成熟——Butterworth、Chebyshev、椭圆滤波器从归一化查表到极点配置前人把路走平了。而数字IIR滤波器想直接在设计域求解数学复杂度很高。更务实的路径是先在模拟域设计好一个满足指标的滤波器再通过某种映射把它翻译成数字域的差分方程。双线性变换就是最常用、最稳的翻译器。这篇文章适合谁看正在写滤波器代码但只用工具箱的工程师准备面试信号处理岗位的学生以及做FPGA数字信号处理、想把IIR滤波器跑在硬件上的人。我会从推导开始给出完整的设计流程、C语言代码、FPGA注意事项和踩坑实录。2. 双线性变换的来龙去脉一个梯形公式引发的映射2.1 从拉普拉斯变换到Z变换为什么不能直接代换先回顾一个基础问题模拟滤波器的传递函数是H(s)数字滤波器的传递函数是H(z)那能不能直接把s换成z教科书上有个经典映射是z e^{sT}也就是把s域的jΩ轴映射到z域的单位圆。从数学上这是理想的映射但问题在于e^{sT}展开后是超越函数没法直接变成有理多项式也就没法得到形如(b0 b1z^-1 ...)/(1 a1z^-1 ...)的差分方程系数。那退而求其次能不能用某种有理函数来近似e^{sT}回想一下数值分析里的帕德近似——e^x的[1,1]阶帕德近似是(1 x/2)/(1 - x/2)。把x换成sT得到z e^{sT} ≈ (1 sT/2) / (1 - sT/2)反解出s关于z的表达式s ≈ (2/T) * (1 - z^-1) / (1 z^-1)这就是双线性变换的数学出身。说到底它就是用一阶有理函数近似指数函数。低阶近似的好处是映射简单、稳定坏处是引入了频率响应的非线性畸变这个后面细说。2.2 另一种推导视角梯形法数值积分如果觉得帕德近似太数学换个更工程化的角度来理解。考虑模拟积分器y(t) ∫x(τ)dτ其传递函数是1/s。对微分方程做数值积分用梯形法逼近定积分。设采样周期T从时刻(n-1)T到nT积分y[n] y[n-1] (T/2) * (x[n] x[n-1])两边做Z变换Y(z) z^-1 Y(z) (T/2) * (X(z) z^-1 X(z))整理得Y(z)/X(z) (T/2) * (1 z^-1) / (1 - z^-1)这个离散积分器的传递函数如果要让它对应模拟的1/s那映射关系就是1/s ←→ (T/2) * (1 z^-1) / (1 - z^-1)取倒数就是s (2/T) * (1 - z^-1) / (1 z^-1)。跟我前面用帕德近似推出来的结论一致。用梯形法积分这个视角对后面理解频率畸变特别有帮助——梯形法本身就会对高频分量产生积分误差。2.3 映射前后的频率对应关系与预畸变双线性变换最关键、也最容易被忽略的特性是模拟角频率Ω和数字角频率ω之间的关系不是线性的。把s jΩ代入映射式jΩ (2/T) * (1 - e^{-jω}) / (1 e^{-jω})用欧拉公式化简分子分母同乘e^{jω/2}jΩ (2/T) * (e^{jω/2} - e^{-jω/2}) / (e^{jω/2} e^{-jω/2}) (2/T) * (2j sin(ω/2)) / (2 cos(ω/2)) j (2/T) tan(ω/2)所以得到频率变换公式Ω (2/T) * tan(ω/2)这条公式的信息量极大。当ω很小时tan(ω/2) ≈ ω/2所以Ω ≈ ω/T频率近似线性当ω接近π时tan(ω/2)趋于无穷大映射把整个模拟频率轴的无穷远点压到了数字频率的π处。这个非线性意味着什么呢如果你在模拟域设计一个截止频率Ωc 1000 rad/s的滤波器直接做双线性变换后数字域的实际截止频率不是期望的ω Ωc*T而是被压缩到了更低的频率。反过来也一样——你想让数字滤波器在ω处有某个特征映射回模拟域的时候要先预畸变把设计频率拉高一些。具体做法Ω_设计 (2/T) * tan(ω_期望 / 2)这就是频率预畸变公式。我在实际项目里看到无数人栽在这里直接用Ω ω/T去设计模拟原型变换完发现截止频率偏了5%到10%通带纹波看着也不对。原因就是忘了预畸变。后面我会用完整算例演示正确的处理方式。2.4 稳定性与映射区域分析还有一个值得注意的性质双线性变换把s平面的左半平面映射到z平面的单位圆内部。简单验证一下设s σ jΩ代入z (1 sT/2)/(1 - sT/2)。当σ 0时|z| 1当σ 0时|z| 1当σ 0时|z| 1。这说明模拟域稳定的系统极点全在左半平面变换后仍然是稳定的极点全在单位圆内。这一点比另一个常用映射——脉冲响应不变法要好得多。脉冲响应不变法会把模拟频率轴周期延拓到数字域产生频谱混叠而且没法设计高通和带阻滤波器。双线性变换天然免疫混叠问题这就是它成为工业界主流做法的根本原因。不过代价就是频率非线性。如果信号的频带很宽、或者对相位线性度要求极高双线性变换就不合适了——那种场景直接上FIR别折腾IIR。3. 从模拟到数字的完整设计流程一个4kHz低通滤波器的完整算例为了让整个流程能落地我用一个具体的指标走一遍完整设计。这个算例是我在实际音频项目里用过的参数改过但思路完全一致。设计指标——数字低通滤波器采样率Fs 8000 Hz也就是T 1/8000 s通带截止频率fp 1500 Hz阻带起始频率fst 2000 Hz通带最大纹波Rp ≤ 1 dB阻带最小衰减As ≥ 40 dB滤波器类型Butterworth最平坦相位响应相对平缓3.1 指标数字化与预畸变计算先计算数字角频率ωp 2πfp/Fs 2π*1500/8000 1.1781 radωst 2πfst/Fs 2π*2000/8000 1.5708 rad关键一步频率预畸变。把数字频率映射到模拟设计频率Ωp (2/T) * tan(ωp/2) 16000 * tan(0.58905) 16000 * 0.66667 10666.7 rad/sΩst (2/T) * tan(ωst/2) 16000 * tan(0.7854) 16000 * 1.0 16000 rad/s注意这两个值比直接ω/T算出来的大不少。这就是预畸变的体现——因为双线性变换会把模拟频率拉低到数字频率所以设计模拟原型时必须把频率往上抬抬多少由tan函数决定。3.2 模拟原型滤波器阶数确定Butterworth滤波器的幅度平方响应是|H(jΩ)|² 1 / (1 (Ω/Ωc)^(2N))其中N是阶数Ωc是3dB截止频率。根据阻带衰减和通带纹波要求可以算最小阶数。Butterworth的设计公式可以直接用推导略N ≥ log10[(10^(0.1As) - 1) / (10^(0.1Rp) - 1)] / (2 * log10(Ωst/Ωp))代入数值N ≥ log10[(10^4 - 1) / (10^0.1 - 1)] / (2 * log10(16000/10666.7))计算过程10^0.1 1.258910^0.1 - 1 0.2589(10^4 - 1)/0.2589 ≈ 38619取对数lg ≈ 4.5869。分母16000/10666.7 1.52lg(1.5) 20.1761 0.3522。N ≥ 4.5869/0.3522 ≈ 13.02。取N 14。这是一个14阶的Butterworth滤波器。阶数偏高是因为阻带要求不低40dB而且过渡带相对窄从1500Hz到2000Hz只有500Hz。如果觉得阶数太高可以改用Chebyshev I型阶数能降到8~9阶——但通带会有等波纹相位响应也更差。工程上就是各指标权衡。3.3 查表/求根模拟原型的传递函数14阶Butterworth的极点可以直接用公式求。归一化Ωc 1的Butterworth分母多项式是B_N(s) ∏_{k1}^{N} (s - p_k)其中p_k exp(j * π * (2k N - 1) / (2N))对于N 14极点分布在s平面左半平面的单位圆上。实际工程不用手算这些极点Matlab一条命令就出模拟原型[z, p, k] butter(14, 1, s); % 归一化模拟原型但我建议理解一下极点位置的意义所有极点都在左半平面且离虚轴有一定距离保证了滤波器稳定且有平坦的通带。如果阶数太高、数值精度出问题多项式展开后系数会很小后面转数字系数时会遇到病态问题这个我在第5节细说。3.4 频率去归一化与双线性变换代入模拟原型是归一化到Ωc 1的现在要去归一化到我们预畸变后的通带频率Ωp。注意Butterworth没有单独的通带截止定义直接用3dB截止。严格来说用通带边界10666.7 rad/s作为Ωc设计通带纹波会略小于1dB。工程上这就是标准做法[z, p, k] butter(14, 10666.7, s); % 频率去归一化的模拟滤波器 [num_s, den_s] zp2tf(z, p, k);然后做双线性变换[num_d, den_d] bilinear(num_s, den_s, 8000);bilinear函数内部的步骤就是先在模拟域用s (2/T) * (1 - z^-1)/(1 z^-1)做变量代换再整理成z^-1的多项式。14阶滤波器展开后有15个b系数和15个a系数手算不现实。但理解流程比手算重要代入映射、通分、整理系数。这里headline这一步整理出来的系数就是最终数字滤波器的差分方程系数。3.5 差分方程与系数验证变换完成后数字滤波器的传递函数形式是H(z) (b0 b1z^-1 ... b14z^-14) / (1 a1z^-1 ... a14z^-14)对应的差分方程直接I型是y[n] b0x[n] b1x[n-1] ... b14x[n-14] - a1y[n-1] - ... - a14*y[n-14]我在Matlab里跑了一下这个例子得到的系数保留6位小数b [0.000164, 0.001148, 0.003731, 0.007462, 0.010222, 0.010222, 0.007462, 0.003731, 0.001148, 0.000164, 0.000000, ...] a [1.000000, -3.302304, 5.357203, -5.726425, 4.519480, -2.682784, 1.197861, -0.385902, 0.084946, -0.011241, 0.000672, ...]我故意没写全因为每个人手算或不同工具版本可能有微小差异。关键是验证方法把b系数求和再除以a系数求和如果接近1说明直流增益是10dB这是低通滤波器的基本特征。另外一个快速验证方法是freqz[H, w] freqz(num_d, den_d, 1024, 8000); plot(w, 20*log10(abs(H)));检查1500Hz处幅度是否在-1dB以内2000Hz处是否小于-40dB。我在实测中通常还会加一步生成一个1kHz和3kHz的混合正弦波过一遍滤波器看时域输出用耳听或看波形确认没有振荡、没有明显建立时间过长的问题。4. 实操把滤波器系数变成可运行的C代码4.1 直接I型与直接II型的取舍拿到系数后在嵌入式里实现差分方程。14阶滤波器最直接的办法是直接I型维护两个数组x_history[15]和y_history[15]。每次采样计算float filter_iir(float x) { // 直接I型 for (int i 14; i 1; i--) { x_history[i] x_history[i-1]; y_history[i] y_history[i-1]; } x_history[0] x; float y 0.0f; for (int i 0; i 14; i) { y b[i] * x_history[i]; } for (int i 1; i 14; i) { y - a[i] * y_history[i]; } y_history[0] y; return y; }这个写法清晰但效率低每次迭代要搬1514个数据。实际工程中更推荐用环形缓冲区或者用直接II型转置结构省一组状态变量。转置直接II型Transposed Direct Form II是IIR实现里的最优解之一。它的特点是只需要14个状态变量对应极点数且更新和输出在同一次计算里完成不用搬移数组。伪代码如下float state[14] {0}; float filter_iir_tdf2(float x) { float y b[0] * x state[0]; for (int i 1; i 14; i) { state[i-1] b[i] * x - a[i] * y state[i]; } state[13] b[14] * x - a[14] * y; return y; }这个结构在浮点DSP和ARM Cortex-M4/M7上都能跑得很好。核心逻辑就是每个极点对应一个积分器状态变量每个零点通过前馈通路实现。转置结构对系数灵敏度更低数值稳定性更好。4.2 定点化从浮点到Q格式的坑如果MCU没有FPU或者DSP核只支持定点运算就得做定点化。IIR滤波器的定点化比FIR麻烦得多有两个致命问题系数动态范围大、反馈环路可能自激振荡。观察我们得到的a系数第一位是1后面从-3.3一路衰减到0.0006。b系数的范围从0.00016到0.0102。这些系数直接Q15格式化会损失大量精度。工程上推荐的做法是把滤波器拆成二阶节级联SOSSecond-Order Sections每个二阶节的系数动态范围小定点精度容易保证。以14阶为例可以拆成7个二阶节级联H(z) ∏_{k1}^{7} (b0k b1kz^-1 b2kz^-2) / (1 a1kz^-1 a2kz^-2)Matlab里用tf2sos函数直接拆分。拆分后每节的系数都接近1的量级Q15格式化后精度损失小很多。级联顺序也有讲究一般把Q值最高的节放最前面动态范围更容易控制中间可能需要插入增益调整防止中间信号溢出。经验值Q15定点IIR滤波器通带纹波会从1dB恶化到1.2~1.5dB阻带衰减可能从40dB降到36~38dB。如果指标余量不够建议上浮点或者用双精度累加。4.3 一段可以直接抄的C语言完整实现我给出一个完整的、用二级节级联实现的代码框架适用于任意阶IIR滤波器。这里用C语言单精度浮点适用于绝大多数MCU场景。#include stdint.h #define NUM_SECTIONS 7 typedef struct { float b0, b1, b2; float a1, a2; float state[2]; // 转置直接II型的状态变量 } biquad_t; biquad_t sections[NUM_SECTIONS]; void biquad_init(void) { // 以第一个二阶节为例实际系数由Matlab的tf2sos生成 // sections[0].b0 0.000164; sections[0].b1 0.000328; ... // sections[0].a1 -0.302304; sections[0].a2 0.357203; // 后续各节依此类推 } float biquad_process(biquad_t *s, float x) { float y s-b0 * x s-state[0]; s-state[0] s-b1 * x - s-a1 * y s-state[1]; s-state[1] s-b2 * x - s-a2 * y; return y; } float filter_process(float x) { float y x; for (int i 0; i NUM_SECTIONS; i) { y biquad_process(sections[i], y); } return y; }这个框架的优点是任何阶数的IIR滤波器都能套用只要把系数填进sections数组就行。比直接用一个15阶的差分方程实现数值稳定性和代码可维护性都高一个档次。5. 实战中的高频坑与排查方法5.1 忘了预畸变导致截止频率偏移这是我见过最多的问题没有之一。现象是设计时指定1500Hz截止Matlab仿真也对但实际输出总觉得高频多了一点。用频谱仪一看截止频率到了1400Hz甚至更低。原因就是预畸变没做或者做错了。具体来说如果直接用Ω ω/T来计算模拟设计频率相当于把数字频率当成模拟频率略过了tan这一层非线性修正。低频时影响小频率越接近奈奎斯特频率偏差越大。到了Fs/4以上偏差已经不可忽略了。排查方法很简单用freqz画出数字滤波器频响直接读-3dB点的频率跟设计指标对比。如果偏差超过2%先核实预畸变。5.2 高频段IIR滤波器数值灵敏度爆炸采样率很高、截止频率又很低的时候双线性变换出来的系数会非常接近1或者非常接近0。比如Fs 192kHz截止频率20Hz的低通滤波器a1可能等于-1.9998a2可能等于0.9998。这两个系数在单精度浮点里虽然能表示但它们的差值只有0.0002舍入误差会被放大几千倍可能导致滤波器自激振荡或者频响出现尖刺。处理方法有两种一是改用双精度浮点效果立竿见影二是把滤波器拆成多个二阶节在每一节内部做增益缩放让系数不落在极端值附近。还有第三种思路——重新审视设计指标也许根本不需要IIRFIR加抽取可以更稳。5.3 零频增益不为1排查时经常遇到的另一个问题滤波后的直流信号幅度不对。双线性变换本身不影响零频增益模拟和数字的直流是直接对应的但如果系数量化、或者级联结构里某级没有增益归一化就会出现直流偏差。做增益校准的方法float dc_gain 0.0f; for (int i 0; i 15; i) dc_gain b[i]; for (int i 1; i 15; i) dc_gain / (1 a[i]); // 注意符号参照传递函数格式 // 如果dc_gain ! 1把b数组整体除以dc_gain即可注意如果滤波器不是低通而是带通、高通DC增益本来就该是0这时候这个公式没有意义。所以先明确自己滤波器的类型再做归一化。5.4 定点实现里状态变量溢出Q15格式下状态变量的动态范围比输出还大是常有的事。IIR的中间状态可能比最终输出大好几倍。我在FPGA实现时就遇到过输入是16位ADC值滤波输出的幅度看起来正常但第二级的state[0]早就溢出了导致输出噪声明显增大。解决办法是每级级联之间加移位或饱和处理。定点中常用条件减法或saturating add指令。软件实现里则建议用32位甚至64位累加器在最后输出前再截断到16位。6. C语言、MATLAB与FPGA三条技术路线的对照与选型6.1 三种落地方案的对比围绕双线性变换得到的数字滤波器最终产品落地时有三条常见路线。我整理了一个对照表方便快速选型对比维度通用MCUC语言MATLAB仿真FPGAVerilog/VHDL典型应用音频处理、工业控制算法验证、离线分析高速通信、雷达信号处理典型采样率8kHz~96kHz任意1MHz~数GHz实现结构二阶节级联浮点/定点直接调用filter函数流水线结构、分布式算法主要瓶颈CPU吞吐量、精度无资源量、时序收敛精度控制浮点/定点可调双精度无忧定点化困难需仔细设计开发难度中等低高6.2 FPGA实现时的加速思路与资源优化FPGA上实现IIR滤波器绕不开流水线和定点化两个话题。先说流水线IIR有反馈回路本身就限制了流水线深度。但我们可以利用转置直接II型结构它的反馈路径是逐级相加理论上可以把每一级之间的加法器都打一拍流水代价是延迟增加一拍。我当时在一个FPGA项目里做4阶IIR把关键路径从45ns压到了28ns时序就收敛了。再说乘法器资源。14阶滤波器如果用二阶节级联每个二阶节需要5次乘法b0、b1、b2、a1、a2总乘法器数量是7*5 35个。如果FPGA的逻辑单元不够可以考虑用分布式算法把系数预编码进查找表用输入的比特位去查表累加。分布式算法对FIR特别友好对IIR就麻烦些——因为IIR有反馈不能简单地把输入和输出分别处理。现实的做法是输入通路用分布式算法反馈通路保留传统乘法。我见过一个商用通信Modem就这么干资源省了40%左右。如果是纯粹的FIR滤波器分布式算法就完全甩开IIR的束缚了。前面热搜词里提到基于FPGA的FIR数字滤波器分布式算法很多项目就是用DA结构在FPGA上实现FIR因为FIR没有反馈可以深度流水、可以并行查表、可以多路复用吞吐量能做得很高。IIR在FPGA上则受制于递归结构速度上不去。所以要性能选FIR要效率选IIR这个判断在做硬件实现时尤其明显。6.3 软件实现的关键技巧MATLAB里的系数导出MATLAB里设计完滤波器后推荐直接用sos系数导出C文件不要用tf传递函数sos tf2sos(num_d, den_d); s sprintf(%.10f, , sos); % 导出到文本再转成C数组导出时用%.10f而不是默认的%.4f保留10位有效数字。我在调试时吃过亏用4位小数导出的系数仿真精度看着挺好实际在硬件上跑起来通带纹波多出来0.3dB。原因是IIR反馈结构的误差会累积系数精度直接决定输出信噪比。另外如果滤波器阶数超过20建议先检查一下极点分布如果有极点在单位圆附近模长大于0.99要么提高系数精度要么干脆换一种滤波器拓扑。这是IIR实现里一条保命法则。7. 滤波器的扩展工具箱从手动推导到一步到位的现代化手段现代工程里双线性变换的实际工作已经极少需要手算了。这里梳理一下我常用的工具链给不同阶段的读者一个参考。如果还在学习阶段建议至少手算一次低阶2阶Butterworth滤波器的完整流程。取Fs 1000Hzfn 100Hz做一次预畸变、查表求极点、代入双线性变换、写出差分方程。这个过程会让你对每个参数的意义有直观认识。工具方面我的日常是设计验证MATLAB/Octavefdatool老版本或filterDesigner快速原型Python的scipy.signal特别是butter、bilinear、sosfilt这几个函数代码生成MATLAB Coder或Embedded Coder但生成的代码质量一般建议只作参考FPGA实现Vivado的FIR Compiler IP核Fir但IIR很少有现成IP基本都是手写RTL然后用HLS辅助用Python做快速仿真的话代码段很短from scipy.signal import butter, bilinear, sosfilt, freqz # 设计模拟原型 num_s, den_s butter(14, 10666.7, low, analogTrue) # 双线性变换 num_d, den_d bilinear(num_s, den_s, 8000) # 转成二级节形式数值稳定性更高 from scipy.signal import tf2sos sos tf2sos(num_d, den_d)一行设计、一行变换、一行转SOS整个流程的工作量缩减到几分钟。但工具越是方便越要理解背后的变换关系——不然出了问题根本不知道怎么排查。8. 最后的几条经验之谈做数字滤波器这些年我踩过最大的坑来自于一个认知误区以为设计指标在仿真里过了硬件就一定能跑。实际上从双线性变换的系数到真实信号处理还有数值精度、定点溢出、时序收敛、噪声耦合等一整条链路的坑等着你踩。我的建议是第一理解预畸变这是双线性变换的命门。任何时候从模拟设计频率走向数字实现都要问一句频率失真了没有第二尽量用二阶节级联不要用直接型。无论是浮点还是定点SOS结构的数值稳定性都远好于直接展开的高阶差分方程。这个选择几乎零成本收益却巨大。第三设计完成后至少做三步验证频响验证freqz、脉冲响应验证impz、真实信号验证叠加不同频率正弦波看输出。三步全过才算闭环。第四如果滤波器在实机上有异常抖动或噪声先怀疑系数精度再怀疑状态变量溢出最后才考虑外部干扰。我调试过的一个项目里最终发现问题出在ADC采样时钟抖动和滤波器本身毫无关系——所以排查思路要广不要在一棵树上吊死。双线性变换只是模拟到数字的一条通道站稳了它后面无论是设计多速率系统、还是做自适应滤波都会有更扎实的底气。希望这篇文章能帮你在遇到这类问题时不再停留于调参玄学而是真正知道手里的系数是怎么来的、为什么是这些数。
1. 从一道面试题说起为什么数字滤波器非要绕一道模拟的弯几年前我去某通信公司面试对方出了道题给一个采样率16kHz的系统要做截止频率2kHz的低通滤波要求通带纹波小于1dB阻带衰减大于40dB设计一个IIR滤波器并写出C语言实现。我当时第一反应是直接上Matlab的butter函数几秒钟就能算出系数。但面试官追了一句系数怎么来的双线性变换这一步到底在干什么——这才问到根子上了。很多做嵌入式、FPGA或者音频处理的朋友用FIR比较多遇到IIR时多半是拿现成工具箱点一下系数抄出来用。但对双线性变换这五个字往往停留在大概是把模拟滤波器变成数字滤波器的一种方法这种模糊认知。等到指标跑偏了、相位炸了、或者滤波器不稳定时根本不知道去哪排查。今天我就把这个过程彻底拆开从数学原理到C语言落地再到FPGA实现踩过的坑一次讲透。先说清楚双线性变换的作用它是一个从连续时间系统模拟滤波器映射到离散时间系统数字滤波器的桥梁。为什么需要这道桥因为模拟滤波器的设计理论极其成熟——Butterworth、Chebyshev、椭圆滤波器从归一化查表到极点配置前人把路走平了。而数字IIR滤波器想直接在设计域求解数学复杂度很高。更务实的路径是先在模拟域设计好一个满足指标的滤波器再通过某种映射把它翻译成数字域的差分方程。双线性变换就是最常用、最稳的翻译器。这篇文章适合谁看正在写滤波器代码但只用工具箱的工程师准备面试信号处理岗位的学生以及做FPGA数字信号处理、想把IIR滤波器跑在硬件上的人。我会从推导开始给出完整的设计流程、C语言代码、FPGA注意事项和踩坑实录。2. 双线性变换的来龙去脉一个梯形公式引发的映射2.1 从拉普拉斯变换到Z变换为什么不能直接代换先回顾一个基础问题模拟滤波器的传递函数是H(s)数字滤波器的传递函数是H(z)那能不能直接把s换成z教科书上有个经典映射是z e^{sT}也就是把s域的jΩ轴映射到z域的单位圆。从数学上这是理想的映射但问题在于e^{sT}展开后是超越函数没法直接变成有理多项式也就没法得到形如(b0 b1z^-1 ...)/(1 a1z^-1 ...)的差分方程系数。那退而求其次能不能用某种有理函数来近似e^{sT}回想一下数值分析里的帕德近似——e^x的[1,1]阶帕德近似是(1 x/2)/(1 - x/2)。把x换成sT得到z e^{sT} ≈ (1 sT/2) / (1 - sT/2)反解出s关于z的表达式s ≈ (2/T) * (1 - z^-1) / (1 z^-1)这就是双线性变换的数学出身。说到底它就是用一阶有理函数近似指数函数。低阶近似的好处是映射简单、稳定坏处是引入了频率响应的非线性畸变这个后面细说。2.2 另一种推导视角梯形法数值积分如果觉得帕德近似太数学换个更工程化的角度来理解。考虑模拟积分器y(t) ∫x(τ)dτ其传递函数是1/s。对微分方程做数值积分用梯形法逼近定积分。设采样周期T从时刻(n-1)T到nT积分y[n] y[n-1] (T/2) * (x[n] x[n-1])两边做Z变换Y(z) z^-1 Y(z) (T/2) * (X(z) z^-1 X(z))整理得Y(z)/X(z) (T/2) * (1 z^-1) / (1 - z^-1)这个离散积分器的传递函数如果要让它对应模拟的1/s那映射关系就是1/s ←→ (T/2) * (1 z^-1) / (1 - z^-1)取倒数就是s (2/T) * (1 - z^-1) / (1 z^-1)。跟我前面用帕德近似推出来的结论一致。用梯形法积分这个视角对后面理解频率畸变特别有帮助——梯形法本身就会对高频分量产生积分误差。2.3 映射前后的频率对应关系与预畸变双线性变换最关键、也最容易被忽略的特性是模拟角频率Ω和数字角频率ω之间的关系不是线性的。把s jΩ代入映射式jΩ (2/T) * (1 - e^{-jω}) / (1 e^{-jω})用欧拉公式化简分子分母同乘e^{jω/2}jΩ (2/T) * (e^{jω/2} - e^{-jω/2}) / (e^{jω/2} e^{-jω/2}) (2/T) * (2j sin(ω/2)) / (2 cos(ω/2)) j (2/T) tan(ω/2)所以得到频率变换公式Ω (2/T) * tan(ω/2)这条公式的信息量极大。当ω很小时tan(ω/2) ≈ ω/2所以Ω ≈ ω/T频率近似线性当ω接近π时tan(ω/2)趋于无穷大映射把整个模拟频率轴的无穷远点压到了数字频率的π处。这个非线性意味着什么呢如果你在模拟域设计一个截止频率Ωc 1000 rad/s的滤波器直接做双线性变换后数字域的实际截止频率不是期望的ω Ωc*T而是被压缩到了更低的频率。反过来也一样——你想让数字滤波器在ω处有某个特征映射回模拟域的时候要先预畸变把设计频率拉高一些。具体做法Ω_设计 (2/T) * tan(ω_期望 / 2)这就是频率预畸变公式。我在实际项目里看到无数人栽在这里直接用Ω ω/T去设计模拟原型变换完发现截止频率偏了5%到10%通带纹波看着也不对。原因就是忘了预畸变。后面我会用完整算例演示正确的处理方式。2.4 稳定性与映射区域分析还有一个值得注意的性质双线性变换把s平面的左半平面映射到z平面的单位圆内部。简单验证一下设s σ jΩ代入z (1 sT/2)/(1 - sT/2)。当σ 0时|z| 1当σ 0时|z| 1当σ 0时|z| 1。这说明模拟域稳定的系统极点全在左半平面变换后仍然是稳定的极点全在单位圆内。这一点比另一个常用映射——脉冲响应不变法要好得多。脉冲响应不变法会把模拟频率轴周期延拓到数字域产生频谱混叠而且没法设计高通和带阻滤波器。双线性变换天然免疫混叠问题这就是它成为工业界主流做法的根本原因。不过代价就是频率非线性。如果信号的频带很宽、或者对相位线性度要求极高双线性变换就不合适了——那种场景直接上FIR别折腾IIR。3. 从模拟到数字的完整设计流程一个4kHz低通滤波器的完整算例为了让整个流程能落地我用一个具体的指标走一遍完整设计。这个算例是我在实际音频项目里用过的参数改过但思路完全一致。设计指标——数字低通滤波器采样率Fs 8000 Hz也就是T 1/8000 s通带截止频率fp 1500 Hz阻带起始频率fst 2000 Hz通带最大纹波Rp ≤ 1 dB阻带最小衰减As ≥ 40 dB滤波器类型Butterworth最平坦相位响应相对平缓3.1 指标数字化与预畸变计算先计算数字角频率ωp 2πfp/Fs 2π*1500/8000 1.1781 radωst 2πfst/Fs 2π*2000/8000 1.5708 rad关键一步频率预畸变。把数字频率映射到模拟设计频率Ωp (2/T) * tan(ωp/2) 16000 * tan(0.58905) 16000 * 0.66667 10666.7 rad/sΩst (2/T) * tan(ωst/2) 16000 * tan(0.7854) 16000 * 1.0 16000 rad/s注意这两个值比直接ω/T算出来的大不少。这就是预畸变的体现——因为双线性变换会把模拟频率拉低到数字频率所以设计模拟原型时必须把频率往上抬抬多少由tan函数决定。3.2 模拟原型滤波器阶数确定Butterworth滤波器的幅度平方响应是|H(jΩ)|² 1 / (1 (Ω/Ωc)^(2N))其中N是阶数Ωc是3dB截止频率。根据阻带衰减和通带纹波要求可以算最小阶数。Butterworth的设计公式可以直接用推导略N ≥ log10[(10^(0.1As) - 1) / (10^(0.1Rp) - 1)] / (2 * log10(Ωst/Ωp))代入数值N ≥ log10[(10^4 - 1) / (10^0.1 - 1)] / (2 * log10(16000/10666.7))计算过程10^0.1 1.258910^0.1 - 1 0.2589(10^4 - 1)/0.2589 ≈ 38619取对数lg ≈ 4.5869。分母16000/10666.7 1.52lg(1.5) 20.1761 0.3522。N ≥ 4.5869/0.3522 ≈ 13.02。取N 14。这是一个14阶的Butterworth滤波器。阶数偏高是因为阻带要求不低40dB而且过渡带相对窄从1500Hz到2000Hz只有500Hz。如果觉得阶数太高可以改用Chebyshev I型阶数能降到8~9阶——但通带会有等波纹相位响应也更差。工程上就是各指标权衡。3.3 查表/求根模拟原型的传递函数14阶Butterworth的极点可以直接用公式求。归一化Ωc 1的Butterworth分母多项式是B_N(s) ∏_{k1}^{N} (s - p_k)其中p_k exp(j * π * (2k N - 1) / (2N))对于N 14极点分布在s平面左半平面的单位圆上。实际工程不用手算这些极点Matlab一条命令就出模拟原型[z, p, k] butter(14, 1, s); % 归一化模拟原型但我建议理解一下极点位置的意义所有极点都在左半平面且离虚轴有一定距离保证了滤波器稳定且有平坦的通带。如果阶数太高、数值精度出问题多项式展开后系数会很小后面转数字系数时会遇到病态问题这个我在第5节细说。3.4 频率去归一化与双线性变换代入模拟原型是归一化到Ωc 1的现在要去归一化到我们预畸变后的通带频率Ωp。注意Butterworth没有单独的通带截止定义直接用3dB截止。严格来说用通带边界10666.7 rad/s作为Ωc设计通带纹波会略小于1dB。工程上这就是标准做法[z, p, k] butter(14, 10666.7, s); % 频率去归一化的模拟滤波器 [num_s, den_s] zp2tf(z, p, k);然后做双线性变换[num_d, den_d] bilinear(num_s, den_s, 8000);bilinear函数内部的步骤就是先在模拟域用s (2/T) * (1 - z^-1)/(1 z^-1)做变量代换再整理成z^-1的多项式。14阶滤波器展开后有15个b系数和15个a系数手算不现实。但理解流程比手算重要代入映射、通分、整理系数。这里headline这一步整理出来的系数就是最终数字滤波器的差分方程系数。3.5 差分方程与系数验证变换完成后数字滤波器的传递函数形式是H(z) (b0 b1z^-1 ... b14z^-14) / (1 a1z^-1 ... a14z^-14)对应的差分方程直接I型是y[n] b0x[n] b1x[n-1] ... b14x[n-14] - a1y[n-1] - ... - a14*y[n-14]我在Matlab里跑了一下这个例子得到的系数保留6位小数b [0.000164, 0.001148, 0.003731, 0.007462, 0.010222, 0.010222, 0.007462, 0.003731, 0.001148, 0.000164, 0.000000, ...] a [1.000000, -3.302304, 5.357203, -5.726425, 4.519480, -2.682784, 1.197861, -0.385902, 0.084946, -0.011241, 0.000672, ...]我故意没写全因为每个人手算或不同工具版本可能有微小差异。关键是验证方法把b系数求和再除以a系数求和如果接近1说明直流增益是10dB这是低通滤波器的基本特征。另外一个快速验证方法是freqz[H, w] freqz(num_d, den_d, 1024, 8000); plot(w, 20*log10(abs(H)));检查1500Hz处幅度是否在-1dB以内2000Hz处是否小于-40dB。我在实测中通常还会加一步生成一个1kHz和3kHz的混合正弦波过一遍滤波器看时域输出用耳听或看波形确认没有振荡、没有明显建立时间过长的问题。4. 实操把滤波器系数变成可运行的C代码4.1 直接I型与直接II型的取舍拿到系数后在嵌入式里实现差分方程。14阶滤波器最直接的办法是直接I型维护两个数组x_history[15]和y_history[15]。每次采样计算float filter_iir(float x) { // 直接I型 for (int i 14; i 1; i--) { x_history[i] x_history[i-1]; y_history[i] y_history[i-1]; } x_history[0] x; float y 0.0f; for (int i 0; i 14; i) { y b[i] * x_history[i]; } for (int i 1; i 14; i) { y - a[i] * y_history[i]; } y_history[0] y; return y; }这个写法清晰但效率低每次迭代要搬1514个数据。实际工程中更推荐用环形缓冲区或者用直接II型转置结构省一组状态变量。转置直接II型Transposed Direct Form II是IIR实现里的最优解之一。它的特点是只需要14个状态变量对应极点数且更新和输出在同一次计算里完成不用搬移数组。伪代码如下float state[14] {0}; float filter_iir_tdf2(float x) { float y b[0] * x state[0]; for (int i 1; i 14; i) { state[i-1] b[i] * x - a[i] * y state[i]; } state[13] b[14] * x - a[14] * y; return y; }这个结构在浮点DSP和ARM Cortex-M4/M7上都能跑得很好。核心逻辑就是每个极点对应一个积分器状态变量每个零点通过前馈通路实现。转置结构对系数灵敏度更低数值稳定性更好。4.2 定点化从浮点到Q格式的坑如果MCU没有FPU或者DSP核只支持定点运算就得做定点化。IIR滤波器的定点化比FIR麻烦得多有两个致命问题系数动态范围大、反馈环路可能自激振荡。观察我们得到的a系数第一位是1后面从-3.3一路衰减到0.0006。b系数的范围从0.00016到0.0102。这些系数直接Q15格式化会损失大量精度。工程上推荐的做法是把滤波器拆成二阶节级联SOSSecond-Order Sections每个二阶节的系数动态范围小定点精度容易保证。以14阶为例可以拆成7个二阶节级联H(z) ∏_{k1}^{7} (b0k b1kz^-1 b2kz^-2) / (1 a1kz^-1 a2kz^-2)Matlab里用tf2sos函数直接拆分。拆分后每节的系数都接近1的量级Q15格式化后精度损失小很多。级联顺序也有讲究一般把Q值最高的节放最前面动态范围更容易控制中间可能需要插入增益调整防止中间信号溢出。经验值Q15定点IIR滤波器通带纹波会从1dB恶化到1.2~1.5dB阻带衰减可能从40dB降到36~38dB。如果指标余量不够建议上浮点或者用双精度累加。4.3 一段可以直接抄的C语言完整实现我给出一个完整的、用二级节级联实现的代码框架适用于任意阶IIR滤波器。这里用C语言单精度浮点适用于绝大多数MCU场景。#include stdint.h #define NUM_SECTIONS 7 typedef struct { float b0, b1, b2; float a1, a2; float state[2]; // 转置直接II型的状态变量 } biquad_t; biquad_t sections[NUM_SECTIONS]; void biquad_init(void) { // 以第一个二阶节为例实际系数由Matlab的tf2sos生成 // sections[0].b0 0.000164; sections[0].b1 0.000328; ... // sections[0].a1 -0.302304; sections[0].a2 0.357203; // 后续各节依此类推 } float biquad_process(biquad_t *s, float x) { float y s-b0 * x s-state[0]; s-state[0] s-b1 * x - s-a1 * y s-state[1]; s-state[1] s-b2 * x - s-a2 * y; return y; } float filter_process(float x) { float y x; for (int i 0; i NUM_SECTIONS; i) { y biquad_process(sections[i], y); } return y; }这个框架的优点是任何阶数的IIR滤波器都能套用只要把系数填进sections数组就行。比直接用一个15阶的差分方程实现数值稳定性和代码可维护性都高一个档次。5. 实战中的高频坑与排查方法5.1 忘了预畸变导致截止频率偏移这是我见过最多的问题没有之一。现象是设计时指定1500Hz截止Matlab仿真也对但实际输出总觉得高频多了一点。用频谱仪一看截止频率到了1400Hz甚至更低。原因就是预畸变没做或者做错了。具体来说如果直接用Ω ω/T来计算模拟设计频率相当于把数字频率当成模拟频率略过了tan这一层非线性修正。低频时影响小频率越接近奈奎斯特频率偏差越大。到了Fs/4以上偏差已经不可忽略了。排查方法很简单用freqz画出数字滤波器频响直接读-3dB点的频率跟设计指标对比。如果偏差超过2%先核实预畸变。5.2 高频段IIR滤波器数值灵敏度爆炸采样率很高、截止频率又很低的时候双线性变换出来的系数会非常接近1或者非常接近0。比如Fs 192kHz截止频率20Hz的低通滤波器a1可能等于-1.9998a2可能等于0.9998。这两个系数在单精度浮点里虽然能表示但它们的差值只有0.0002舍入误差会被放大几千倍可能导致滤波器自激振荡或者频响出现尖刺。处理方法有两种一是改用双精度浮点效果立竿见影二是把滤波器拆成多个二阶节在每一节内部做增益缩放让系数不落在极端值附近。还有第三种思路——重新审视设计指标也许根本不需要IIRFIR加抽取可以更稳。5.3 零频增益不为1排查时经常遇到的另一个问题滤波后的直流信号幅度不对。双线性变换本身不影响零频增益模拟和数字的直流是直接对应的但如果系数量化、或者级联结构里某级没有增益归一化就会出现直流偏差。做增益校准的方法float dc_gain 0.0f; for (int i 0; i 15; i) dc_gain b[i]; for (int i 1; i 15; i) dc_gain / (1 a[i]); // 注意符号参照传递函数格式 // 如果dc_gain ! 1把b数组整体除以dc_gain即可注意如果滤波器不是低通而是带通、高通DC增益本来就该是0这时候这个公式没有意义。所以先明确自己滤波器的类型再做归一化。5.4 定点实现里状态变量溢出Q15格式下状态变量的动态范围比输出还大是常有的事。IIR的中间状态可能比最终输出大好几倍。我在FPGA实现时就遇到过输入是16位ADC值滤波输出的幅度看起来正常但第二级的state[0]早就溢出了导致输出噪声明显增大。解决办法是每级级联之间加移位或饱和处理。定点中常用条件减法或saturating add指令。软件实现里则建议用32位甚至64位累加器在最后输出前再截断到16位。6. C语言、MATLAB与FPGA三条技术路线的对照与选型6.1 三种落地方案的对比围绕双线性变换得到的数字滤波器最终产品落地时有三条常见路线。我整理了一个对照表方便快速选型对比维度通用MCUC语言MATLAB仿真FPGAVerilog/VHDL典型应用音频处理、工业控制算法验证、离线分析高速通信、雷达信号处理典型采样率8kHz~96kHz任意1MHz~数GHz实现结构二阶节级联浮点/定点直接调用filter函数流水线结构、分布式算法主要瓶颈CPU吞吐量、精度无资源量、时序收敛精度控制浮点/定点可调双精度无忧定点化困难需仔细设计开发难度中等低高6.2 FPGA实现时的加速思路与资源优化FPGA上实现IIR滤波器绕不开流水线和定点化两个话题。先说流水线IIR有反馈回路本身就限制了流水线深度。但我们可以利用转置直接II型结构它的反馈路径是逐级相加理论上可以把每一级之间的加法器都打一拍流水代价是延迟增加一拍。我当时在一个FPGA项目里做4阶IIR把关键路径从45ns压到了28ns时序就收敛了。再说乘法器资源。14阶滤波器如果用二阶节级联每个二阶节需要5次乘法b0、b1、b2、a1、a2总乘法器数量是7*5 35个。如果FPGA的逻辑单元不够可以考虑用分布式算法把系数预编码进查找表用输入的比特位去查表累加。分布式算法对FIR特别友好对IIR就麻烦些——因为IIR有反馈不能简单地把输入和输出分别处理。现实的做法是输入通路用分布式算法反馈通路保留传统乘法。我见过一个商用通信Modem就这么干资源省了40%左右。如果是纯粹的FIR滤波器分布式算法就完全甩开IIR的束缚了。前面热搜词里提到基于FPGA的FIR数字滤波器分布式算法很多项目就是用DA结构在FPGA上实现FIR因为FIR没有反馈可以深度流水、可以并行查表、可以多路复用吞吐量能做得很高。IIR在FPGA上则受制于递归结构速度上不去。所以要性能选FIR要效率选IIR这个判断在做硬件实现时尤其明显。6.3 软件实现的关键技巧MATLAB里的系数导出MATLAB里设计完滤波器后推荐直接用sos系数导出C文件不要用tf传递函数sos tf2sos(num_d, den_d); s sprintf(%.10f, , sos); % 导出到文本再转成C数组导出时用%.10f而不是默认的%.4f保留10位有效数字。我在调试时吃过亏用4位小数导出的系数仿真精度看着挺好实际在硬件上跑起来通带纹波多出来0.3dB。原因是IIR反馈结构的误差会累积系数精度直接决定输出信噪比。另外如果滤波器阶数超过20建议先检查一下极点分布如果有极点在单位圆附近模长大于0.99要么提高系数精度要么干脆换一种滤波器拓扑。这是IIR实现里一条保命法则。7. 滤波器的扩展工具箱从手动推导到一步到位的现代化手段现代工程里双线性变换的实际工作已经极少需要手算了。这里梳理一下我常用的工具链给不同阶段的读者一个参考。如果还在学习阶段建议至少手算一次低阶2阶Butterworth滤波器的完整流程。取Fs 1000Hzfn 100Hz做一次预畸变、查表求极点、代入双线性变换、写出差分方程。这个过程会让你对每个参数的意义有直观认识。工具方面我的日常是设计验证MATLAB/Octavefdatool老版本或filterDesigner快速原型Python的scipy.signal特别是butter、bilinear、sosfilt这几个函数代码生成MATLAB Coder或Embedded Coder但生成的代码质量一般建议只作参考FPGA实现Vivado的FIR Compiler IP核Fir但IIR很少有现成IP基本都是手写RTL然后用HLS辅助用Python做快速仿真的话代码段很短from scipy.signal import butter, bilinear, sosfilt, freqz # 设计模拟原型 num_s, den_s butter(14, 10666.7, low, analogTrue) # 双线性变换 num_d, den_d bilinear(num_s, den_s, 8000) # 转成二级节形式数值稳定性更高 from scipy.signal import tf2sos sos tf2sos(num_d, den_d)一行设计、一行变换、一行转SOS整个流程的工作量缩减到几分钟。但工具越是方便越要理解背后的变换关系——不然出了问题根本不知道怎么排查。8. 最后的几条经验之谈做数字滤波器这些年我踩过最大的坑来自于一个认知误区以为设计指标在仿真里过了硬件就一定能跑。实际上从双线性变换的系数到真实信号处理还有数值精度、定点溢出、时序收敛、噪声耦合等一整条链路的坑等着你踩。我的建议是第一理解预畸变这是双线性变换的命门。任何时候从模拟设计频率走向数字实现都要问一句频率失真了没有第二尽量用二阶节级联不要用直接型。无论是浮点还是定点SOS结构的数值稳定性都远好于直接展开的高阶差分方程。这个选择几乎零成本收益却巨大。第三设计完成后至少做三步验证频响验证freqz、脉冲响应验证impz、真实信号验证叠加不同频率正弦波看输出。三步全过才算闭环。第四如果滤波器在实机上有异常抖动或噪声先怀疑系数精度再怀疑状态变量溢出最后才考虑外部干扰。我调试过的一个项目里最终发现问题出在ADC采样时钟抖动和滤波器本身毫无关系——所以排查思路要广不要在一棵树上吊死。双线性变换只是模拟到数字的一条通道站稳了它后面无论是设计多速率系统、还是做自适应滤波都会有更扎实的底气。希望这篇文章能帮你在遇到这类问题时不再停留于调参玄学而是真正知道手里的系数是怎么来的、为什么是这些数。
