无限脉冲响应数字滤波器的设计

无限脉冲响应数字滤波器的设计第一页，共八十六页，2022年，8月28日

6.1数字滤波器的基本概念

1.数字滤波器的分类数字滤波器从实现的网络结构或者从单位脉冲响应分类，可以分成无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。它们的系统函数分别为：第二页，共八十六页，2022年，8月28日理想低通、高通、带通、带阻滤波器幅度特性第三页，共八十六页，2022年，8月28日2数字滤波器的技术要求我们通常用的数字滤波器一般属于选频滤波器。假设数字滤波器的传输函数H(ejω)用下式表示：第四页，共八十六页，2022年，8月28日如图所示。滤波器的主要性能指标（以低通滤波器为例）如下：（1）频带：包括通带、过渡带、阻带。（2）误差：通带和阻带误差过渡带通带阻带第五页，共八十六页，2022年，8月28日3.数字滤波器设计方法概述

IIR滤波器和FIR滤波器的设计方法是很不相同的。IIR滤波器设计方法有两类，经常用的一类设计方法是借助于模拟滤波器的设计方法进行的。其设计步骤是：先设计模拟滤波器得到传输函数Ha(s)，然后将Ha(s)按某种方法转换成数字滤波器的系统函数H(z)。第六页，共八十六页，2022年，8月28日

6.2模拟滤波器的设计

模拟滤波器的理论和设计方法已发展得相当成熟，且有若干典型的模拟滤波器供我们选择，如巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用。第七页，共八十六页，2022年，8月28日1.巴特沃斯低通滤波器的设计方法巴特沃斯低通滤波器的幅度平方函数|Ha(jΩ)|2用下式表示：(6.2.6)第八页，共八十六页，2022年，8月28日第九页，共八十六页，2022年，8月28日2、切比雪夫滤波器切比雪夫滤波器有两种类型：切比雪夫Ⅰ型滤波器的特点是在通带内有等波纹变化，阻带内单调下降。切比雪夫Ⅱ滤波器的特点是在通带内单调下降，阻带内有等波纹变化。第十页，共八十六页，2022年，8月28日波纹系数为0.5dB的切比雪夫Ⅰ低通滤波器的幅频响应平方函数第十一页，共八十六页，2022年，8月28日阻带衰减10dB的切比雪夫Ⅱ低通滤波器的幅频响应平方函数第十二页，共八十六页，2022年，8月28日切比雪夫Ⅰ滤波器与巴特沃兹滤波器过渡带的比较第十三页，共八十六页，2022年，8月28日在具有相同阶数

的情况下，切比雪夫滤波器与巴特沃兹滤波器相比，具有较窄的过渡带第十四页，共八十六页，2022年，8月28日3、椭圆滤波器椭圆滤波器的特点是幅度频率响应在通带和阻带内均为等波纹的。且与上述两种滤波器相比，过渡带的下降斜度更大。一般来说，对于指定的滤波器指标（阶数与波纹），椭圆滤波器能以最低的阶数实现。椭圆滤波器的幅频响应平方函数如下：Jacobian椭圆函数第十五页，共八十六页，2022年，8月28日通带波纹0.3dB、阻带波纹10dB的椭圆低通滤波器幅频响应平方函数第十六页，共八十六页，2022年，8月28日6.3用脉冲响应不变法设计IIR

数字低通滤波器

为了保证转换后的H(z)稳定且满足技术要求，对转换关系提出两点要求：

(1)因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的。

(2)数字滤波器的频率响应模仿模拟滤波器的频响，s平面的虚轴映射z平面的单位圆，相应的频率之间成线性关系。第十七页，共八十六页，2022年，8月28日

设模拟滤波器的传输函数为Ha(s),相应的单位冲激响应是ha(t)

设模拟滤波器Ha(s)只有单阶极点，且分母多项式的阶次高于分子多项式的阶次，将Ha(s)用部分分式表示：(6.3.1)

式中si为Ha(s)的单阶极点。将Ha(s)进行逆拉氏变换得到ha(t)：(6.3.2)一、由Ha(s)求相应的H(z)第十八页，共八十六页，2022年，8月28日

式中u(t)是单位阶跃函数。对ha(t)进行等间隔采样，采样间隔为T，得到：(6.3.3)对上式进行Z变换，得到数字滤波器的系统函数H(z)：(6.3.4)第十九页，共八十六页，2022年，8月28日二、S平面和Z平面之间的映射第二十页，共八十六页，2022年，8月28日

上式表示采样信号的拉氏变换与相应的序列的Z变换之间的映射关系可用下式表示：(6.3.6)第二十一页，共八十六页，2022年，8月28日(6.3.6)式可称为标准映射关系。下面进一步分析这种映射关系。设按照(6.3.6)式，得到:因此得到:(6.3.10)第二十二页，共八十六页，2022年，8月28日我们由三方面说明这个关系

σ=0,r=1σ<0,r<1σ>0,r>1

另外，注意到z=esT是一个周期函数，可写成为任意整数第二十三页，共八十六页，2022年，8月28日z=esT,s平面与z平面之间的映射关系第二十四页，共八十六页，2022年，8月28日图6.3.2脉冲响应不变法的频率混叠现象第二十五页，共八十六页，2022年，8月28日

例若模拟滤波器的系统函数为试用脉冲响应不变法求相应的数字滤波器的系统函数。（设采样周期为

）第二十六页，共八十六页，2022年，8月28日将系统函数用部分分式法展开模拟滤波器的极点S1=-0.5,S2=-2因此，相应数字滤波器的系统函数为第二十七页，共八十六页，2022年，8月28日

6.4用双线性变换法设计IIR数字

低通滤波器

正切变换实现频率压缩：(6.4.1)

式中T仍是采样间隔，当Ω1从π/T经过0变化到π/T时，Ω则由∞经过0变化到+∞，实现了s平面上整个虚轴完全压缩到s1平面上虚轴的±π/T之间的转换。这样便有(6.4.2)第二十八页，共八十六页，2022年，8月28日

再通过转换到z平面上，得到：(6.4.3)(6.4.4)第二十九页，共八十六页，2022年，8月28日

下面分析模拟频率Ω和数字频率ω之间的关系。图6.4.1双线性变换法的映射关系第三十页，共八十六页，2022年，8月28日

令s=jΩ,z=ejω，并代入(6.4.3)式中，有(6.4.5)图6.4.2双线性变换法的频率变换关系第三十一页，共八十六页，2022年，8月28日图6.4.3双线性变换法幅度和相位特性的非线性映射第三十二页，共八十六页，2022年，8月28日

设第三十三页，共八十六页，2022年，8月28日

表6.4.1系数关系表第三十四页，共八十六页，2022年，8月28日第三十五页，共八十六页，2022年，8月28日

例试分别用脉冲响应不变法和双线性不变法将图所示的RC低通滤波器转换成数字滤波器。解首先按照图写出该滤波器的传输函数Ha(s)为

利用脉冲响应不变法转换，数字滤波器的系统函数H1(z)为第三十六页，共八十六页，2022年，8月28日

利用双线性变换法转换，数字滤波器的系统函数H2(z)为H1(z)和H2(z)的网络结构分别如图6.4.5(a),(b)所示。图6.4.5例图——H1(z)和H2(z)的网络结构

(a)H1(z);(b)H2(z)第三十七页，共八十六页，2022年，8月28日

下面我们总结利用模拟滤波器设计IIR数字低通滤波器的步骤。

(1)确定数字低通滤波器的技术指标：通带截止频率ωp、通带衰减αp、阻带截止频率ωs、阻带衰减αs。

(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。

如果采用双线性变换法，边界频率的转换关系为第三十八页，共八十六页，2022年，8月28日图例6.4.1图——数字滤波器H1(z)和H2(z)的幅频特性第三十九页，共八十六页，2022年，8月28日(3)按照模拟低通滤波器的技术指标设计模拟低通滤波器。

(4)将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。例6.4.2设计低通数字滤波器，要求在通带内频率低于0.2πrad时，容许幅度误差在1dB以内；在频率0.3π到π之间的阻带衰减大于15dB。指定模拟滤波器采用巴特沃斯低通滤波器。试分别用脉冲响应不变法和双线性变换法设计滤波器。第四十页，共八十六页，2022年，8月28日

解

(1)用脉冲响应不变法设计数字低通滤波器。①数字低通的技术指标为

ωp=0.2πrad,αp=1dB;ωs=0.3πrad,αs=15dB②模拟低通的技术指标为

T=1s,Ωp=0.2πrad/s,αp=1dB;Ωs=0.3πrad/s,αs=15dB第四十一页，共八十六页，2022年，8月28日③设计巴特沃斯低通滤波器。先计算阶数N及3dB截止频率Ωc。第四十二页，共八十六页，2022年，8月28日

取N=6。为求3dB截止频率Ωc,将Ωp和αp代入(6.2.17)式，得到Ωc=0.7032rad/s，显然此值满足通带技术要求，同时给阻带衰减留一定余量，这对防止频率混叠有一定好处。根据阶数N=6，查表，得到归一化传输函数为

为去归一化，将p=s/Ωc代入Ha(p)中，得到实际的传输函数Ha(s),第四十三页，共八十六页，2022年，8月28日④用脉冲响应不变法将Ha(s)转换成H(z)。首先将Ha(s)进行部分分式，并按照(6.3.11)式、(6.3.12)式，或者(6.3.13)式和(6.3.14)式，得到：第四十四页，共八十六页，2022年，8月28日图6.4.7例图——用脉冲响应不变法设计的数字低通滤波器的幅度特性第四十五页，共八十六页，2022年，8月28日(2)用双线性变换法设计数字低通滤波器。①数字低通技术指标仍为

ωp=0.2πrad,αp=1dB;ωs=0.3πrad,αs=15dB②模拟低通的技术指标为第四十六页，共八十六页，2022年，8月28日③设计巴特沃斯低通滤波器。阶数N计算如下：

取N=6。为求Ωc,将Ωs和αs代入(6.2.18)式中，得到Ωc=0.7662rad/s。这样阻带技术指标满足要求，通带指标已经超过。第四十七页，共八十六页，2022年，8月28日

根据N=6，查表得到的归一化传输函数Ha(p)与脉冲响应不变法得到的相同。为去归一化，将p=s/Ωc代入Ha(p)，得实际的Ha(s)，④用双线性变换法将Ha(s)转换成数字滤波器H(z)：第四十八页，共八十六页，2022年，8月28日图6.4.8例图——用双线性变换法设计的数字低通滤波器的幅度特性第四十九页，共八十六页，2022年，8月28日

6.5数字高通、带通和带阻滤波器的设计

例如高通数字滤波器等。具体设计步骤如下：

(1)确定所需类型数字滤波器的技术指标。

(2)将所需类型数字滤波器的技术指标转换成所需类型模拟滤波器的技术指标，转换公式为

第五十页，共八十六页，2022年，8月28日(3)将所需类型模拟滤波器技术指标转换成模拟低通滤波器技术指标(具体转换公式参考本章6.2节)。

(4)设计模拟低通滤波器。

(5)将模拟低通通过频率变换，转换成所需类型的模拟滤波器。

(6)采用双线性变换法，将所需类型的模拟滤波器转换成所需类型的数字滤波器。第五十一页，共八十六页，2022年，8月28日

例6.5.1设计一个数字高通滤波器，要求通带截止频率ωp=0.8πrad，通带衰减不大于3

dB，阻带截止频率ωs=0.44πrad,阻带衰减不小于15dB。希望采用巴特沃斯型滤波器。解

(1)数字高通的技术指标为

ωp=0.8πrad,αp=3dB;ωs=0.44πrad,αs=15dB第五十二页，共八十六页，2022年，8月28日(2)模拟高通的技术指标计算如下：令T=1，则有(3)模拟低通滤波器的技术指标计算如下：第五十三页，共八十六页，2022年，8月28日

将Ωp和Ωs对3dB截止频率Ωc归一化，这里Ωc=Ωp,(4)设计归一化模拟低通滤波器G(p)。模拟低通滤波器的阶数N计算如下：第五十四页，共八十六页，2022年，8月28日

查表，得到归一化模拟低通传输函数G(p)为

为去归一化，将p=s/Ωc代入上式得到：(5)将模拟低通转换成模拟高通。将上式中G(s)的变量换成1/s，得到模拟高通Ha(s)：第五十五页，共八十六页，2022年，8月28日(6)用双线性变换法将模拟高通H

(s)转换成数字高通H(z)：实际上(5)、(6)两步可合并成一步，即第五十六页，共八十六页，2022年，8月28日

例设计一个数字带通滤波器，通带范围为0.3πrad到0.4πrad，通带内最大衰减为3dB，0.2πrad以下和0.5πrad以上为阻带，阻带内最小衰减为18dB。采用巴特沃斯型模拟低通滤波器。解

(1)数字带通滤波器技术指标为通带上截止频率

ωu=0.4πrad

通带下截止频率

ωl=0.3πrad第五十七页，共八十六页，2022年，8月28日

阻带上截止频率

ωs2=0.5πrad

阻带下截止频率

ωs1=0.2πrad

通带内最大衰减αp=3dB，阻带内最小衰减αs=18dB。第五十八页，共八十六页，2022年，8月28日(2)模拟带通滤波器技术指标如下：设T=1，则有(通带中心频率)(带宽)第五十九页，共八十六页，2022年，8月28日

将以上边界频率对带宽B归一化，得到

ηu=3.348,ηl=2.348;ηs2=4.608,ηs1=1.498;η0=2.804(3)模拟归一化低通滤波器技术指标：归一化阻带截止频率归一化通带截止频率λp=1αp=3dB,αs=18dB第六十页，共八十六页，2022年，8月28日(4)设计模拟低通滤波器：查表，得到归一化低通传输函数G(p),第六十一页，共八十六页，2022年，8月28日(5)将归一化模拟低通转换成模拟带通：

(6)通过双线性变换法将Ha(s)转换成数字带通滤波器H(z)。下面将(5)、(6)两步合成一步计算：第六十二页，共八十六页，2022年，8月28日

将上式代入(5)中的转换公式，得将上面的p等式代入G(p)中，得第六十三页，共八十六页，2022年，8月28日

例设计一个数字带阻滤波器，通带下限频率ωl=0.19π,阻带下截止频率ωs1=0.198π，阻带上截止频率ωs2=0.202π，通带上限频率ωu=0.21π，阻带最小衰减αs=13dB，ωl和ωu处衰减αp=3dB。采用巴特沃斯型。解

(1)数字带阻滤波器技术指标：

ωl=0.19πrad,ωu=0.21πrad,αp=3dB;ωs1=0.198πrad,ωs2=0.202πrad,αs=13dB第六十四页，共八十六页，2022年，8月28日(2)模拟带阻滤波器的技术指标：设T=1，则有阻带中心频率平方为

Ω20=ΩlΩu=0.421阻带带宽为

B=Ωu-Ωl=0.07rad/s第六十五页，共八十六页，2022年，8月28日

将以上边界频率对B归一化：

ηl=8.786,ηu=9.786,ηs1=9.186,ηs2=9.386;η20=ηlηu=85.98(3)模拟归一化低通滤波器的技术指标：按照(6.2.48)式，有

λp=1,αp=3dB第六十六页，共八十六页，2022年，8月28日(4)设计模拟低通滤波器：(5)将G(p)转换成模拟阻带滤波器Ha(s)：第六十七页，共八十六页，2022年，8月28日(6)将Ha(s)通过双线性变换，得到数字阻带滤波器H(z)。第六十八页，共八十六页，2022年，8月28日6.6IIR数字滤波器的直接设计法1.零极点累试法称为零极点累试法。在确定零极点位置时要注意：

(1)极点必须位于z平面单位圆内，保证数字滤波器因果稳定；

(2)复数零极点必须共轭成对，保证系统函数有理式的系数是实的。第六十九页，共八十六页，2022年，8月28日

图6.6.1例图(a)零极点分布；(b)幅度特性第七十页，共八十六页，2022年，8月28日2.在频域利用幅度平方误差最小法直接设计IIR数字滤波器设IIR滤波器由K个二阶网络级联而成，系统函数用H(z)表示，(6.6.1)

式中，A是常数；ai,bi,ci,di是待求的系数；Hd(ejω)是希望设计的滤波器频响。如果在(0，π)区间取N点数字频率ωi,i=1,2,:,N，在这N点频率上，比较|Hd(ejω)|和|H(ejω)|，写出两者的幅度平方误差E为(6.6.2)第七十一页，共八十六页，2022年，8月28日

而在(6.6.1)式中共有(4K+1)个待定的系数，求它们的原则是使E最小。下面我们研究采用(6.6.1)式网络结构，如何求出(4K+1)系数。按照(6.6.2)式，E是(4K+1)个未知数的函数，用下式表示：

上式θ表示4K个系数组成的系数向量。为推导公式方便，令(6.6.3)第七十二页，共八十六页，2022年，8月28日

为选择A使E最小，令(6.6.4)

设θk是θ的第k个分量(ak或bk或ck或dk)，(6.6.5)

因为,式中H*i表示对Hi函数共轭。(6.6.6)第七十三页，共八十六页，2022年，8月28日

将上式具体写成对ak,bk,ck,dk的偏导，得到：(6.6.7)式中，k=1,2,3,:,K；i=1,2,3,:,N。同理求得(6.6.8)(6.6.9)(6.6.10)第七十四页，共八十六页，2022年，8月28日

由于系统函数是一个有理函数，极、零点均以共轭成对的形式存在，对于极点z1,一定有下面关系：(6.6.11)第七十五页，共八十六页，2022年，8月28日

图6.6.2例图(a)要求的幅度特性；(b)k=1，2时的幅度特性

例6.6.2设计低通数字滤波器，其幅度特性如图6.6.2(a)所示。截止频率ωs=0.1πrad。第七十六页，共八十六页，2022年，8月28日

解考虑到通带和过渡带的重要，在0~0.2π区间，每隔0.01π取一点ωi值，在0.2π~π区间每隔0.1π取一点ωi值，并增加一点过渡带，在ω=0.1π处

|Hd(ejω)|=0.5。

1.0,ω=0,0.01π,0.02π,:,0.09π0.5,ω=0.1π0.0,ω=0.11π,0.12π,:,0.19π0.0,ω=0.2π,0.3π,:,πN=29，取k=1，系统函数为第七十七页，共八十六页，2022年，8月28日

待求的参数是A，a1,b1,c1,d1。设初始值θ=(0000.25)T经过90

次迭代，求得E=1.2611，系统函数零、极点位置为零点0.67834430±j0.73474418;

极点0.75677793±j1.3213916

为使滤波器因果稳定，将极点按其倒数搬入单位圆内，再进行62次优化迭代，求得结果为

零点0.82191163±j0.56961501;

极点0.89176390±j0.19181084;Ag=0.11733978,E=0.56731第七十八页，共八十六页，2022年，8月28日

误差函数用下式表示：(6.6.12)3.在时域直接设计IIR数字滤波器设我们希望设计的IIR数字滤波器的单位脉冲响应为hd(n)，要求设计一个单位脉冲响应h(n)充分逼近hd(n)。下面我们介绍这种设计方法。设滤波器是因果性的，系统函数为(6.6.13)第七十九页，共八十六页，2022年，8月28日

式中a0=1，未知系数ai和bi共有N+M+1个，取h(n)的一段，0≤n≤p-1，使其充分逼近hd(n)，用此原则求解M+N+1个系数。将(6.6.13)式改写为令p=M+N+1,则(6.6.14)第八十页，共八十六页，2022年，8月28日

令上面等式两边z的同幂次项的系数相等，可得到N+M+1个方程：

h(0)=