设计并应用一个IIR-ButterWorth-Filter的例子

设计一个butter-worth filter

为何是巴特沃斯滤波器

首先我们需要知道我们为什么要设计一个巴特沃斯滤波器，巴特沃斯滤波器有一个重要的特点就是最大平坦通带​​（Maximally Flat Passband）：在通带内（0 ≤ ω ≤ ωₐ）幅度响应尽可能平坦（无波纹）

如图所示

这个函数是幅度平方函数，即

考虑一阶滤波器

step 1 定义幅度平方响应

$$|H(j\omega ){|}^2=\frac{1}{1+(\frac{\omega }{{\omega }_c}{)}^{2}n}$$

其中n是滤波器阶数，$w_c$指的是截止频率

step 2 解析延拓到s域

用 $$s=j\omega$$替换，得到s域表达式:

$H(s)H(-s)=\frac{1}{(\frac{s}{j}{)}^2}=\frac{1}{1-s^2}$

step 3 获得极点

求分母的根，得到 $s=\pm 1$

我们选择左边平面极点

$p=-1$

step 4 构造传递函数

其实很容易观察得到，$H(s)=\frac{1}{s+1}$

二阶巴特沃斯滤波器推导

Step 1 定义幅度平方响应

幅度平方函数的一般形式：
$$|H(j\omega ){|}^2=\frac{1}{1+{\left(\frac{\omega }{{\omega }_c}\right)}^{2n}}$$

对于二阶滤波器（n=2）：
$$|H(j\omega ){|}^2=\frac{1}{1+{\left(\frac{\omega }{{\omega }_c}\right)}^4}$$

---

Step 2 解析延拓到s域

用 $s=j\omega$ 替换，得到s域表达式：
$$H(s)H(-s)=\frac{1}{1+{\left(\frac{s}{j{\omega }_c}\right)}^4}$$

将 $j^4=1$ 代入：
$$H(s)H(-s)=\frac{1}{1+{\left(\frac{s}{{\omega }_c}\right)}^4}$$

---

Step 3 获得极点

求分母的根（极点）：
$$1+{\left(\frac{s}{{\omega }_c}\right)}^4=0$$

解得：
$$s_k={\omega }_c\cdot e^{j\frac{\pi }{4}(2k+1)}(k=0,1,2,3)$$

具体极点位置（取 ${\omega }_c=1$ 归一化）：
$$\begin{array}{rl}{s}_0& =e^{j\frac{\pi }{4}}=\frac{\sqrt{2}}{2}+j\frac{\sqrt{2}}{2}\\ s_1& =e^{j\frac{3\pi }{4}}=-\frac{\sqrt{2}}{2}+j\frac{\sqrt{2}}{2}\\ s_2& =e^{j\frac{5\pi }{4}}=-\frac{\sqrt{2}}{2}-j\frac{\sqrt{2}}{2}\\ s_3& =e^{j\frac{7\pi }{4}}=\frac{\sqrt{2}}{2}-j\frac{\sqrt{2}}{2}\end{array}$$

选择左半平面极点（$s_1$ 和 $s_2$）以保证系统稳定。

---

Step 4 构造传递函数

将左半平面极点组合成传递函数：
$$H(s)=\frac{{\omega }_c^2}{(s-s_1)(s-s_2)}$$

代入具体值（${\omega }_c=1$）：
$$\begin{array}{rl}H(s)& =\frac{1}{\left(s+\frac{\sqrt{2}}{2}-j\frac{\sqrt{2}}{2}\right)\left(s+\frac{\sqrt{2}}{2}+j\frac{\sqrt{2}}{2}\right)}\\ & =\frac{1}{s^2+\sqrt{2}s+1}\end{array}$$

最终标准形式：
$$H(s)=\frac{{\omega }_c^2}{s^2+\sqrt{2}{\omega }_{c}s+{\omega }_c^2}$$

---

关键特性

1. -40 dB/decade 高频衰减
2. 最大平坦通带（Butterworth特性）
3. 截止频率处增益为 $-3\text{dB}$（即 $|H(j{\omega }_c)|=\frac{1}{\sqrt{2}}$）

对比一阶二阶巴特沃斯滤波器

详细推导过程

1.连续域原型（模拟滤波器）

二阶归一化巴特沃斯低通滤波器（截止频率ωₐ=1 rad/s）：

$H_{analog}(s)=\frac{1}{s^2+\sqrt{2}s+1}$ 显然极点位置 $s=-\frac{-\sqrt{2}}{2}\pm j\frac{\sqrt{2}}{2}$

-3db的点:$\omega =1rad/s$

2.频率预畸变（关键代码变量 C）

双线性变换会导致频率扭曲，需对截止频率预补偿：

${\omega }_c=2\pi f_c$

预畸变公式：

$$C=\frac{1}{tan(\frac{{\omega }_{c}T}{2})}=\frac{1}{tan(\frac{\pi f_c}{f_s})}$$

3. 双线性变换

双线性变换公式

$s=\frac{2}{T}\cdot \frac{1-z^{-1}}{1+z^{-1}}$ 其中T为采样间隔

双线性变换会将模拟频率(${\omega }_a$)和数字频率(${\omega }_d$)的非线性映射

${\omega }_a=\frac{2}{T}tan(\frac{{\omega }_{d}T}{2})$

这种映射会导致低频段近似相等，高频段扭曲严重

为了保证数字滤波器的截止频率${\omega }_d={\omega }_c$准确对应模拟滤波器中的${\omega }_a={\omega }_c$，需要对模拟截止频率进行预畸变：

${\omega }_a={\omega }_c=\frac{2}{T}tan\frac{{\omega }_{c}T}{2}$

从而引入预畸变系数
$$C=\frac{1}{tan(\frac{{\omega }_{c}T}{2})}=\frac{1}{tan(\frac{\pi f_c}{f_s})}$$

//C++ 中如此表示
double C = 1 / std::tan(M_PI * center_freq / (double)PublicVar::ins().sample_rate);

修正后的双线性变换 $s=C\cdot \frac{1-z^{-1}}{1+z^{-1}}$

可以尝试一下

4：代入双线性变换

将 $s=C\cdot \frac{1-z^{-1}}{1+z^{-1}}$ 代入连续传递函数：
$$H(z)=\frac{1}{{\left(C\cdot \frac{1-z^{-1}}{1+z^{-1}}\right)}^2+\sqrt{2}\left(C\cdot \frac{1-z^{-1}}{1+z^{-1}}\right)+1}$$

5：通分整理分母

分母部分展开：
$$\begin{array}{rl}\text{分母}& =C^2{\left(\frac{1-z^{-1}}{1+z^{-1}}\right)}^2+\sqrt{2}C\left(\frac{1-z^{-1}}{1+z^{-1}}\right)+1\\ & =\frac{C^2(1-z^{-1}{)}^2+\sqrt{2}C(1-z^{-1})(1+z^{-1})+(1+z^{-1}{)}^2}{(1+z^{-1}{)}^2}\end{array}$$

6：展开分子

分子是 $(1+z^{-1}{)}^2$，分母展开为：
$$\begin{array}{rl}\text{分子}& =C^2(1-2z^{-1}+z^{-2})+\sqrt{2}C(1-z^{-2})+(1+2z^{-1}+z^{-2})\\ & =(C^2+\sqrt{2}C+1)+(-2C^2+2)z^{-1}+(C^2-\sqrt{2}C+1)z^{-2}\end{array}$$

7：合并传递函数

组合分子分母：
$$H(z)=\frac{(1+z^{-1}{)}^2}{(C^2+\sqrt{2}C+1)+(-2C^2+2)z^{-1}+(C^2-\sqrt{2}C+1)z^{-2}}$$

8：归一化分母

将分母首项系数归一化为1：
$$H(z)=\frac{\frac{1}{C^2+\sqrt{2}C+1}(1+2z^{-1}+z^{-2})}{1+\frac{-2C^2+2}{C^2+\sqrt{2}C+1}{z}^{-1}+\frac{C^2-\sqrt{2}C+1}{C^2+\sqrt{2}C+1}{z}^{-2}}$$

9：系数对应

定义系数：
$$\begin{array}{rl}{b}_0& =\frac{1}{C^2+\sqrt{2}C+1}\\ b_1& =2b_0\\ b_2& =b_0\\ a_1& =\frac{2-2C^2}{C^2+\sqrt{2}C+1}\\ a_2& =\frac{C^2-\sqrt{2}C+1}{C^2+\sqrt{2}C+1}\end{array}$$

最终差分方程：
$$y[n]=b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]-a_{1}y[n-1]-a_{2}y[n-2]$$

python代码示例

import mathclass TFZ_coefficients:def __init__(self):self.b0 = 0.0self.b1 = 0.0self.b2 = 0.0self.a1 = 0.0self.a2 = 0.0class PublicVar:_instance = None@classmethoddef ins(cls):if cls._instance is None:cls._instance = cls()return cls._instancedef __init__(self):self.sample_rate = 44100.0  # Default value, adjust as neededdef MakeCoeffLowPass(center_freq):ret = TFZ_coefficients()C = 1 / math.tan(math.pi * center_freq / float(PublicVar.ins().sample_rate))SqC = math.pow(C, 2)MultC = math.sqrt(2) * Cc = 1.0 / (1.0 + MultC + SqC)ret.b0 = cret.b1 = 2.0 * cret.b2 = cret.a1 = 2.0 * c * (1.0 - SqC)ret.a2 = c * (1.0 - MultC + SqC)return ret