双线性变换公式推导到C语言实现：数字滤波器设计全解析

网上搜“双线性变换”，十篇有八篇上来就甩给你一个替换式：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[n+1] - 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⁻¹)/(1+z⁻¹)，代入得”，读者永远不知道这个“令”为什么成立。自己推一遍之后，以后再看到任何双线性变换的变体（例如带频率预畸变的版本）都能一眼看穿。

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 预畸变的标准操作

解决办法是在设计模拟滤波器时，先做一次频率预畸变。步骤是：

1. 确定数字滤波器目标截止频率 ω_d，单位是弧度/采样点，通常 ω_d = 2π·fc/fs。
2. 用公式 Ω_p = (2/T)·tan(ω_d/2) 反算出模拟滤波器设计要用的截止频率。
3. 用这个 Ω_p 设计模拟滤波器 H(s)。
4. 再用双线性变换数字化。

这里的逻辑是：双线性变换会把模拟频率 Ω_p 压缩成数字频率 ω_d，所以提前把设计频率抬高到 Ω_p，数字化后就正好落在 ω_d 上。

下表是采样率 1kHz、目标截止频率不同取值时，预畸变前后的模拟设计频率对比：

看到最后两行的差异了吗？频率越高，误差越大。做预畸变之后，截止频率才能精确落在目标点上。

这里必须澄清一个常见误解：预畸变只是把设计频点对准了，并不是让整条幅频曲线跟模拟原型完全重合。双线性变换的频率压缩是全局性的，通带内的形状在接近奈奎斯特频率时仍然会被压缩变形，只是截止频点这个关键节点被校正回来了。如果你的滤波器通带要求非常严格，覆盖范围又很宽，那么双线性变换可能不是最佳选择，后续可以考虑更高阶的数值方法或者直接设计数字滤波器，不走模拟原型这条路。

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.8：

s = 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,900
• a² = 649.8² ≈ 422,240

展开后按 z 的同幂次合并。z⁰ 项：

d² + √2·a·d + a² = 4,000,000 + 1,837,900 + 422,240 = 6,260,140

z⁻¹ 项：

-2d² + 2a² = -8,000,000 + 844,480 = -7,155,520

z⁻² 项：

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⁻²)

各项系数为：

注意符号约定：差分方程里通常写成 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，资源消耗比并联实例化小一个数量级。

我个人在实际项目里的体会是，双线性变换最大的价值不是那个公式本身，而是它把“模拟域积累几十年的滤波器设计经验”和“数字域的实现便利”桥接了起来。推导过程看懂之后，再遇到任何模拟滤波器数字化的问题，第一反应不再是去查表抄系数，而是能自己推、自己验、自己改。按这个思路走一次完整流程，你的数字滤波器设计水平会上一个台阶。