信号处理基础 | 巴特沃斯IIR低通滤波器的设计：从频率响应到稳定模拟原型，再到离散实现

本文是二阶 Butterworth 低通滤波器系列文章的第 1 篇，主要介绍 Butterworth 数字 IIR 低通滤波器，是怎么从“频率指标”一步步设计出来的？ 本文示例的二阶 Butterworth 低通滤波器的设计参数如下：

滤波器类型：Low-pass
设计方法：Butterworth
阶数：2
归一化 3-dB 截止频率：0.1
频率归一化方式：1.0 = Nyquist frequency

生成的频响、阶跃响应和零极点图如下：

本文和仿真代码同步发在个人博客网站 VoxWorking · 声学与音频知识平台，网站定期更新音频处理、语音处理和声学相关知识，欢迎大家浏览交流~

1. 数字滤波器设计的整体路线

IIR 数字滤波器的一种经典设计路线是：

频率响应指标
    ↓
模拟低通原型 H_a(s)
    ↓
频率变换或截止频率缩放
    ↓
双线性变换
    ↓
数字滤波器 H(z)
    ↓
差分方程或二阶节实现

本例是低通滤波器，且已经确定为 Butterworth，因此路线可以更具体地写成：

Butterworth 幅度平方响应
    ↓
取 N = 2
    ↓
求出模拟原型极点
    ↓
只选择 s 平面左半平面极点，得到稳定 H_a(s)
    ↓
根据数字截止频率做预畸变
    ↓
用双线性变换把 H_a(s) 变成 H(z)

这里需要注意：（1）幅度响应和模拟原型 H(s) 是两个不同的东西，我们是要根据期望的幅度响应推导出相应的模拟原型；（2）模拟域s和数字域z的稳定性条件不一样：模拟滤波器看 s 平面左半平面稳定，数字滤波器看 z 平面单位圆内稳定。

2. 期望的幅度响应

Butterworth 滤波器的设计目标是通带尽可能平坦，其模拟低通原型幅度平方响应为：

|H_a(jΩ)|^2 = 1 / (1 + (Ω / Ωc)^(2N))

其中：

H_a(s) ：模拟滤波器传递函数
Ω      ：模拟角频率
Ωc     ：模拟截止角频率
N      ：滤波器阶数

由上述的幅度平方响应，我们因为可以看出来，当 Ω < Ωc 时, |H_a(jΩ)|^2越接近于1；当 Ω > Ωc 时且Ω越来越大时, |H_a(jΩ)|^2会越来越小；这正是我们期望的低通滤波的效果。当 Ω = Ωc 时：

|H_a(jΩc)|^2 = 1 / 2
|H_a(jΩc)| = 1 / sqrt(2)

换算成 dB：

20log10(1 / sqrt(2)) ≈ -3.0103 dB

所以 Butterworth 的截止频率通常就是 -3 dB 频率。本例是二阶滤波器，让N = 2得到：

|H_a(jΩ)|^2 = 1 / (1 + (Ω / Ωc)^4)

这个式子就是二阶低通幅度响应的表示，其说明了不同频率应该被放大或衰减多少，但它不是可以直接实现的系统函数。

3. 系统函数 H(s)

幅度平方响应只描述幅度，不描述相位。可是具体系统函数需要不仅决定幅度，也决定相位、极点位置和稳定性。幅度平方响应可写成

|H_a(jΩ)|^2 = H_a(jΩ) H_a(jΩ)^*

如果系统系数是实数，则：

H_a(jΩ)^* = H_a(-jΩ)

所以：

|H_a(jΩ)|^2 = H_a(jΩ) H_a(-jΩ)

把s = jΩ看成复频域变量，就可以写成：

H_a(s)H_a(-s)

4. 从幅度平方响应得到候选极点

我们从二阶 Butterworth 幅度平方响应开始：

|H_a(jΩ)|^2 = 1 / (1 + (Ω / Ωc)^4)

因为：

s = jΩ
Ω = s / j

所以：

Ω / Ωc = s / (jΩc)

代入后得到：

H_a(s)H_a(-s) = 1 / (1 + (s / jΩc)^4)

由于：

j^4 = 1

所以：

(s / jΩc)^4 = s^4 / Ωc^4

因此分母可以写成：

1 + s^4 / Ωc^4

要求候选极点，就令分母为零：

1 + s^4 / Ωc^4 = 0

也就是：

s^4 = -Ωc^4

把 -1 写成复指数形式：

-1 = e^(jπ), e^(j3π), e^(j5π), e^(j7π), ...

于是四个根为：

s = Ωc e^(jπ/4)
s = Ωc e^(j3π/4)
s = Ωc e^(j5π/4)
s = Ωc e^(j7π/4)

这四个极点都在半径为 Ωc 的圆上。它们分布在 s 平面中，其中两个在左半平面，两个在右半平面。

5. 基于稳定性只选择左半平面的极点

模拟连续时间系统的稳定条件是，所有极点必须位于 s 平面左半平面，即Re(s) < 0，四个候选极点中：

Ωc e^(jπ/4)   ：右半平面
Ωc e^(j3π/4)  ：左半平面
Ωc e^(j5π/4)  ：左半平面
Ωc e^(j7π/4)  ：右半平面

因此稳定二阶 Butterworth 模拟原型选择：

p1 = Ωc e^(j3π/4)
p2 = Ωc e^(j5π/4)

利用欧拉公式：

e^(jθ) = cosθ + j sinθ

可以得到：

p1 = -Ωc/sqrt(2) + jΩc/sqrt(2)
p2 = -Ωc/sqrt(2) - jΩc/sqrt(2)

这两个极点是一对共轭复数，并且实部都小于 0，所以对应稳定模拟系统。

6. 由稳定极点得到模拟原型 H_a(s)

如果传递函数有两个极点p1, p2，那么它的分母可以写成：

(s - p1)(s - p2)

代入：

p1 = -Ωc/sqrt(2) + jΩc/sqrt(2)
p2 = -Ωc/sqrt(2) - jΩc/sqrt(2)

得到：

s - p1 = s + Ωc/sqrt(2) - jΩc/sqrt(2)
s - p2 = s + Ωc/sqrt(2) + jΩc/sqrt(2)

所以：

(s - p1)(s - p2)
= (s + Ωc/sqrt(2) - jΩc/sqrt(2))
  (s + Ωc/sqrt(2) + jΩc/sqrt(2))

利用：

(a - jb)(a + jb) = a^2 + b^2

令：

a = s + Ωc/sqrt(2)
b = Ωc/sqrt(2)

得到：

(s - p1)(s - p2)
= (s + Ωc/sqrt(2))^2 + (Ωc/sqrt(2))^2

展开：

= s^2 + sqrt(2)Ωc s + Ωc^2/2 + Ωc^2/2
= s^2 + sqrt(2)Ωc s + Ωc^2

因此二阶 Butterworth 模拟低通原型为：

H_a(s) = Ωc^2 / (s^2 + sqrt(2)Ωc s + Ωc^2)

分子取 Ωc^2，是为了让直流增益为 1：

H_a(0) = Ωc^2 / Ωc^2 = 1

至此，我们完成了从“幅度响应目标”到“稳定模拟低通原型”的推导。

7. 双线性变换：模拟域转离散域

现在我们有了稳定模拟低通原型（模拟域系统函数）：

H_a(s) = Ωc^2 / (s^2 + sqrt(2)Ωc s + Ωc^2)

现在我们需要把这个系统函数从模拟域s转到离散域z，这其中会用到双线性变换，如下：

s = 2/T · (1 - z^-1) / (1 + z^-1)

由于篇幅限制和更容易了解滤波器设计过程，我们暂时先不推导为什么这个双线性变换是这个样子，暂时记住它是基于“ 用梯形积分法把连续时间系统离散化”的思路，进行模拟域和离散域之间转换的。

由于双线性变换的存在，模拟域和离散域之间的频率轴存在非线性的关系，令数字频率z = e^(jω),并代入双线性变换公式，得到：

s = 2/T · (1 - e^(-jω)) / (1 + e^(-jω))

化简得到

s = j · 2/T · tan(ω/2)

由于模拟频率轴是s = jΩ，所以有

Ω = 2/T · tan(ω/2)

所以

Ωc = 2/T · tan(ωc/2)

如果令采样周期归一化T = 1，则Ωc = 2 tan(ωc / 2)，其中ωc = πWn 所以Ωc = 2 tan(πWn / 2)

为了让公式简洁，令K = tan(πWn / 2), 则Ωc = 2K, 然后把

s = 2(1 - z^-1)/(1 + z^-1)
Ωc = 2K

代入模拟模型函数，得到

H(z) = K²(1 + z^-1)² / 
       [(1 - z^-1)² + sqrt(2)K(1 - z^-1)(1 + z^-1) + K²(1 + z^-1)²]

展开分母并合并同类项可得到分母：

[1 + sqrt(2)K + K²] + [2(K² - 1)] z^-1 + [1 - sqrt(2)K + K²] z^-2

为了让数字滤波器分母首项 a0 = 1，需要把分子分母都除以常数项：

D = 1 + sqrt(2)K + K²

所以经过以上推导，离散差分方程：

          b0 + b1 z^-1 + b2 z^-2
H(z) = --------------------------------
          1 + a1 z^-1 + a2 z^-2

中的各个系数值是：

b0 = K² / D
b1 = 2K² / D
b2 = K² / D

a0 = 1
a1 = 2(K² - 1) / D
a2 = (1 - sqrt(2)K + K²) / D

本例 Wn = 0.1，得到：

b = [0.0200833655642112, 0.0401667311284225, 0.0200833655642112]
a = [1.0, -1.5610180758007182, 0.6413515380575632]

对应差分方程为：

y[n] = b0 x[n] + b1 x[n-1] + b2 x[n-2]
       - a1 y[n-1] - a2 y[n-2]

这就是最终可以在代码中运行的数字 IIR 低通滤波器。

8. 小结

本文讲解了Butterworth IIR 数字低通滤波器的设计过程：

下篇文章我们会用幅频、相频、群延迟、阶跃响应、零极点图讲解设计结果。 一起加油吧 ~ ~ ~