数字信号处理6：IIR滤波器设计

IIR滤波器设计

1. 简介

IIR滤波器（无限冲激响应滤波器）相比FIR（有限冲激响应滤波器），是一种在运算效率方面非常高效的滤波器，相比于FIR滤波器，在实现同等的滤波效果时，IIR滤波器所需的阶数要比FIR滤波器低很多。FIR滤波器的阶数一方面取决于通带与阻带之间的增益（低阶滤波器应用高增益，会导致吉布斯效应），一方面取决于滤波精度（FIR滤波器的精度为$F_s/N$，$F_s$为采样率，$N$为滤波器阶数）。而IIR滤波器阶数一方面也是由通带和阻带之间的增益差决影响的，一方面则是由通带与阻带之间的频率的比值决定。但是IIR滤波器有两个缺点：相位响应不是线性的；阶数过高时难以稳定。

通常，定义一个线性时不变系统的传输函数如下：

$$\begin{gathered} y[n]+a[1]y[n-1]+a[2]y[n-2]+\cdots +a[M]y[n-M] \\ =b[0]x[n]+b[1]x[n-1]+\cdots +b[N]x[n-N]\text{\hspace{1em}(1)} \end{gathered}$$

对于FIR滤波器来说，输出只与输入有关，即$a[n]=0$，而IIR滤波器的输出不仅与输入有关，还与过去的输出有关，即$a[n]$不总是为0。

2. 设计步骤简明

IIR滤波器设计相比于FIR滤波器，所需的计算和步骤都更为复杂，很多书上开始讲这一点都从艰深的步骤开始，而对于我这样的数学不好的新手来说则是一头雾水。实际上在一开始掌握IIR滤波器设计的大体流程对于实用学习来说是非常有益的。接下来我就从工程实用的角度来大概说明IIR滤波器设计的几个重要方法和步骤。

1. 与FIR滤波器可以直接从频域设计频率响应然后通过傅里叶变换得到时域滤波器不同，IIR滤波器无法从时域直接设计响应，而是通过经典模拟滤波器转换而来。几种经典模拟滤波器分别是巴特沃兹滤波器、切比雪夫I型和II型滤波器、椭圆滤波器。而这几种滤波器都是有参考表的，可以根据滤波器阶数直接确定模拟滤波器系数。
2. 以上所得的模拟滤波器都是低通形式的，因此要采用一些方法将其转换为带通、带阻和高通滤波器。
3. 将上述转换后的不同形式的模拟滤波器转换成数字滤波器。转换的方法有冲激不变法、双线性变换法。最常用的是双线性变换法，这也是这篇博客所讲的方法。
4. 若要像FIR滤波器那样，针对不同频段有不同的增益，那么就需要将IIR滤波器进行串联或并联，后面我会讲一篇音频EQ的例子，它就是使用了IIR滤波器并联。

以上基本就是IIR滤波器的总结，其中出现了几个大家会在教材里经常见到的关键字。不过IIR滤波器和模拟滤波器相关，因此接下来首先要对拉普拉斯变换和Z变换有一个大概的认识。

3. 拉普拉斯变换和Z变换

这里只是简单让大家对这两种变换有个印象，不会具体讲理论。需要的同学可以查看其他教材。

3.1 拉普拉斯变换

定义拉普拉斯变换如下：

$$X(s)={\int }_{-\infty }^{+\infty }x(t)e^{-st}\text{\hspace{1em}(2)}$$

其中，$s=\sigma +j\Omega$。

拉普拉斯变换是描述模拟域的。如果将$\sigma$和$w$分别作为平面坐标系的x轴和y轴，这个平面就是我们定义的s平面。

可以看到，拉普拉斯变换实际上就是为将信号幅值增益为原来的$e^{\sigma }$倍。如果$\sigma =0$，则特化为连续信号的傅里叶变换，落在s平面上就是虚轴。

一个信号的拉普拉斯变换是有收敛域的。如果x(t)的拉普拉斯变换是有理的，如果x(t)是右边信号，那么其收敛域ROC就位于s平面上最右边极点的右边。如果x(t)是左边信号，那么其收敛域就位于s平面上最左边极点的左边。 我们实际处理的信号都是有理的且为右边信号，因此适用于第一种情况。
关于更详细的拉普拉斯变换的知识，可以参考奥本海默的《信号与系统》。

3.2 Z变换

定义Z变换为

$$X[z]=\sum _{n=-\infty }^{+\infty }x[n]z^{-n}\text{\hspace{1em}(3)}$$

其中，$z=r{e}^{jw}$。
若将$r$看作$z$的幅值，$w$看作$z$的相角，则两者的关系可看作是一个极坐标表示的平面，称为z平面。

和拉普拉斯变换一致，可以看作是将信号$x[n]$的幅值增大到$r$倍。如果$r=0$，那么就特化离散信号的傅里叶变换，落在z平面的单位圆上。

基于Z变换，我们可以将数字域LTI系统的转换方程用两个多项式的比来表示。
首先，基于以下推导：

$$\begin{gathered} Z(x[n-k])=\sum _{n=-\infty }^{+\infty }x[n-k]z^{-n}=\sum _{n=-\infty }^{+\infty }x[n-k]z^{-(n-k)}{z}^{-k} \\ =z^{-k}\sum _{n=-\infty }^{+\infty }x[n-k]z^{-(n-k)}=z^{-k}X[z] \end{gathered}$$

应用到式1

$$\begin{gathered} y[n]+a[1]y[n-1]+a[2]y[n-2]+\cdots +a[M]y[n-M] \\ =b[0]x[n]+b[1]x[n-1]+\cdots +b[N]x[n-N] \\ \Downarrow \\ Y[z]+a[1]Y[Z]z^{-1}+a[2]Y[z]z^{-2}+\cdots +a[M]Y[z]z^{-M} \\ =b[0]X[z]+b[1]X[Z]z^{-1}+\cdots +b[N]X[z]z^{-N} \end{gathered}$$

令

$$H[z]=\frac{Y[z]}{X[z]}=\frac{{\sum }_{k=0}^{N}b[k]z^{-k}}{1+{\sum }_{k=1}^{M}a[k]z^{-k}}\text{\hspace{1em}(4)}$$

为LTI系统的转换方程。
也可以写成

$$H[z]=\frac{Y[z]}{X[z]}=\frac{{\sum }_{k=0}^{N}b[k]z^{-k}}{{\sum }_{k=0}^{M}a[k]z^{-k}}$$

此时$a[0]=1$。

可以将其变换为：

$$H[z]=g\frac{(z-z_0)(z-z_1)\dots (z-z_N)}{(z-p_0)(z-p_1)\dots (z-z_M)}\text{\hspace{1em}(5)}$$

其中，$z_k$为零点，$p_k$为极点。同拉普拉斯变换一样，Z变换也是有收敛域的。如果x[n]是右边序列，且x[n]的Z变换X[z]是有理的，那么其收敛域就位于模最大的极点的外边。

4. 双线性变换法

4.1 模拟域与数字域的映射

不论是冲激响应不变法还是双线性变换法，都是为了将模拟域滤波器转换到数字域滤波器。冲激响应不变法使用了一种较为直观的方式，直接对模拟滤波器的冲激响应进行采样，按照上一节讲的，直观地使$z=e^{sT}$，$T$为采样周期。但是这样一来就会引入频谱混叠，因为这样的转换关系使得s平面的$j\Omega$轴每隔$\frac{2\pi }{T}$便会映射到z平面上的单位圆一周。 因此冲激响应不变法不宜构造高通滤波器。

因此，要使用另一种变换关系，这种关系应保证：

1. s平面的整个$j\Omega$轴映射到z平面上的单位圆一周。
2. 若$G(s)$是稳定的，那么$H[z]$也是稳定的。
3. 这种映射是可逆的，既能由$G(s)$得到$H[z]$，也能由$H[z]$得到$G(s)$。
4. 如果$G(j0)=1$，那么$H[e^{j0}]$也应等于1。

双线性变换便是满足上面四个条件的一种变换形式：

$$s=\frac{2}{T}\frac{z-1}{z+1}\text{\hspace{1em}(6)}$$

其反向转换为：

$$z=\frac{1+\frac{T}{2}s}{1-\frac{T}{2}s}\text{\hspace{1em}(7)}$$

对于s平面特定的一点$s={\alpha }_0+j{\Omega }_0$，根据式7得：

$$z=\frac{(1+\frac{T}{2}{\alpha }_0)+j\frac{T}{2}{\Omega }_0}{(1-\frac{T}{2}{\alpha }_0)-j\frac{T}{2}{\Omega }_0}$$

其模为

$$|z{|}^2=\frac{(1+\frac{T}{2}{\alpha }_0{)}^2+(\frac{T}{2}{\Omega }_0{)}^2}{(1-\frac{T}{2}{\alpha }_0{)}^2+(\frac{T}{2}{\Omega }_0{)}^2}$$

因此，对于不同的${\alpha }_0$，有

$$|z{|}^2=\{\begin{array}{ll}>1,& {\alpha }_0>0\\ =1,& {\alpha }_0=0\\ <1,& {\alpha }_0<0\end{array}$$

因此我们可以看到，s平面上的虚轴$j\Omega$映射到z平面单位圆上，$\alpha >0$的右半平面映射到z平面的单位圆外，$\alpha <0$的左半平面映射到z平面的单位圆内。

由s平面的虚轴($s=j\Omega$)和z平面的单位圆($z=e^{jw}$)之间的映射关系，由式6不难推出$w$与$\Omega$之间的映射关系。

$$j\Omega =\frac{T}{2}(\frac{1-e^{-jw}}{1+e^{-jw}})=j\frac{T}{2}\tan (\frac{w}{2})$$

或

$$\Omega =\frac{T}{2}\tan (\frac{w}{2})\text{\hspace{1em}(8)}$$

或反之

$$w=2\arctan \frac{\Omega T}{2}\text{\hspace{1em}(9)}$$

式8和9便是模拟滤波器和数字滤波器之间的频率转换关系。注意，模拟域和数字域之间的频率映射并非线性关系，因此，我们在将数字域的截止频率映射到模拟域时，需要通过式9进行换算，这一操作被称为预畸。

4.2 模拟滤波器的性能指标

4.2.1 基本参数定义

前面已经提到过，IIR滤波器设计的原型都是模拟低通滤波器，因此此处的性能指标都是模拟低通滤波器的性能指标。

上图是一个一般的低通滤波器的频响。注意，之后的参数中，下标s代表阻带（stop），下标p代表通带（pass）。${\delta }_p$和${\delta }_s$称为波纹，通常波纹以单位为dB的峰值通带波纹${\alpha }_p$和最小阻带衰减${\alpha }_s$的的形式来给出：

$${\alpha }_p=-20\lg (1-{\delta }_p)\mathrm{d}\mathrm{B}$$

$${\alpha }_s=-20\lg ({\delta }_s)\mathrm{d}\mathrm{B}$$

滤波器的频响性能指标以$a(\Omega )$给出，称为损益函数或衰减函数，定义为增益的负值，单位为$\mathrm{d}\mathrm{B}$，即$-20\lg |H_a(j\Omega )|$。

以上性能参数还可以以归一化形式给出。

假定通带幅度最大值为1，而通带波纹记为$\frac{1}{\sqrt{1+{\epsilon }^2}}$，是通带波纹最小值。因此，最大通带增益或最小通带损益是0dB。最大阻带波纹用$\frac{1}{A}$表示，从而最小阻带衰减为$-20\lg (1/A)$。
在模拟滤波器理论中还定义了另外两个参数。第一个被称为过渡比或选择性参数，定义为通带边界频率${\Omega }_p$和阻带边界频率${\Omega }_s$之比，通常用$k$表示。

$$k=\frac{{\Omega }_p}{{\Omega }_s}$$

对于低通滤波器，$k<1$。

第二个参数称为分辨参数，用$k_1$表示

$$k_1=\frac{\epsilon }{\sqrt{A^2-1}}$$

之前我们在FIR滤波器设计中见过吉布斯现象，如果通带和阻带截止频率过于接近而且两者之间增益差过大，滤波器阶数过低，都会导致出现吉布斯现象，也就是通带和阻带内的波纹。也就是通带阻带的截止频率、增益差、通带或阻带内波纹以及滤波器阶数是成一个制约关系的。通常我们在设计IIR滤波器时，都是先确定要求的性能参数：通带内波纹大小、通带和阻带的截止频率以及阻带损益。通过以上三个指标，我们便能获知符合要求的滤波器阶数，然后从模拟域构建对应阶数的模拟滤波器。

对于常用的四种模拟滤波器原型：巴特沃兹、切比雪夫I型和II型、椭圆滤波器，其对应的阶数都有不同的计算方式。由于推导过程过于复杂，因此直接给出。注意，计算出的滤波器阶数不一定是正数，向上取整即可。

4.2.2 巴特沃兹滤波器阶数公式

$$N=\frac{\lg (1/k_1)}{\lg (1/k)}$$

巴特沃兹低通滤波器的频率响应通常如下：

4.2.3 切比雪夫I型和II型滤波器阶数公式

$$N=\frac{\mathrm{arccosh}(1/k_1)}{\mathrm{arcosh}(1/k)}$$

其中，$\mathrm{arcosh}$是反双曲余弦。

低通切比雪夫I型滤波器的频响如下：

低通切比雪夫II型滤波器的频响则是反了过来，波纹在阻带里：

4.2.4 椭圆滤波器阶数公式

$$N\approx \frac{2\lg (4/k_1)}{\lg (1/\rho )}$$

其中，$\rho$的计算方式如下：

$$\begin{gathered} k^{\prime }=\sqrt{1-k^2} \\ {\rho }_0=\frac{1-\sqrt{k^{\prime }}}{2(1+\sqrt{k^{\prime }})} \\ \rho ={\rho }_0+2({\rho }_0{)}^5+15({\rho }_0{)}^9+150({\rho }_0{)}^{13} \end{gathered}$$

低通椭圆滤波器的频响如下：

由上式，我们可以看到这几种原型模拟滤波器的频响特征，在实际应用时可根据需要选取。另外有个比较重要的特性就是不同的滤波器在达到相同的损益指标时，所需的阶数是不同的。从直观上来说，相比于巴特沃兹滤波器在通带和阻带内完全平坦，切比雪夫滤波器和椭圆滤波器的波纹则是逐渐劣化的。在牺牲了波纹特性后，它们就可以用较小的阶数达到相同的损益值，其中，椭圆滤波器所需的阶数是最小的。

5. 三种模拟滤波器原型

前面已经说了那么多，现在是时候看看三种模拟滤波器的原型了。关于如何设计模拟滤波器以及它们是怎么来的，这个不属于我们要讲的范畴，对这部分刨根问底的同学需要自己去看一下相关专业书籍。

巴特沃兹滤波器的频率响应为

$$|H_a(j\Omega ){|}^2=\frac{1}{1+(\Omega /{\Omega }_c{)}^{2N}}$$

其中，${\Omega }_c$被称为3dB截止频率，因为当$\Omega =0$时，$|H_a(j\Omega ){|}^2=1$，当$\Omega ={\Omega }_c$时，$|H_a(j\Omega ){|}^2=1/2$，由于这是幅度的平方，因此这里降低的dB值为$10\lg (1/2)\approx -3\mathrm{d}\mathrm{B}$。$N$为滤波器阶数，从这我们也可以看到阶数是如何影响滚降的。通带截止频率和阻带截止频率离得越近，所需要滤波器阶数越高。

切比雪夫I型滤波器的频率响应为

$$|H_a(j\Omega ){|}^2=\frac{1}{1+{\epsilon }^2{T}_N^2(\Omega /{\Omega }_p)}$$

其中，$T_N(\Omega )$是N阶切比雪夫多项式。

$$T_N(\Omega )=\{\begin{array}{l}\cos (N\arccos \Omega ),|\Omega |\le 1\\ \cosh (N\mathrm{arccosh}\Omega ),|\Omega |>1\end{array}$$

切比雪夫II型滤波器的频率响应为

$$|H_a(j\Omega ){|}^2=\frac{1}{1+{\epsilon }^2[\frac{T_N({\Omega }_s/{\Omega }_p)}{T_N({\Omega }_s/\Omega )}{]}^2}$$

椭圆滤波器的频率响应如下

$$|H_a(j\Omega ){|}^2=\frac{1}{1+{\epsilon }^2{R}_N^2(\Omega /{\Omega }_p)}$$

其中，$R_N(\Omega )$是满足性质$R_N(1/\Omega )=1/R_N(\Omega )$的N阶有理函数，其分子的根落在区间$0<\Omega <1$内，分母的根落在$1<\Omega <\infty$内。

6. 从模拟低通滤波器构造模拟高通、带通、带阻滤波器

我们可以简单验证，以上三种模拟滤波器都是低通滤波器。为了得到高通、带通、带阻滤波器，我们需要对低通滤波器进行转换。有一些比较简单的直觉上的方法，比如要构造一个高通滤波器，只要让其频率为低通滤波器频率的倒数即可；带通滤波器则是由一个高通滤波器和一个低通滤波器串联得到；带阻滤波器可以由一个低通滤波器和高通滤波器并联得到。

为了避免模拟原型低通滤波器$H_{LP}(s)$的拉普拉斯变量与目标滤波器$H_D(s)$之间混淆，这里使用$s$表示原型低通滤波器$H_{LP}(s)$的拉普拉斯变量，用$\hat{s}$表示目标滤波器$H_D(\hat{s})$的拉普拉斯变量，两者的频率分别用$\Omega$和$\hat{\Omega }$表示。

6.1 模拟高通滤波器设计

一个通带边界频率为${\Omega }_p$的原型低通滤波器$H_{LP}(s)$可以用谱变换

$$s=\frac{{\Omega }_p{\hat{\Omega }}_p}{\hat{s}}$$

变换成一个通带边界频率为${\hat{\Omega }}_p$的高通滤波器$H_{HP}(\hat{s})$。在虚轴上，上面的变化被简化为

$$\Omega =-\frac{{\Omega }_p{\hat{\Omega }}_p}{\hat{\Omega }}$$

上面的映射表明低通滤波器在正频率范围$0\le \Omega \le {\Omega }_p$的通带映射到高通滤波器在负频率范围$-\infty <\hat{\Omega }\le -{\hat{\Omega }}_p$内。而我们知道频谱是对称的，因此低通滤波器在负频率上的通带$-{\Omega }_p\le \Omega \le 0$会映射到${\hat{\Omega }}_p<\hat{\Omega }\le \infty$内。阻带也是同理。

6.2 模拟带通滤波器设计

模拟带通滤波器变换如下：

$$s={\Omega }_p\frac{{\hat{s}}^2+{\hat{\Omega }}_o^2}{\hat{s}({\hat{\Omega }}_{p2}-{\hat{\Omega }}_{p1})}$$

在虚轴上，简化为

$$\Omega =-{\Omega }_p\frac{{\hat{\Omega }}_o^2-{\hat{\Omega }}^2}{\hat{\Omega }({\hat{\Omega }}_{p2}-{\hat{\Omega }}_{p1})}$$

其中，${\hat{\Omega }}_o$称为带通滤波器的通带中心频率。由上式可知，$\Omega =0$映射到${\hat{\Omega }}_o$。${\hat{\Omega }}_{p2}$是带阻滤波器的上通带边界频率，${\hat{\Omega }}_{p1}$是带阻滤波器的下通带边界频率。

6.3 模拟带阻滤波器设计

模拟带阻滤波器的变换如下：

$$s={\Omega }_s\frac{\hat{s}({\hat{\Omega }}_{s2}-{\hat{\Omega }}_{s2})}{{\hat{s}}^2+{\hat{\Omega }}_o^2}$$

在虚轴上简化为

$$\Omega ={\Omega }_s\frac{\hat{\Omega }({\hat{\Omega }}_{s2}-{\hat{\Omega }}_{s1})}{{\hat{\Omega }}_o^2-{\hat{\Omega }}^2}$$

其中，${\hat{\Omega }}_{s2}$为上阻带边界频率，${\hat{\Omega }}_{s1}$为下阻带边界频率。

以上是从模拟域进行变换的，实际上也可以从数字域进行变换。这里不展开讲，因为模拟域变换用得比较多。

7. IIR滤波器设计步骤总结

以上基本就是设计IIR滤波器的基本要素，这里将步骤总结一下：

1. 确定数字滤波器的性能指标，比如通带阻带频率、损益等。
2. 由第四节双线性变换法我们知道，数字域到模拟域的频率映射不是线性的，因此，我们要通过式8将从数字域映射到模拟域，获得同类型的模拟滤波器的频率指标。这一步称为预畸。
3. 通过第六节的变换关系，将目标模拟滤波器的指标转换为低通原型滤波器的指标。
4. 通过第四节的性能指标，计算出模拟原型低通滤波器的阶数。
5. 通过第五节讲的几种原型低通滤波器的频率响应，设计出符合性能要求的低通滤波器。
6. 通过第六节的变换关系，将设计好的模拟原型低通滤波器变换为目标类型的模拟滤波器。
7. 通过式9，将目标类型的模拟滤波器转换到数字域，即得到了目标数字滤波器。

8. 使用Matlab设计IIR滤波器

这里以切比雪夫II型滤波器为例。

假设我们的采样率是48000Hz，那么截止频率就是24000Hz。要设计一个通带边界频率为1000，阻带边界频率为1500的低通滤波器，阻带损益为40dB，代码如下：

clear all
close all

sampleRate = 48000;
maxFreq = sampleRate / 2;

Wp = 1000;
Ws = 1500;

lossGain = 40;

% 使用cheby2ord计算阶数和被重新调整过的边界频率。注意
% 这里频率参数统一是归一化的。第三个参数代表通带波纹。
[N, Wn] = cheb2ord(Wp / maxFreq, Ws / maxFreq, 1, lossGain);

% 使用cheby2计算滤波器系数。它是形如
% y[n]=b[0]x[n]+b[1]x[n−1]+⋯+b[N]x[n−N]+
% a[1]y[n−1]+a[2]y[n−2]+⋯+a[M]y[n−M]
% 的系数。
[b, a] = cheby2(N, lossGain, Wn, 'low');

freqz(b, a);

其频率响应和相位响应为

可见其完全符合预期，且相位响应不是线性的。

打印各数据如下：

N =

     7


Wn =

    0.0625


b =

    0.0051   -0.0238    0.0411   -0.0224   -0.0224    0.0411   -0.0238    0.0051


a =

    1.0000   -6.2280   16.6627  -24.8213   22.2308  -11.9698    3.5872   -0.4615

a为y的系数，b为x的系数，可以看到a的第一项实际上总是1。

9. 有限字长效应

IIR滤波器有个非常显著的缺点就是它的稳定性较差，FIR滤波器可以设计出任何阶数，但是IIR的滤波器阶数如果过高，就会导致不稳定。这是因为在程序中，数据类型都以二进制进行存储，而二进制存储浮点数，必然会有误差存在。因此，IIR滤波器的阶数不宜过高。
通过另一种方法可以规避有限字长效应，那就是将高阶IIR滤波器拆分成多个低阶IIR滤波器的串联。这一点从式5可以看出来。这样以来，就可以降低误差带来的影响。在matlab中，tf2sos函数可以将高阶传输函数转换为2阶传输函数。

其调用方式为

[sos, g] = tf2sos(b, a)

其中，g就是式5中的系数g。
sos的形式如下以及对应传输函数的形态如下，更详细的用法可以查看matlab的帮助文档。