Butterworth型IIR滤波器的设计与C语言实现

仲夏夜的乙醇

双非CS专业菜鸡一枚

前言

最近在公司实习，帮忙解决了一个之前4个人都解决不了的技术难点——MEP（运动诱发电位）信号的滤波，我是上来接手的第五人。我设计了 巴特沃夫滤波器 的参数，以及 滤波器的级联 方式。最后是我拿着信号发生器调出不同频率的正弦波做测试，在公司研发中心的技术总监面前演示了滤波的效果。现在想把自己做项目学到的东西记录下来，万一以后还会用到呢。

1、求取n阶低通模拟滤波器的极点

在计算前，先确定IIR滤波器的阶数 n ，以及截止频率 f_{c} 。如果要设计的Butterworth型IIR滤波器是带通或者是带阻的，就要先确定截止频率的下限 f_{L} 和上限 f_{H} 。

一般地，对于n阶Butterworth型低通模拟滤波器，其n个极点 p_{i} ,i\in[0,n-1] 均匀地分布在s平面单位圆的左半侧。即 p_{i}= e^{j\pi\frac{n+1+2i}{2n}} 。

该低通模拟滤波器无零点。

2、双线性变换

利用 双线性变换 可以将低通模拟滤波器转化为低通数字滤波器，进而获得滤波器的传递函数（本文的步骤6）。这里先描述 低通滤波器 的处理方法，高通、带通、带阻见下文。

2.1 预畸变（Prewarp）

记采样周期 T_s 为采样频率 f_s 的倒数，即 T_s=\frac{1}{f_s} ，预畸变的原始公式为： \Omega = \frac{2}{T_s}tan\frac{\omega}{2} 。其中 \Omega 为模拟角频率， \omega 为数字角频率。预畸变的公式可 由双线性变换的公式推导 得来。

实际编程计算时，在原始公式的基础上把 T_s 两边同乘，让等式左边包含 T_s ，方便下一步做双线性变换的时候整体带进去计算。同时令 \frac{\omega}{2} =\frac{f_c}{f_s} ，即可得预畸变的系数 w=\Omega T_s=2tan\frac{f_c}{f_s} 。

为什么可以这样代换呢？在数字信号处理的理论中， \omega 是归一化后的，即都除以基准频率（信号的奈奎斯特频率）。于是有 \frac{\omega}{2} = \frac{f_c / \frac{f_s}{2}}{2} =\frac{f_c}{f_s} 。

随后对每个极点都乘以 w ，改变极点在s平面的位置。

2.2 双线性变换（Tustin变换）

将模拟滤波器极点所处的s平面映射至z平面，按下式进行： z = \frac{2+sT_s}{2-sT_s} 。分子、分母各是一个线性变换，因此得名“双线性变换”。

实际编程计算时，预畸变后的零、极点相当于上式中的 sT_s 。记 x = sT_s ，直接按照 z = \frac{2+x}{2-x} 进行变换。变换后随后按照步骤6，从滤波器的零、极点导出滤波器的传递函数。对于 低通模拟滤波器 而言并 没有零点 ，在实际转化为数字滤波器的时候 需要在z平面额外加入n个坐标为（-1,0）的零点 ，保证传递函数分子分母多项式的次数相同。

双线性变换的公式有一个很简单的推导过程：从 z平面和s平面之间的映射 z=e^{sT_s} 出发，改写为： z=\frac{e^{sT_s/2}}{e^{-sT_s/2}} ，利用 e^x=1+x+x^2+... 这一泰勒展开，展开至一阶项，代入分子分母得 z=\frac{1+sT_s/2}{1-sT_s/2} 。该式即为双线性变换的原始公式。

3、低通转高通滤波器的处理

与低通滤波器类似，也需要进行预畸变。但原来低通模拟滤波器的极点都要取倒数，即新的极点 p'_i=\frac{1}{p_i} 。此外，新的极点都要做一次双线性变换。下文的带通、带阻也类似。

与低通滤波器不同的是，双线性变换前在s平面额外引入n个位于坐标原点的零点，与极点一同进行双线性变换。

4、低通转带通滤波器的处理

这个过程稍微有点复杂。首先利用截止频率下、上限分别算出两个预畸变的系数 w_L=2tan\frac{f_L}{f_s} ， w_H=2tan\frac{f_H}{f_s} 。随后产生新的2n个极点，新的极点 p'_k,k\in[0,2n-1] 由以下两个公式得出：

p'_i=\frac{w_H-w_L}{2}p_i(1+\sqrt{1-\frac{w_Lw_H}{{(\frac{w_H-w_L}{2}p_i)}^{2}}}) ， p'_{n+i}=\frac{w_H-w_L}{2}p_i(1-\sqrt{1-\frac{w_Lw_H}{{(\frac{w_H-w_L}{2}p_i)}^{2}}})

至于零点的分布，在原有低通模拟滤波器的基础上引入n个位于坐标原点的零点。在双线性变换后，再添加n个坐标为（-1,0）的零点。

5、低通转带阻滤波器的处理

与带通类似，首先分别算出两个预畸变的系数 w_L=2tan\frac{f_L}{f_s} ， w_H=2tan\frac{f_H}{f_s} 。随后产生新的2n个极点，新的极点 p'_k,k\in[0,2n-1] 由以下两个公式得出：

p'_i=\frac{w_H-w_L}{2p_i}(1+\sqrt{1-\frac{w_Lw_H}{{(\frac{w_H-w_L}{2p_i})}^{2}}}) ， p'_{n+i}=\frac{w_H-w_L}{2p_i(}1-\sqrt{1-\frac{w_Lw_H}{{(\frac{w_H-w_L}{2p_i})}^{2}}})

至于零点的分布，在原有低通模拟滤波器的基础上引入2n个零点，前n个位于（0， \sqrt{w_Lw_H} ），后n个位于（0， -\sqrt{w_Lw_H} ）。

6、计算滤波器的传递函数

假设数字滤波器有 M 个零点 z_i ， N 个极点 p_i 。一个传递函数可以写成：

H(z)=\frac{\prod_{i=0}^{M-1}{(x-z_i)}}{\prod_{i=0}^{N-1}{(x-p_i)}}

对这个式子展开成多项式的形式，分别获得传递函数分子多项式的系数 b_i 和分母多项式的系数 a_i 。在本文的滤波器中， M=N ：

H(z)=\frac{\sum_{i=0}^{N}{b_ix^{N-i}}}{\sum_{i=0}^{N}{a_ix^{N-i}}}

做完这一步，一个数字滤波器就设计得差不多了，剩下来就要将传递函数的系数归一化啦~

7、计算增益

滤波器在频率 f 处的增益为 \left| H(e^{j2\pi\frac{f}{f_s}}) \right| 。根据滤波器的不同，分四种情况讨论：

7.1 低通滤波器的增益

将 z=e^{j2\pi\cdot0}=1 代入上一步得出的传递函数，所得结果再取绝对值，即为低通滤波器的直流增益。

7.2 高通滤波器的增益

高通滤波器的增益指的是在奈奎斯特频率处的增益。由奈奎斯特频率的定义可知，奈奎斯特频率的值为采样频率的一半，即 \frac{f_s}{2} 。将 z=e^{j2\pi\frac{f_{s}/2}{f_{s}}}=e^{j\pi}=-1 代入传递函数计算，所得结果再取绝对值，即为高通滤波器的增益。

7.3 带通滤波器的增益

带通滤波器的增益指的是在中心频率 f_{C}=\frac{f_L+f_H}{2} 处的增益。将 z=e^{j2\pi\frac{(f_L+f_H)/2}{f_{s}}}=e^{j\pi\frac{f_L+f_H}{f_{s}}} 代入传递函数计算，所得结果再取绝对值，即为带通滤波器的增益。

7.4 带阻滤波器的增益

分别计算带阻滤波器的直流增益和位于奈奎斯特频率处的增益，取二者乘积的平方根（这里不是取绝对值了，因为负负得正），即为带阻滤波器的增益。

把上面四种不同的情况的增益都用 G 表示，以便计算IIR的系数。

8、计算IIR系数

我们知道，一个n阶IIR的输出为： y_i=\sum_{j=0}^{N}{b_jx_{i-j}} - \sum_{j=1}^{N}{a_jy_{i-j}} 。该式的 b_j 和 a_j 系数（用j做下标是为了和上面传递函数的系数区分）可由传递函数分子分母多项式系数导出。

（1） a_j 系数：传递函数分母多项式系数 a_i 除以分母首项系数 a_0 ，即为 a_j 系数

（2） b_j 系数：传递函数分子多项式系数 b_i 除以分母首项系数 a_0 ，再除以增益 G ，即为 b_j 系数

9、直接I型IIR滤波器的C语言实现

大多数IIR滤波器的实现都是直接I型的， 直接套用IIR输出的定义式 即可。具体的系数 可使用MATLAB的工具箱进行设计 ，也可以按照上文的方法计算求解。值得注意的是，对于过往的 x_i 和 y_i 可以使用一个环形数组进行存储，降低时间复杂度。本文给出的代码以三阶直接I型IIR滤波器为例，导入MATLAB生成的系数可改造为任意阶数的直接I型IIR滤波器。

#define N 3  //IIR滤波器的阶数
double a[N + 1] = {1,1,2,3}; //对应前面的a[j]系数
double b[N + 1] = {1,2,3,4}; //对应前面的b[j]系数
double old_x[N-1] = {0};
double old_y[N-1] = {0};
int pos = 0;
double iir_step(double input)
    int i,p;
    double output = b[0] * input;
    for (i = 1; i <= N; i++) {
        p = (pos + N - i) % N;
        output += b[i] * old_x[p] - a[i] * old_y[p];
    if (N > 1) {
        old_x[pos] = input;
        old_y[pos] = output;
        pos = (pos + 1) % N;
    return output;
}

10、直接II型IIR滤波器的C语言实现

n阶的直接II型IIR滤波器按照n+1个差分方程进行计算：

y_i=b_0x_i+d_1

d_1=b_1x_i+a_1y_i+d_2

（中间省略n-3个差分方程）

d_{n-1}=b_{n-1}x_i+a_{n-1}x_i+d_n

d_n=b_nx_i+a_ny_i

直接II型IIR滤波器框图可以查阅这篇文章：飞控中的IIR二阶滤波器 - 知乎 (zhihu.com)

#define N 3  //IIR滤波器的阶数
double a[N + 1] = { 1,1,2,3 }; //对应前面的a[j]系数
double b[N + 1] = { 1,2,3,4 }; //对应前面的b[j]系数
double d[N] = {0};
double iir_step(double input)