数字信号处理(六)IIR数字滤波器的设计

数字滤波器

什么是数字滤波器？

指输入输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。是离散系统的通用名称。

数字滤波器的分类：

经典滤波器，特点是其输入信号中有用的频率成分和期望滤除的频率成分各占不同的频带，通过一个合适的选频滤波器滤除某个频带或频率成分的干扰，得到纯净信号，达到滤波的目的

现代滤波器，用于信号与干扰的频带互相重叠的情况，现代滤波器是根据随机信号的一些统计特性，在某种最佳准则下，最大限度的抑制干扰，恢复信号，从而达到滤波的目的，如维纳滤波器、卡尔曼滤波器等

经典数字滤波器从滤波特性上可分为低通、高通、带通、带阻和全通

数字滤波器从单位脉冲响应长度分类，可分为无限长单位脉冲响应(IIR)滤波器和有限长单位脉冲响应(FIR)滤波器

离散线性时不变系统的差分方程为：$y(n)={\sum }_{i=0}^M{b}_{i}x(n-i)-{\sum }_{k=1}^N{a}_{k}y(n-k)$，那么其系统函数可表示为：$H(z)=\frac{{\sum }_{i=0}^M{b}_i{z}^{-i}}{1+{\sum }_{k=1}^N{a}_k{z}^{-k}}$

• 当N≥1，$a_k$中至少有一个非零系数时，该系统中存在反馈回路，其所对应的滤波器称作无限长脉冲响应滤波器（infinite impulse response filter），简称IIR滤波器。N是IIR滤波器的阶数，表示系统中反馈环的个数。一般假定IIR数字滤波器满足M≤N，这时将系统称为N阶的IIR数字滤波器。
• 当$a_k$均为零系数时，其所对应的滤波器称作有限长脉冲响应滤波器（finite impulse response filter），简称FIR滤波器

数字滤波器技术指标

常用的数字滤波器一般属于选频滤波器，假设数字滤波器的频率响应函数$H(e^{jw})$用下式表示：
$$H(e^{jw})=|H(e^{jw})|e^{j\theta (w)}$$

• $|H(e^{jw})|$称为幅频特性函数，表示信号通过该滤波器后各频率成分振幅衰减情况
• $\theta (\omega )$为相频特性函数，反映各频率成分通过滤波器后在时间上的延时情况
• 一般选频滤波器的技术指标由幅频特性给出

实际滤波器的通带和阻带中都允许有一定的误差容限，即通带不是完全平的，阻带不是绝对衰减到零，通带和阻带之间还有一定的过渡带

数字低通滤波器的幅频响应曲线

通带内和阻带内允许的衰减一般用dB数表示，归一化后通带衰减${\alpha }_p$和阻带衰减${\alpha }_s$分别定义为：
$${\alpha }_p=-20lg|H_a(j{\omega }_p)|dB$$

$${\alpha }_s=-20lg|H_a(j{\omega }_s)|dB$$

当幅度下降到$\sqrt{2}/2$时，$\omega ={\omega }_c$，此时${\alpha }_p=3dB$，称为3dB通带截止频率

${\omega }_p,{\omega }_s,{\omega }_c$统称为边界频率，它们是数字滤波器的重要参数

IIR滤波器设计方法

直接法： 直接在频域或者时域中设计数字滤波器

间接法： 借助于模拟滤波器的设计原型进行设计的，又称滤波器的函数模型设计法，较常用。

IIR滤波器的函数模型设计法(间接法)

在IIR滤波器设计过程中，通常利用模拟滤波器来设计数字滤波器

借助模拟滤波器的设计方法的基本步骤：

• 将数字滤波器的技术指标转换成模拟滤波器的技术指标

• 按转换后技术指标、设计模拟低通滤波器的$H_a(s)$

• 然后将$H_a(s)$按某种方法转换成数字滤波器的系统函数$H(z)$

• 如果不是低通，则必须先将其转换成低通模拟滤波器的技术指标

• 设计过程中用到的两种变换：

  ①频带变换： 将低通转换为高通、带通、带阻滤波器

  ②变换域变换： 将模拟滤波器转换成数字滤波器，从s->z或者$H_a(s)->H(z)$的变换

IIR滤波器的间接设计法是以模拟滤波器为基础的，有若干典型的模拟滤波器原型可供选择：

• 巴特沃思滤波器：具有单调下降的幅频特性
• 切比雪夫滤波器：幅频特性在通带或者阻带有等波纹特性，可以提高选择性
• 椭圆滤波器：选择性相对前两种是最好的，单通带和阻带内均呈现等波纹幅频特性，相位特性的非线性也稍严重
• 贝塞尔滤波器

模拟低通滤波器的技术指标

模拟低通滤波器的技术指标有${\alpha }_p,{\Omega }_p,{\alpha }_s,{\Omega }_s$，其中${\Omega }_p和{\Omega }_s$分别称为通带截止频率和阻带截止频率，${\alpha }_p$是通带内的最大允许衰减，${\alpha }_s$是阻带内最小允许衰减。归一化的${\alpha }_p,{\alpha }_s$分别表示如下：
$${\alpha }_p=-10lg|H_a(j{\Omega }_p){|}^{2}dB$$

$${\alpha }_s=-10lg|H_a(j{\Omega }_s){|}^{2}dB$$

$H_a(j{\Omega }_c)=\sqrt{2}/2$时，$-10lg|H_a(j{\Omega }_c){|}^2=3dB$，故称${\Omega }_c$为模拟低通滤波器的3dB截止频率

模拟低通滤波器的上述技术指标给定后，需要设计滤波器的系统函数$H_a(s)$，具体步骤如下：

模拟滤波器原型介绍

1、巴特沃斯模拟低通滤波器

介绍

特点：通带和阻带的幅度响应都是平的

N阶巴特沃斯(Butterworth)模拟低通滤波器的幅频响应的平方函数为：
$$|H_a(j\Omega ){|}^2=\frac{1}{1+(\Omega /{\Omega }_c{)}^{2N}}$$
其中${\Omega }_C$是巴特沃斯低通滤波器的3dB截止频率
$$|H_a(j0){|}^2=1,|H_a(j\infty ){|}^2=0$$

巴特沃斯模拟低通滤波器的性质：

①幅频响应平方函数$|H_a(j\Omega ){|}^2$是 $\Omega$的单调减函数 ：
$$\frac{d|H_a(j\Omega ){|}^2}{d\Omega }=-\frac{2N(\Omega /{\Omega }_c{)}^{2N-1}}{[1+(\Omega /{\Omega }_c{)}^{2N}{]}^2}<0$$
②对于所有的N，当$\Omega =0$时，$|H_a(j\Omega ){|}^2=1$

③$\Omega ={\Omega }_c$时，$|H_a(j\Omega )|=\sqrt{2}/2$，即在$|H_a(j\Omega )|=1$处有3dB的衰减

④$\Omega >{\Omega }_c$时，随$\Omega$增大，$|H_a(j\Omega )|$迅速下降，${\Omega }_c$衰减的速度与阶数N有关，N愈大，衰减速度愉快，过渡带愈窄，N->$\infty$时，趋向于理想低通滤波器

⑤在$\Omega =0$处各阶导数存在且等于0，因此在该点得到最大值

巴特沃斯滤波器的系统函数：

$\because$
$$H_a(s)H_a(-s)=|H_a(j\Omega ){|}^2=\frac{1}{1+(\Omega /{\Omega }_c{)}^{2N}}{|}_{\Omega =s/j}=\frac{1}{1+(\frac{s}{j{\Omega }_c}{)}^{2N}}=\frac{(j{\Omega }_c{)}^{2N}}{s^{2N}+(j{\Omega }_c{)}^{2N}}$$
$\therefore$

​ $H_a(s)H_a(-s)$的极点为$p_k={\Omega }_c[e^{j\frac{\pi }{2N}(2k+N-1)}]，k=1,2,...,2N$

由于$H_a(s)$的因果稳定性，要求极点必须落在s平面左半平面，即极点满足$\frac{\pi }{2}<\pi (\frac{2k+N-1}{2N})<\frac{3\pi }{2}$

化简得：$\frac{1}{2}<k<N+\frac{1}{2}$，即k=1,2,…,N

因此，N阶巴特沃斯模拟低通滤波器N个极点
$$p_k={\Omega }_c[e^{j\frac{\pi }{2N}(2k+N-1)}]，k=1,2,...,N$$
其系统函数为
$$H_a(s)=\frac{{\Omega }_c^N}{{\prod }_{k=1}^N(s-p_k)}$$

巴特沃斯滤波器的设计步骤

1、根据模拟滤波器的设计指标${\alpha }_p,{\Omega }_p,{\alpha }_s,{\Omega }_s$确定滤波器的阶数N
$$N=\frac{lg[(10^{0.1{\alpha }_p}-1)/(10^{0.1{\alpha }_s}-1)]}{2lg({\Omega }_p/{\Omega }_s)}$$
取N为比计算结果大的最小整数，就是巴特沃斯模拟低通滤波器的阶数

2、确定滤波器的3dB截止频率
$${\Omega }_{cp}=\frac{{\Omega }_p}{\sqrt[2N]{10^{0.1{\alpha }_p}-1}},{\Omega }_{cs}=\frac{{\Omega }_s}{\sqrt[2N]{10^{0.1{\alpha }_s}-1}}$$
实际设计时，${\Omega }_c$可在${\Omega }_{cp}\le {\Omega }_c\le {\Omega }_{cs}$范围内选择

3、求出N个极点，$p_k={\Omega }_c[e^{j\frac{\pi }{2N}(2k+N-1)}]，k=1,2,...,N$，得到滤波器的系统函数$H_a(s)$。

2、切比雪夫低通滤波器

切比雪夫模拟低通滤波器的幅频响应在一个频带中具有等波纹特性

• $Chebyshev-Ifilters$ 通带等波纹响应
• $Chebyshev-IIfilters$ 阻带等波纹响应

相同指标下，切比雪夫滤波器比巴特沃斯滤波器阶数低

介绍

切比雪夫I型滤波器的平方幅度响应为：
$$|H_a(j\Omega ){|}^2=\frac{1}{1+{ϵ}^2{C}_N^2(\Omega /{\Omega }_p)}$$
$ϵ$是小于1的正数，表示滤波器通带内幅度波动的程度，$ϵ$愈大，幅度波动也愈大

${\Omega }_p$是通带截止频率，$\lambda =\Omega /{\Omega }_p$称为对${\Omega }_p$的归一化频率

$C_N(x)$是N阶切比雪夫多项式：
$$C_N(x)=\{\begin{array}{ll}cos[Narccos(x)]& |x|\le 1\\ cosh[Narccosh(x)]& |x|>1\end{array}$$

• 当$|x|\le 1$时，$|C_N(x)|\le 1$。$C_N(x)$在-1和+1之间振荡，振荡的次数与N成正比，具有等波纹性
• 当$|x|\le 1$时，${ϵ}^2{C}_N^2(x)$在0到${ϵ}^2$之间波动，函数$1+{ϵ}^2{C}_N^2(\Omega /{\Omega }_p)$的倒数即是平方幅度响应$|H_a(j\Omega ){|}^2$，其在$[0,{\Omega }_p]$上$x=(\Omega /{\Omega }_p)<1$，因此具有等波纹波动，最大值为1，最小值为$\frac{1}{1+{ϵ}^2}$
• 当$|x|>1$时，$|C_N(x)|>1$，且$|C_N(x)|$随着$|x|$的增加单调上升。在$[{\Omega }_p,\infty )$上$x=(\Omega /{\Omega }_p)>1$，$1+{ϵ}^2{C}_N^2(\Omega /{\Omega }_p)$单调上升，因此$|H_a(j\Omega ){|}^2$单调减小
• N相同时，切比雪夫滤波器比巴特沃斯滤波器具有较窄的过渡带

切比雪夫滤波器的设计步骤

1、确定模拟滤波器的设计指标${\Omega }_p,{\Omega }_s,{\alpha }_p,{\alpha }_s$

2、确定参数$ϵ$
$${ϵ}^2=10^{0.1{\alpha }_p}-1$$
3、求出滤波器的阶数N，取N为比计算结果大的最小整数，就是切比雪夫模拟低通滤波器的阶数
$$N=arccosh(\frac{10^{0.1{\alpha }_s}-1}{10^{0.1{\alpha }_p}-1}{)}^{1/2}/arccosh({\Omega }_s/{\Omega }_p)$$
4、由下式求归一化切比雪夫I型低通滤波器的极点
$$p_{nk}={\sigma }_k+j{\Omega }_k,k=1,2,...,N,其中$$

$${\sigma }_k=-sinh(\beta )sin\frac{(2k-1)\pi }{2N},{\Omega }_k=-cosh(\beta )cos\frac{(2k-1)\pi }{2N},\beta =\frac{arsinh(1/ϵ)}{N}$$

5、求归一化切比雪夫I型低通滤波器的系统函数$H_n(s)$：

N为奇数时，$H_n(s)=\frac{sinh\beta }{1+sinh\beta }{\prod }_{k=1}^{N/2}\frac{{\sigma }_k^2+{\Omega }_k^2}{s^2-2{\sigma }_{k}s+({\sigma }_k^2+{\Omega }_k^2)}$

N为偶数时，$H_n(s)=\frac{1}{\sqrt{1+{ϵ}^2}}{\prod }_{k=1}^{N/2}\frac{{\sigma }_k^2+{\Omega }_k^2}{s^2-2{\sigma }_{k}s+({\sigma }_k^2+{\Omega }_k^2)}$

3、椭圆模拟低通滤波器

介绍

通带和阻带均具有等波纹响应，幅频响应特性与FIR等波纹滤波器类似

幅频响应平方函数$|H_a(j\Omega ){|}^2=\frac{1}{1+{ϵ}^2{U}_N^2(\Omega /{\Omega }_c)}$

其中$U_N(.)$是N阶雅克比椭圆函数

对比

一、当阶数相同时，对相同的通带最大衰减$a_p$，和阻带最小衰减$a_s$：

1、巴特沃斯滤波器具有单调下降的幅频特性，过渡带最宽。两种类型的切比雪夫滤波器的过渡带宽度相等，比巴特沃斯滤波器的过渡带窄，但比椭圆滤波器的过渡带宽

2、切比雪夫I型滤波器在通带具有等波纹幅频特性，过渡带和阻带是单调下降的幅频特性

3、切比雪夫II型滤波器的通带幅频响应几乎与巴特沃斯滤波器相同，阻带是等波纹幅频特性

4、椭圆滤波器的过渡带最窄，通带和阻带均是等波纹幅频特性

二、复杂性：

在满足相同的滤波器幅频响应指标条件下，巴特沃斯滤波器阶数最高，椭圆滤波器的阶数最低，而且阶数差别较大。所以，就满足滤波器幅频响应指标而言，椭圆滤波器的性价比最高，应用较广泛。

模拟到数字滤波器的转换

模拟滤波器$H_a(s)$到数字滤波器$H(z)$的转换通过s到z的复值映射来实现，需要满足：

• 因果稳定的模拟滤波器变成数字滤波器，仍应是因果稳定的，保证s平面的左半平面要映射到z平面的单位圆内部
• 数字滤波器的频率响应 应模仿模拟滤波器的频率特性，即s平面的虚轴映射为z平面的单位圆，相应的频率之间呈线性关系

$H_a(s)$从s平面转换到z平面的典型方法有：脉冲响应不变法、双线性变换法

脉冲响应不变法

基本原理： 使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应$h_a(t)$

具体步骤：

以理想采样信号${\hat{h}}_a(t)$为桥梁，推导从模拟滤波器转换到数字滤波器时，s平面和z平面之间的映射关系：

设${\hat{h}}_a(t)={\sum }_{n=-\infty }^{\infty }{h}_a(t)\delta (t-nT)$，拉氏变换得到：
$${\hat{H}}_a(s)={\int }_{-\infty }^{\infty }{\hat{h}}_a(t)e^{-st}dt={\int }_{-\infty }^{\infty }[\sum _{n=-\infty }^{\infty }{h}_a(t)\delta (t-nT)]e^{-st}dt=\sum _n{h}_a(nT)e^{-snT}$$
与${\hat{H}}_a(s)={\sum }_{n}h(n)e^{-snT}={\sum }_{n}h(n)z^{-n}{|}_{z=e^{sT}}=H(z){|}_{z=e^{sT}}$比较得：

理想采样信号的拉氏变换与相应的采样序列的Z变换之间的映射关系为：$z=e^{sT}$，即脉冲响应不变法对应的s和z平面之间的映射关系。

利用$H_a(s)$直接变换成$H(z)$：
$$H(z)=H_a(s){|}_{\frac{1}{s-s_k}=\frac{1}{1-e^{s_{k}T}{z}^{-1}}},k=1,2,...,N,其中s_k是模拟低通滤波器H_a(s)的极点$$
脉冲响应不变法设计IIT数字滤波器的步骤：

1、确定采样间隔T

2、根据采样间隔T，将给定的数字滤波器频率指标转换为模拟滤波器的频率指标：${\Omega }_p={\omega }_p/T,{\Omega }_s={\omega }_s/T$

3、根据指标${\alpha }_p,{\Omega }_p,{\alpha }_s,{\Omega }_s$，设计模拟滤波器$H_a(s)$

4、用脉冲响应不变法，将模拟滤波器$H_a(s)$变换为数字滤波器$H(z)$：
$$H(z)=H_a(s){|}_{\frac{1}{s-s_k}=\frac{1}{1-e^{s_{k}T}{z}^{-1}}},k=1,2,...,N,其中s_k是模拟低通滤波器H_a(s)的极点$$

双线性变换法

脉冲响应不变法是从s平面到z平面的多值映射，会产生频谱混叠。

双线性变换法采用非线性频率压缩的方法，将整个频率轴压缩到$\pm \pi /T$之间，再用$z=e^{sT}$转换到z平面上，由于s平面与z平面是一一对应的关系，从而有效地避免了脉冲响应不变法中多值映射所引起的频谱混叠现象；双线性变换法将s的左半平面映射到z平面的单位圆内，因此也满足了对系统因果稳定性的要求。

双线性变换法介绍：

设模拟滤波器的系统函数是$H_a(s),s=j\Omega$，经过非线性频率压缩后用${\hat{H}}_a(s_1),s_1=j{\Omega }_1$表示

用正切变换实现频率压缩：$\Omega =\frac{2}{T}tan(\frac{1}{2}{\Omega }_{1}T)$，T是采样间隔

当${\Omega }_1从-\pi /T$经过0变化到$\pi /T$时，$\Omega$则由$-\infty$经过0变化到+$\infty$，实现了s平面上整个虚轴完全压缩到$s_1$平面上虚轴的$\pm \pi /T$之间的转换

用正切变换实现频率压缩：$\Omega =\frac{2}{T}tan(\frac{1}{2}{\Omega }_{1}T)$：

根据欧拉公式的$sin\theta =\frac{e^{j\theta }-e^{-j\theta }}{2j},cos\theta =\frac{e^{j\theta }+e^{-j\theta }}{2}$

所以$jtan\theta =\frac{e^{j\theta }-e^{-j\theta }}{e^{j\theta }+e^{-j\theta }}$

因此
$$j\Omega =\frac{2}{T}jtan(\frac{1}{2}{\Omega }_{1}T)=\frac{2}{T}\frac{e^{j{\Omega }_{1}T/2}-e^{-j{\Omega }_{1}T/2}}{e^{j{\Omega }_{1}T/2}+e^{-j{\Omega }_{1}T/2}}=\frac{2}{T}\frac{1-e^{-j{\Omega }_{1}T}}{1+e^{-j{\Omega }_{1}T}}$$

$$s=j\Omega =\frac{2}{T}\frac{1-e^{-s_{1}T}}{1+e^{-s_{1}T}}=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}},其中z=e^{s_{1}T}$$

展开得$\frac{T}{2}sz+\frac{T}{2}s-z+1=0$

如上图所示，由s平面到$s_1$平面的非线性频率压缩，使${\hat{H}}_a(s_1)$带限于$\pi /T$，因此，再用脉冲响应不变法从$s_1$平面转换到z平面不可能产生频谱混叠现象

因此，模拟系统$H_a(s)$和对应的数字系统H(z)的关系为：$H(z)=H_a(s){|}_{s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}$

模拟频率$\Omega$与数字频率$\omega$之间的关系：$\Omega =\frac{2}{T}tan\frac{\omega }{2}$， 即s平面上的模拟频率$\Omega$与z平面上的数字频率$\omega$成非线性正切关系

双线性变换法设计IIR数字滤波器的步骤：

1、确定参数T。为简单起见，可取T=2。通常取${\Omega }_{c}T=1$

2、确定数字低通滤波器的技术指标：通带截止频率${\omega }_p$、通带衰减${\alpha }_p$、阻带截止频率${\omega }_p$、阻带衰减${\alpha }_p$。采用双线性变换法模数频率间的转换关系，将数字低通滤波器的指标转换成模拟低通滤波器的技术指标：$\Omega =\frac{2}{T}tan(\frac{1}{2}\omega )$

3、按照模拟低通滤波器的技术指标设计模拟低通滤波器

4、将模拟滤波器$H_a(s)$，通过下式从s平面转换到z平面，得到数字低通滤波器系统函数H(z)：
$$H(z)=H_a(s){|}_{s=\frac{2}{T}.\frac{1-z^{-1}}{1+z^{-1}}}$$

数字滤波器的频带变换

IIR滤波器间接设计法的两种思路：

从原型低通数字滤波器到各种类型的数字滤波器的频带变换也称作Z平面变换法

设原型低通数字滤波器的系统函数为$H_L(z)$，所需类型数字滤波器的系统函数为$H(z)$，这种变换是将$H_L(z)$的$z_L$平面映射到$H(z)$的z平面

用$z_L和z$来表示变换前后的z平面，假设从$z_L$平面到z平面的映射关系为：
$$z_L^{-1}=G(z^{-1}),G(Z^{-1})必须是z^{-1}的有理函数，这样H(z)才是可实现的$$

$$\begin{gathered} H(z)=H_L(z){|}_{z_L^{-1}=G(z^{-1})},其中，z_L平面的单位圆必须映射到z平面的单位圆上 \\ {z}_L平面的单位圆内部必须映射到z平面的单位圆内部 \end{gathered}$$

如果能找到适当的从单位圆周到单位圆周的数字变换把原型低通数字的通带频带变换成其他类型数字滤波器所要求的通带频带，就等于完成了这种滤波器的设计

设$w_L和w$分别是$z_L和z$的频率变量，则在各自单位圆$z_L=e^{j{w}_L}和z=e^{jw}$上有：
$$|z_L^{-1}|=|G(z^{-1})|=|G(e^{-jw})|=1，z_L是单位圆$$
因此，函数G的幅频特性对所有的频率都为1，是一个全通函数

全通滤波器一般形式为：
$$\begin{gathered} z_L^{-1}=G(z^{-1})=\pm \prod _{k=1}^N\frac{z^{-1}-p_k^{\ast }}{1-p_k{z}^{-1}}， \\ 其中N全通函数的阶数，p_k是G(z^{-1})的极点且必须在单位圆内， \end{gathered}$$
选择合适的N和$p_k$就可以得到各种类型的变换