网上搜“双线性变换”十篇有八篇上来就甩给你一个替换式s (2/T)(1 - z⁻¹)/(1 z⁻¹)。然后就是“代入即可、整理可得、最后得到”三连。公式谁都会抄问题是这个式子到底怎么来的为什么偏偏是它而不是 s (z-1)/T 这种看起来更直觉的替换模拟滤波器设计得好好的为什么非要绕这么大一圈去做数字化这篇文章我就把这笔账完整算一遍从微分方程讲到梯形积分再讲到频率畸变和预畸变最后用一个二阶巴特沃斯低通的实际案例把数字滤波器从模拟滤波器到C语言实现的全过程走一遍。适合正在学数字信号处理、被教材推导跳过细节卡住的人也适合想自己动手写滤波器却不知道系数怎么来的工程师。1. 从模拟到数字为什么不能直接把 s 换成某个 z 的表达式1.1 数字域和模拟域的根本差异模拟滤波器设计时我们通常工作在s域传递函数是 H(s)变量 s σ jΩ频率响应就是让 s 沿着虚轴走看 H(jΩ) 的值。模拟世界里的电感、电容、电阻它们的阻抗随频率变化构成了天然的滤波结构。设计出来的 H(s) 是连续函数输入输出都是连续时间信号。数字滤波器就不一样了它处理的是采样序列工作在z域传递函数是 H(z)变量 z 对应离散时间系统的单位延迟。z e^{sT} 这个关系把s平面虚轴上的点映射到z平面的单位圆上s平面左半平面映射到单位圆内部右半平面映射到单位圆外部。稳定性判据也随之变化模拟滤波器要求极点全部落在左半平面数字滤波器要求极点全部落在单位圆内。所以模拟滤波器数字化本质上是找一种映射关系把s域的设计成果搬到z域同时尽量保住频率选择性。这里最容易踩的坑是以为一切映射都可以随便把s替换成某个z的函数。实际上一种映射能否用至少要满足两个条件第一稳定的模拟系统映射后仍然是稳定的数字系统也就是s平面左半平面要映射到z平面单位圆内部第二虚轴 jΩ 要映射到单位圆否则频率响应的概念就变了。1.2 两条“看似合理”的直接替换路线都栽在这里先说前向欧拉法。微分的离散化最直觉就是前向差分dy/dt ≈ (y[n1] - y[n]) / T对应z变换后 y(z)(z - 1)/T ...因此 s (z - 1)/T。看着很顺但一旦替换进去s平面左半平面的点映射到哪了设 s σ jΩ那么 z 1 sT实部 Re(z) 1 σT。当 σ 0 时只要 |σ| 足够大z 的实部可以变成负数但模长呢|z|² (1 σT)² (ΩT)²。这个式子可以大于1也就是说一个原本稳定的模拟极点数字化后可能跑到单位圆外面去。这意味着什么你精心设计的模拟滤波器数字化之后居然不稳定了。我当年第一次验证这个结论时很震撼因为直观上差分近似不应该带来这么严重的后果但数学就是这么不讲情面。再说后向欧拉法。改成后向差分dy/dt ≈ (y[n] - y[n-1]) / T对应 s (z - 1) / (zT) (1 - z⁻¹)/T。这次稳定性没问题了s平面左半平面映射到z平面某个单位圆内部的区域数字系统一定稳定。但代价是模拟虚轴 jΩ 并没有映射到z平面单位圆上而是映射到单位圆内部一条曲线上。这导致数字滤波器的频率响应跟模拟原型的频率响应在形状上对不上失真比前向欧拉还难修正。把这两种情况摆在一起结论很明显直接替换不是不行而是总得在稳定性和频率保真之间做二选一两者都想要就得换思路。2. 双线性变换的来历它其实是一次梯形积分2.1 从传递函数回到微分方程要理解双线性变换为什么长成那样不能只盯着 s 和 z 的替换式得回到微分方程本身。以一阶低通系统为例模拟传递函数是H(s) a / (s a)它对应的微分方程是dy/dt a·y a·x其中 x(t) 是输入y(t) 是输出。这个方程物理含义很清晰输出变化率加上 a 倍输出等于 a 倍输入直流增益是1。现在要数字化就是对 dy/dt 做近似。前向欧拉和后向欧拉都是单点近似精度只有一阶而梯形法用的是两点的平均斜率精度高一阶稳定性和失真特性也更好。这就是双线性变换的核心思想来源用梯形法做数值积分。2.2 梯形法代入差分方程出炉梯形法把连续积分离散化的公式是y[n] y[n-1] (T/2)·(y[n] y[n-1])意思是这一段区间的积分面积用前后两点导数的平均值乘以区间长度来近似。把刚才的微分方程变形y a·x - a·y代入y[n] y[n-1] (T/2)·(a·x[n] - a·y[n] a·x[n-1] - a·y[n-1])两边整理。把含 y[n] 的项移到左边含 y[n-1] 的项移到右边(1 aT/2)·y[n] (1 - aT/2)·y[n-1] (aT/2)·(x[n] x[n-1])这就是一阶系统用双线性变换后得到的差分方程。注意 x[n-1] 和 y[n-1] 的交叉项这是梯形法的特征跟前向欧拉那种单边结构完全不同。到这里离散化已经做完了还没有出现 s 和 z 的替换式但系数结构已经体现出来数字滤波器的系数里有两个非平凡项一个来自当前时刻一个来自上一时刻这种“两边都取平均”的做法就是双线性变换名字的由来——它把连续时间积分变成了矩形加三角形面积几何上就是双线性两个线性函数夹出来的面积近似。2.3 把结果整理成 s 到 z 的替换式有了差分方程取z变换看看能不能整理出 H(z)。对差分方程两边做z变换(1 aT/2)·Y(z) (1 - aT/2)·z⁻¹·Y(z) (aT/2)·(1 z⁻¹)·X(z)移项后得到H(z) Y(z)/X(z) (aT/2)·(1 z⁻¹) / (1 aT/2 - (1 - aT/2)·z⁻¹)这个形式跟模拟的 H(s) a/(s a) 对不上但如果强行提出一个 s 的表达式让分母变成 1 (T/2)·s·(...) 的形式对比一下就发现s (2/T)·(1 - z⁻¹)/(1 z⁻¹)代入 H(s) a/(s a)得到的分式正好就是上面那个 H(z)。这才是双线性变换替换式的真正来源它不是凭空发明的一个映射而是“先对微分方程做梯形法离散化再反解出 s 和 z 的关系”。推导到这一步公式就不再是魔术而是一个操作流程模拟传递函数 H(s) 是微分方程的频域表达梯形法是离散化的数值手段两者结合就是双线性变换。我建议你在纸上自己走一遍一阶推导哪怕只有一次。很多教材把推导省略成“令 s 2/T·(1-z⁻¹)/(1z⁻¹)代入得”读者永远不知道这个“令”为什么成立。自己推一遍之后以后再看到任何双线性变换的变体例如带频率预畸变的版本都能一眼看穿。3. 频率畸变与补偿预畸变的设计流程3.1 映射关系到底压扁了什么双线性变换虽然解决了稳定性问题但付出了成本频率轴被非线性压缩了。把 s jΩ 代入替换式jΩ (2/T)·(1 - e^{-jω})/(1 e^{-jω})化简后得到频率映射关系ω_digital 2·arctan(Ω_analog·T/2)反过来Ω_analog (2/T)·tan(ω_digital/2)低频时tan(x) ≈ x所以 ω_digital ≈ Ω_analog·T频率是近似线性的。但频率一高就不对了。Ω 趋向无穷大时ω_digital 只趋向 π也就是数字频率的奈奎斯特频率 fs/2。换句话说整个模拟频率轴从 0 到无穷大被压缩到了数字频率 0 到 fs/2 的有限区间里。这个压缩效果在频响曲线上表现为原本设计好的模拟滤波器截止频率是 fc数字化之后实际截止频率会往低处偏。你以为是100Hz的截止可能做出来是97Hz极端情况可能偏到只剩原目标的三分之一。我见过不少人第一次遇到这个问题时以为是计算错误反复检查系数其实没有错就是频率轴压缩造成的。越靠近奈奎斯特频率压缩越严重。举个例子采样率 1kHz目标数字截止频率 450Hz 的滤波器如果不做预畸变直接用模拟截止频率 2π×450 去设计数字化后实际截止频率大约只有 304Hz。从450偏到304这个误差在滤波器的通带设计里是完全不可接受的。3.2 预畸变的标准操作解决办法是在设计模拟滤波器时先做一次频率预畸变。步骤是确定数字滤波器目标截止频率 ω_d单位是弧度/采样点通常 ω_d 2π·fc/fs。用公式 Ω_p (2/T)·tan(ω_d/2) 反算出模拟滤波器设计要用的截止频率。用这个 Ω_p 设计模拟滤波器 H(s)。再用双线性变换数字化。这里的逻辑是双线性变换会把模拟频率 Ω_p 压缩成数字频率 ω_d所以提前把设计频率抬高到 Ω_p数字化后就正好落在 ω_d 上。下表是采样率 1kHz、目标截止频率不同取值时预畸变前后的模拟设计频率对比目标数字截止频率 fc (Hz)数字角频率 ω_d (rad/采样)预畸变后模拟角频率 Ω_p (rad/s)不预畸变直接用 Ω ω_d/T 的后果 (数字化后实际截止 Hz)1000.6283649.8约 972001.25661426.5约 1853001.88502453.2约 2574002.51334108.7约 3144502.82746496.4约 304看到最后两行的差异了吗频率越高误差越大。做预畸变之后截止频率才能精确落在目标点上。这里必须澄清一个常见误解预畸变只是把设计频点对准了并不是让整条幅频曲线跟模拟原型完全重合。双线性变换的频率压缩是全局性的通带内的形状在接近奈奎斯特频率时仍然会被压缩变形只是截止频点这个关键节点被校正回来了。如果你的滤波器通带要求非常严格覆盖范围又很宽那么双线性变换可能不是最佳选择后续可以考虑更高阶的数值方法或者直接设计数字滤波器不走模拟原型这条路。4. 手把手设计一个二阶巴特沃斯低通数字系数全程推导4.1 设计指标与归一化原型我用一个具体例子把整个流程走一遍。设计目标采样率 fs 1000Hz截止频率 fc 100Hz二阶巴特沃斯低通滤波器。先算出数字角频率ω_d 2π·fc/fs 2π×100/1000 0.6283 rad采样间隔 T 1/fs 0.001s。预畸变后的模拟截止频率Ω_p (2/T)·tan(ω_d/2) 2000 × tan(0.31415) ≈ 2000 × 0.3249 649.8 rad/s二阶巴特沃斯低通模拟原型的归一化传递函数是H(s_n) 1 / (s_n² √2·s_n 1)其中 s_n 是归一化复频率。要做截止频率为 Ω_p 的滤波器做频率去归一化s_n s/Ω_p代入得到H(s) Ω_p² / (s² √2·Ω_p·s Ω_p²)这个形式很标准分子、分母都是二阶且分母首项系数为1。4.2 代入双线性变换的完整计算现在代入双线性变换替换式。令 d 2/T 2000记 a Ω_p 649.8s d·(1 - z⁻¹)/(1 z⁻¹)把 H(s) 的分子分母都展开。分子是 a²分母是 s² √2·a·s a²。代入后整体乘以 (1 z⁻¹)² 消去分母得到H(z) a²·(1 z⁻¹)² / [d²(1 - z⁻¹)² √2·a·d·(1 - z⁻²) a²(1 z⁻¹)²]一步步展开分母各项。这里我建议你亲自手算一遍我可以把关键数值给出来d² 4,000,000√2·a·d 1.4142 × 649.8 × 2000 ≈ 1,837,900a² 649.8² ≈ 422,240展开后按 z 的同幂次合并。z⁰ 项d² √2·a·d a² 4,000,000 1,837,900 422,240 6,260,140z⁻¹ 项-2d² 2a² -8,000,000 844,480 -7,155,520z⁻² 项d² - √2·a·d a² 4,000,000 - 1,837,900 422,240 2,584,340所以H(z) 422,240·(1 2z⁻¹ z⁻²) / (6,260,140 - 7,155,520z⁻¹ 2,584,340z⁻²)把分子分母同时除以 6,260,140得到标准的数字滤波器传递函数形式H(z) (b0 b1·z⁻¹ b2·z⁻²) / (1 a1·z⁻¹ a2·z⁻²)各项系数为系数数值b00.06745b10.13490b20.06745a1-1.1430a20.4128注意符号约定差分方程里通常写成 y[n] b0·x[n] b1·x[n-1] b2·x[n-2] - a1·y[n-1] - a2·y[n-2]这里的 a1、a2 是分母多项式的系数代入负号后就是正数参与运算。很多人的代码跑出来结果不对十有八九是这个符号约定搞反了。4.3 差分方程与C语言实现对应的差分方程y[n] 0.06745·x[n] 0.13490·x[n-1] 0.06745·x[n-2] 1.1430·y[n-1] - 0.4128·y[n-2]这个结构是最常见的二阶IIR也叫biquad。在C语言里实现非常简单typedef struct { float x[3]; // 输入历史 float y[3]; // 输出历史 } biquad_t; float biquad_process(biquad_t *f, float xn) { float yn 0.06745f * xn 0.13490f * f-x[1] 0.06745f * f-x[2] 1.1430f * f-y[1] - 0.4128f * f-y[2]; // 移位历史 f-x[2] f-x[1]; f-x[1] xn; f-y[2] f-y[1]; f-y[1] yn; return yn; }这段代码每次调用只处理一个样本点非常适合嵌入到采样中断或者音频回调里。浮点精度足够时直接算如果是MCU上想省资源可以转成Q15定点但要格外小心系数量化后的稳定性后面第5.3节细说。4.4 验证直流增益和几个关键频点推导完系数之后第一件事是验证直流增益。z 1 对应直流代入 H(z)H(1) (0.06745 0.13490 0.06745) / (1 - 1.1430 0.4128) 0.2698 / 0.2698 1直流增益是1说明这个低通滤波器直流信号可以无损通过符合设计预期。再来验证100Hz处的增益。需要一个频率响应的数值计算方法对二阶系统可以直接用复指数代入。下面这个Python函数不依赖scipy只用numpy就能算幅频响应import numpy as np def biquad_magnitude_db(b, a, fs, freq_hz): w 2 * np.pi * freq_hz / fs z np.exp(-1j * w) num b[0] b[1]*z**-1 b[2]*z**-2 den a[0] a[1]*z**-1 a[2]*z**-2 return 20 * np.log10(np.abs(num / den)) b [0.06745, 0.13490, 0.06745] a [1.0, -1.1430, 0.4128] fs 1000 for f in [0, 50, 100, 200, 400, 500]: print(f, biquad_magnitude_db(b, a, fs, f))跑到100Hz附近幅值大约是 -3dB也就是 0.707 倍左右这正是截止频率的定义400Hz以上衰减明显加快500Hz处虽然理论上是数字频率π但响应已经降到很低的电平。到这里整个设计闭环就完成了系数的正确性由直流增益和截止频点两个角度的验证背书可以放心使用。5. 双线性变换不擅长什么替代方案与工程注意点5.1 与脉冲响应不变法的取舍双线性变换不是唯一的模拟滤波器数字化手段。脉冲响应不变法impulse invariance的思路是直接对模拟冲激响应采样h[n] T·h_a(nT)然后做z变换。它最大的优点是频率映射是线性的ω Ω·T频率轴没有压缩数字滤波器的通带形状和模拟原型基本一致特别适合低通和带通这种带限场合。但致命弱点是混叠如果模拟滤波器在奈奎斯特频率以上还有明显的频响残留采样后这些高频分量会折叠回低频段破坏设计指标。高通和带阻滤波器几乎不能直接用脉冲响应不变法因为它们的频响在奈奎斯特频率处不衰减。双线性变换没有混叠问题因为频率轴被压缩到有限范围但代价就是频率非线性。实际选型逻辑很清晰要保留模拟滤波器的精确过渡带形状且信号是带限的用脉冲响应不变法要求绝对无混叠、实现简单、频率响应要求主要集中在低频段的用双线性变换。工程上后者的应用面要广得多尤其是开关电源里的环路滤波器、音频均衡器、传感器信号调理基本都是双线性变换的天下。5.2 设计高阶滤波器时的注意事项如果只是拿双线性变换设计二阶滤波器上面已经够用了。但实际工程经常碰到四阶、六阶甚至更高阶的需求。这时候直接展开成一个高阶差分方程是灾难性的高阶IIR的系数对定点量化极其敏感一个小数点后几位的舍入误差就可能导致极点移出单位圆系统从稳定变成振荡。正确做法是分解成二阶节biquad级联。以四阶巴特沃斯为例归一化原型可以分解成两个二阶节每个二阶节有自己的一组 b、a 系数。级联时把第一个二阶节的输出作为第二个二阶节的输入。这里有两个细节容易被忽略第一各节的增益分配要做均衡。模拟原型的直流增益可能不是1数字化后每一节的分子系数大小也不一样级联时如果第一节输出已经接近满幅第二节再放大就溢出了。一般是把所有节的直流增益乘起来等于总增益每一节内部保证不会超过目标电平。我用过一个简单的经验先把每节的直流增益算出来按总增益要求分配给各节再逐节验证峰值。第二二阶节的排序会影响数值性能。通常把Q值较低的节放在前面Q值较高的节放在后面这样能减少量化噪声放大。不过这个是细化优化初版可以直接按分数分解的顺序排跑通后再调整。5.3 FPGA落地时的系数定点与结构选择热搜词里有“基于FPGA的FIR数字滤波器”和“分布式算法”这里说清楚一个关键区别双线性变换设计出来的是IIR滤波器而分布式算法DA通常用来高效实现FIR。FIR没有反馈结构可以用查表加移位累加做乘累加天然适合FPGA。IIR因为有反馈回路不能直接用分布式算法的标准形式只能分解成二阶节后逐节实现每个biquad用乘累加器完成。在FPGA上实现双线性变换系数时第一件要处理的事是系数量化。比如 b1 0.13490在二进制里是无限循环小数定点化之后必然有误差。IIR滤波器的极点位置对系数精度极其敏感尤其是靠近单位圆的极点。设计流程是先用Matlab或浮点C模型确认理想系数再定点仿真最后上板实测。每一步的系数都必须统一来源否则你在仿真里验证的滤波特性和板子上跑出来的可能根本不是一个滤波器。结构选择上Direct Form II Transposed转置直接II型是IIR在FPGA上的常见选择因为它只需要一个状态变量累加链减少了寄存器和布线压力。实现细节上有三点提示一是状态变量的位宽要比输入输出宽一般至少宽4~8比特给中间累积量留出余量二是每个biquad的输出不要急着截断先全精度累加再统一量化三是如果系统时钟高于采样率很多可以用时分复用方式让一个乘法器轮流处理多个biquad资源消耗比并联实例化小一个数量级。我个人在实际项目里的体会是双线性变换最大的价值不是那个公式本身而是它把“模拟域积累几十年的滤波器设计经验”和“数字域的实现便利”桥接了起来。推导过程看懂之后再遇到任何模拟滤波器数字化的问题第一反应不再是去查表抄系数而是能自己推、自己验、自己改。按这个思路走一次完整流程你的数字滤波器设计水平会上一个台阶。