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

第6章无限脉冲响应数字滤波器的设计(IIR)6.1数字滤波器的基本概念6.2模拟滤波器的设计6.3用脉冲响应不变法设计IIR数字低通滤波器6.4用双线性变换法设计IIR数字低通滤波器6.5数字高通、带通和带阻滤波器的设计

6.1数字滤波器的基本概念

一、数字滤波器的分类

数字滤波器是指输入输出均为数字信号，通过一定运算关系改变输入信号所含频率成份的相对比例或滤除某些频率成份的器件。1.从功能上分：低通、带通、高通、带阻滤波器。2.从实现方法上分：FIR、IIR滤波器3.从设计方法上来分：Chebyshev(切比雪夫）,Butterworth（巴特沃斯）等4.从处理信号分：经典滤波器、现代滤波器等等。

无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器的系统函数分别为：图5.1.1理想低通、高通、带通、带阻滤波器幅度特性二、性能指标我们在进行滤波器设计时，需要确定其性能指标。因为理想滤波器物理上是不可实现的。（由于从一个频带到另一个频带之间的突变）要物理可实现：应从一个带到另一个带之间设置一个过渡带且在通带和止带内也不应该严格为1或零。应给以较小容限。通常用的数字滤波器一般属于选频滤波器。假设数字滤波器的传输函数H(ejω)用下式表示：1、低通滤波器的性能指标fswsfpwpδ21-δ11ApAsfw|H(ejw)|或|H(f)δ1:通带的容限δ2：阻带容限通带截止频率：fp(wp)又称为通带上限频率。通带衰减：Ap阻带截止频率：fs(ws)又称阻带下限截止频率。阻带衰减：As2、高通滤波器的性能指标fswsfpwp1ApAsfw|H(ejw)|或|H(f)通带截止频率：fp(wp)又称为通带下限频率。通带衰减：Ap阻带截止频率：fp(ws)又称阻带上限截止频率。阻带衰减：As3、带通滤波器的性能指标fs1ws1fp1wp11ApAsfw|H(ejw)|或|H(f)通带截止频率：上限截止频率fp2(wp2)，下限截止频率fp1(wp1)。通带衰减：Ap阻带截止频率：上限截止频率fs2(ws2)，下限截止频率fs1(ws1)。阻带衰减：Asfp2wp2fs2ws24、带阻滤波器的性能指标fs1ws1fp1wp11ApAsfw|H(ejw)|或|H(f)通带截止频率：上限截止频率fp2(wp2)，下限截止频率fp1(wp1)。通带衰减：Ap阻带截止频率：上限截止频率fs2(ws2)，下限截止频率fs1(ws1)。阻带衰减：Asfp2wp2fs2ws2

通带内和阻带内允许的衰减一般用dB数表示，通带内允许的最大衰减用表示，阻带内允许的最小衰减用表示，如将|H(ej0)|归一化为1，(6.1.3)和(6.1.4)式则表示成：5、通常具体技术指标

因为，DF是一种具有频率选择性的离散线性系统。它是在确定信号与随机信号的数字处理中有着广泛的应用。所以，数字滤波器的设计是确定其系统函数并实现的过程。三、数字滤波器设计方法概述

滤波器设计的步骤1.根据任务，确定性能指标。2.用因果系统的线性时不变系统函数去逼近。3.用有限精度算法实现这个系统函数。（包括选择运算结构、选择合适的字长、有效数字处理方法。）4.用适当的软、硬件技术实现包括采用：通用计算机软件、数字滤波器硬件、或者二者结合。四、H(z)如何推导出（1）根据提出对滤波器的性能要求、频率特性（低、高、带通、带阻）来设计系统H(z).(2)根据时域波形提出要求来设计-->单位冲激响应h(n)或g(n)的形状。（3）有时也直接给出H(z).(但要求因果稳定).即间接设计法和直接设计法。五、确定DF的采用的结构及运算结构的好坏确定DF的采用的结构将会影响DF的精度、稳定性、经济性及运算速度等很多重要性质。1.计算复杂性一个运算结构应含有最少的乘法器和最少的延时器。乘法器最费时间，乘法器少，运算速度快；延时器最费存储单元，延时器少，存储器用的少，计算少。2.有限存储器的长度的影响与运算结构有关。即有时会希望使用一种运算结构，虽然它的乘法器和延时器并不是最少的，但它对存储器的有限字长效应是最不敏感的。IIR滤波器和FIR滤波器的设计方法是很不相同的。IIR滤波器设计方法有两类，经常用的一类设计方法是助于模拟滤波器的设计方法进行的。其设计步骤是：1、数字滤波器的技术指标转换成模拟滤波器指标2、模拟滤波器设计得到传输函数Ha(s)3、映射实现：从模拟滤波器再转换成数字滤波器H(z)本章介绍IIR滤波器的设计方法。

6.2模拟滤波器的设计

模拟滤波器的理论和设计方法已发展得相当成熟，且有若干典型的模拟滤波器供我们选择，如巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用。

图5-3各种理想滤波器的幅频特性

一、模拟低通滤波器的设计指标及逼近方法模拟低通滤波器的设计指标有,Ωp,和Ωs。其中Ωp和Ωs分别称为通带截止频率和阻带截止频率，是通带Ω(=0~Ωp)中的最大衰减系数，是阻带Ω≥Ωs的最小衰减系数，和一般用dB数表示。对于单调下降的幅度特性，可表示成：

如果Ω=0处幅度已归一化到1，即|Ha(j0)|=1

和表示为

以上技术指标用图6.2.2表示。图中Ωc称为3dB截止频率，因(6.2.3)(6.2.4)图低通滤波器的幅度特性

滤波器的技术指标给定后，需要设计一个传输函数Ha(s)，希望其幅度平方函数满足给定的指标和，一般滤波器的单位冲激响应为实数，因此说明：Ha(s)必须稳定，因此极点必须在s平面的左半平面，Ha(-s)的极点在右半平面。二、巴特沃斯低通滤波器的设计方法

1.幅度平方函数Butterworth低通滤波器具有通带最大平坦的幅度特性，是一全极点型滤波器，且极点均匀分布上Ωc的园上，并且与虚轴对称。其特点：在通带内，幅频特平坦，随着频率的升高而单调下降。其幅度平方函数为

其中N为整数，表示滤波器的阶次，Ωc定义为截止频率，为振幅响应衰减到-3dB处的频率;2、Butterworth滤波器阶数N与幅度响应的关系

图5-5巴特沃斯幅度特性和N的关系当N增大时，滤波器的特性曲线变得陡峭，则更接近理想矩形幅度特性。

将幅度平方函数|Ha(jΩ)|2写成s的函数：jΩ=s代入得：

此式表明幅度平方函数有2N个极点，极点sk用下式表示：K=0,1,2,…(2N-1)3、Butterworth滤波器的极点分布图6-6三阶巴特沃斯滤波器极点分布2N个极点等间隔分布在半径为Ωc的圆上（Butterworth圆）

为形成稳定的滤波器，2N个极点中只取s平面左半平面的N个极点构成Ha(s)，而右半平面的N个极点构成Ha(-s)。Ha(s)的表示式为设N=3，极点有6个，它们分别为取s平面左半平面的极点s0,s1,s2组成Ha(s)：

由于各滤波器的幅频特性不同，为使设计统一，将所有的频率归一化。这里采用对3dB截止频率Ωc归一化，归一化后的Ha(s)表示为

式中，s/Ωc=jΩ/Ωc。

令λ=Ω/Ωc，λ称为归一化频率；令p=jλ，p称为归一化复变量，这样归一化巴特沃斯的传输函数为(6.2.10)4、归一化的Butterworth滤波器的系统函数

式中，pk为归一化极点，用下式表示：(6-12)

将极点表示式(6-12)代入式，得到的Ha(p)的分母是p的N阶多项式，用下式表示：

(6-11)(6-11)5.确定阶数N再将幅度平方代入可得：

由(6.2.14)和(6.2.15)式得到：令,则N由下式表示：

用上式求出的N可能有小数部分，应取大于等于N的最小整数。关于3dB截止频率Ωc，如果技术指标中有给出，

总结以上，低通巴特沃斯滤波器的设计步骤如下：(1)根据技术指标Ωp,αp,Ωs和αs，求出滤波器的阶数N。(2)，求出归一化极点pk，将pk代入(6.2.11)式，得到归一化传输函数Ha(p)。

(3)将Ha(p)去归一化。将p=s/Ωc代入Ha(p)，得到实际的滤波器传输函数Ha(s)。表6-1巴特沃斯归一化低通滤波器参数例6.2.1已知通带截止频率fp=5kHz，通带最大衰减αp=2dB，阻带截止频率fs=12kHz，阻带最小衰减αs=30dB，按照以上技术指标设计巴特沃斯低通滤波器。

解(1)确定阶数N。(2)其极点为归一化传输函数为

上式分母可以展开成为五阶多项式，或者将共轭极点放在一起，形成因式分解形式。这里不如直接查表5-1简单，由N=5，直接查表得到：

极点：-0.3090±j0.9511,-0.8090±j0.5878;-1.0000

式b0=1.0000,b1=3.2361,b2=5.2361,b3=5.2361,b4=3.2361(3)为将Ha(p)去归一化，先求3dB截止频率Ωc。得到：将Ωc代入，得到：将p=s/Ωc代入Ha(p)中得到：MATLAB程序%li6.2.1fp=5kHz,Rp=2,fs=12kHz,As=30wp=2*pi*5000,ws=2*pi*12000,Rp=2,As=30;ksp=((10^(0.1*Rp)-1)/(10^(0.1*As)-1))^0.5Asp=ws/wpN=-(log(ksp)/log(Asp))N=ceil(N)wc=wp/((10^(Rp/10)-1)^(1/(2*N)))[z,p,k]=buttap(N)%求butterworth滤波器的归一化原形型，零极点增益形式%确定butturworth滤波器的系统函数N=5wc=3.3147e+004z=[]p=-0.3090+0.9511i

-0.3090-0.9511i

-0.8090+0.5878i

-0.8090-0.5878i

-1.0000

运行结果如下：%li6.2.1fp=5kHz,Rp=2,fs=12kHz,As=30wp=2*pi*5000,ws=2*pi*12000,Rp=2,As=30;[N,wc]=buttord(wp,ws,Rp,As,‘s’)

%‘s’表示模拟滤波器，缺省为数字滤波器[B,A]=butter(N,wc,'s');k=0:511;fk=0:14000/512:14000;wk=2*pi*fk;[h,w]=freqs(B,A,wk);subplot(2,2,1);plot(fk/1000,20*log10(abs(h)));title(‘buttorth滤波器的幅频特性')subplot(2,2,2);plot(fk/1000,20*log10(angle(h)));title('butterworth滤波器的相频特性')k=1N=5wc=3.7792e+004三、切贝雪夫低通滤波器Chebyshev

Butterworth滤波器频率特性，无论在通带与阻带都随频率而单调变化，因此如果在通带边缘满足指标，则在通带内肯定会有富裕量，也就是会超过指标的要求，因而并不经济，所以更有效的方法是将指标的精度要求均匀地分布在通带内，或均匀分布在阻带内，或同时均匀在通带与阻带内，这时就可设计出阶数较低的滤波器。这种精度均匀分布的办法可通过选择具有等波纹特性的逼近函数来完成。Chebyshev滤波器的分类

在一个频带中，通带或阻带具有这种等纹特性可分为：（1）ChebyshevI型：在通带中是等波纹的，在阻带内是单调的；（2）ChebyshevII型：在通带中是单调的，在阻带内是等波纹的；由应用的要求，决定采用哪种型式的Chebyshev滤波器（1）ChebyshevI型幅频特性和零极点图（N=3）N=3ChebyshevI型（2）ChebyshevII型幅频特性和零极点图（N=3）N=3ChebyshevII型，其设计思想同ChebyshevI型，在此课程中我们就不作介绍。

1ChebyshevI型幅度平方函数我们这里仅介绍切比雪夫Ⅰ型滤波器的设计方法。图分别画出阶数N为奇数与偶数时的切比雪夫Ⅰ型滤波器幅频特性。其幅度平方函数用A2(Ω)表示：图6.2.7切比雪夫Ⅰ型滤波器幅频特性ε为小于1的正数，表示通带内幅度波动的程度，ε愈大，波动幅度也愈大。Ωp称为通带截止频率。2CN(x):N阶Chebyshev多项式

(6.2.20)

图5-8示出了阶数N=0,4,5时的切比雪夫多项式特性。

由图可见：(1)切比雪夫多项式的过零点在|x|≤1的范围内；(2)当|x|<1时，|CN(x)|≤1,在|x|<1范围内具有等波纹性；(3)当|x|>1时，CN(x)是双曲线函数，随x单调上升。图6-8N=0,4,5切比雪夫多项式曲线3、通带等波纹振荡4、确定通带内波纹值ε

图切比雪夫Ⅰ型与巴特沃斯低通的A2(Ω)曲线5、确定阶数N和截止频率Ωc

阶数N等于通带内最大和最小值个数的总和。可由幅频特性中看出N阶数。且当：N=奇数，则Ω=0处有一最大值，N=偶数，则Ω=0处有一最小值。N=3和N=5N=4和N=6

设阻带的起始点频率(阻带截止频率)用Ωs表示，在Ωs处的A2(Ωs)用(6.2.19)式确定：

令λs=Ωs/Ωp，由λs>1，有可以解出ch3dB截止频率用Ωc表示，按照(6.2.19)式，有通常取λc>1，因此上式中仅取正号，得到3dB截止频率计算公式：

以上Ωp,ε和N确定后，可以求出滤波器的极点，并确定Ha(p)，p=s/Ωp。求解的过程请参考有关资料。下面仅介绍一些有用的结果。设Ha(s)的极点为si=σi+jΩi，可以证明：式中6.求滤波器的极点，确定Ha(p)上式是一个椭圆方程，长半轴为Ωpchξ(在虚轴上)，短半轴为Ωpshξ(在实轴上)。令bΩp和aΩp分别表示长半轴和短半轴，可推导出：切比雪夫滤波器的极点分布在一个椭圆上图三阶切比雪夫滤波器的极点分布

设N=3，平方幅度函数的极点分布如图6.2.8所示(极点用X表示)。为稳定，用左半平面的极点构成Ha(p)，即

式中c是待定系数。根据幅度平方函数(6.2.19)式可导出：c=ε·2N-1，代入(6.2.32)式，得到归一化的传输函数为

按照以上分析，下面介绍切比雪夫Ⅰ型滤波器设计步骤。1)确定技术要求αp,Ωp,αs和Ωsαp是Ω=Ωp时的衰减系数，αs是Ω=Ωs时的衰减系数，它们为7.求去归一化后的传输函数(6.2.33b)(6.2.34)(6.2.35)

这里αp就是前面定义的通带波纹δ，见(6.2.21)式。归一化频率2)求滤波器阶数N和参数ε

由(6.2.19)式，得到：

将以上两式代入(6.2.34)式和(6.2.35)式，得到：令(6.2.36)(6.2.37)

这样，先由(6.2.36)式求出k-11，代入(6.2.37)式，求出阶数N，最后取大于等于N的最小整数。

按照(6.2.22)式求ε，这里αp=δ。ε2=100.1δ-13)求归一化传输函数Ha(p)

为求Ha(p)，先按照(6.2.27)式求出归一化极点pk,k=1,2,:,N。

将极点pk代入(6.2.33)式，得到：4)将Ha(p)去归一化，得到实际的Ha(s)，即(6.2.38)(6.2.39)

例6.2.2设计低通切比雪夫滤波器，要求通带截止频率fp=3kHz，通带最大衰减αp=0.1dB，阻带截止频率fs=12kHz，阻带最小衰减αs=60dB。

解：(1)滤波器的技术要求：(2)求阶数N和ε：(3)求Ha(p):由(6.2.38)式求出N=5时的极点pi，代入上式，得到：(4)将Ha(p)去归一化，得到：例6.2.2%li6.2.2Rp=0.1;As=60;wp=2*pi*3000;ws=2*pi*12000;

%function[b,a]=afd_chb1(wp,ws,Rp,As);A=10^(As/20)ep=sqrt(10^(Rp/10)-1)C=wp;R=ws/wpg=sqrt(A*A-1)/epN=ceil(log10(g+sqrt(g*g-1))/log10(R+sqrt(R*R-1)))A=1000ep=0.1526R=4g=6.5522e+003N=5%LI6.2.3Rp=0.1;As=60;wp=2*pi*3000;ws=2*pi*12000;[N,wc]=cheb1ord(wp,ws,Rp,As,'s')[B,A]=cheby1(N,Rp,wp,'s');k=0:511;fk=0:12000/512:12000;wk=2*pi*fk;[h,w]=freqs(B,A,wk);plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|')N=5wc=1.8850e+004％cheb2Rp=0.1;As=60;wp=2*pi*3000;ws=2*pi*12000;[N,wc]=cheb2ord(wp,ws,Rp,As,'s')[B,A]=cheby2(N,As,wc,'s');k=0:511;fk=0:12000/512:12000;wk=2*pi*fk;[h,w]=freqs(B,A,wk);plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|')N=5wc=6.4185e+004%li6.2.4Rp=0.1;As=60;wp=2*pi*3000;ws=2*pi*12000;[N,wc]=ellipord(wp,ws,Rp,As,'s')[B,A]=ellip(N,Rp,As,wc,'s');k=0:511;fk=0:14000/512:14000;wk=2*pi*fk;[h,w]=freqs(B,A,wk);subplot(1,2,1);plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|')subplot(122);plot(fk/1000,20*log10(angle(h)));gridonxlabel('kHz');ylabel('angle')

例Rp=0.5;As=50;wp=0.3;ws=0.4;%确定butterworth滤波器的阶数和截止频率[N,Wn]=buttord(wp,ws,Rp,As)%确定butturworth滤波器的系统函数[b,a]=butter(N,Wn,'low');printsys(b,a,'s')[h,w]=freqz(b,a,501);subplot(2,2,1);plot(w/pi,abs(h));title('butterworth滤波器的幅频特性')subplot(2,2,2);plot(w/pi,angle(h));title('butterworth滤波器的相频特性')%确定Chebyshev1型滤波器的阶数和截止频率[Nche1,Wcche1]=cheb1ord(wp,ws,Rp,As)%确定Chebyshev1型滤波器的系统函数[bche1,ache1]=cheby1(Nche1,Rp,Wcche1);[hche1,wche1]=freqz(bche1,ache1,501);printsys(bche1,ache1,'s')figuresubplot(2,2,1);plot(wche1/pi,abs(hche1));title('Chebyshev1滤波器的幅频特性')subplot(2,2,2);plot(wche1/pi,angle(hche1));title('Chebyshev1滤波器的相频特性')%确定椭圆滤波器的阶数和截止频率[Nelli,Wcelli]=ellipord(wp,ws,Rp,As)%确定椭圆滤波器的系统函数[belli,aelli]=ellip(Nelli,Rp,As,Wcelli);[helli,welli]=freqz(belli,aelli,501);printsys(belli,aelli,'s')subplot(2,2,3);plot(welli/pi,abs(helli));title('椭圆滤波器的幅频特性')subplot(2,2,4);plot(welli/pi,angle(helli))title('椭圆滤波器的相频特性')N=20Wn=0.3176Nche1=9Wcche1=0.3000Nelli=6Wcelli=0.3000四、模拟滤波器的频率变换——模拟高通、带通、带阻滤波器的设计

为了防止符号混淆，先规定一些符号如下：

1低通到高通的频率变换λ和η之间的关系为

上式即是低通到高通的频率变换公式，如果已知低通G(jλ)，高通H(jη)则用下式转换：(6.2.41)(6.2.40)图6.2.9低通与高通滤波器的幅度特性

模拟高通滤波器的设计步骤如下：

(1)确定高通滤波器的技术指标：通带下限频率Ω′p，阻带上限频率Ω′s，通带最大衰减αp，阻带最小衰减αs。(2)确定相应低通滤波器的设计指标：按照(6.2.40)式

将高通滤波器的边界频率转换成低通滤波器的边界频率，各项设计指标为：①低通滤波器通带截止频率Ωp=1/Ω′p；②低通滤波器阻带截止频率Ωs=1/Ω′s；③通带最大衰减仍为αp，阻带最小衰减仍为αs。(3)设计归一化低通滤波器G(p)。(4)求模拟高通的H(s)。将G(p)按照(6.2.41)式，转换成归一化高通H(q)，为去归一化，将q=s/Ωc代入H(q)中，得

(6.2.42)

解①高通技术要求：fp=200Hz,αp=3dB;fs=100Hz,αs=15dB

归一化频率②低通技术要求：例6.2.3设计高通滤波器,fp=200Hz,fs=100Hz，幅度特性单调下降，fp处最大衰减为3dB，阻带最小衰减αs=15dB。③设计归一化低通G(p)。采用巴特沃斯滤波器，故④求模拟高通H(s)：

%li6.2.5Rp=0.1;As=40;wp=2*pi*1000;ws=2*pi*4000;[N,wc]=buttord(wp,ws,Rp,As,'s')[B,A]=butter(N,wc,'s');k=0:511;fk=0:12000/512:12000;wk=2*pi*fk;[h,w]=freqs(B,A,wk);subplot(121);plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|');axis([0,10,-80,0])Rp=0.1;As=40;wp=1;ws=4;[N,wc]=buttord(wp,ws,Rp,As,'s')[B,A]=butter(N,wc,'s');subplot(122);wph=2*pi*4000;[Bh,Ah]=lp2hp(B,A,wph);k=0:511;fk=0:12000/512:12000;wk=2*pi*fk;[h,w]=freqs(Bh,Ah,wk);plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|');axis([0,10,-80,0])%高通IIRRp=0.1;As=40;ws=2*pi*1000;wp=2*pi*4000;[N,wc]=buttord(wp,ws,Rp,As,'s')[BH,AH]=butter(N,wc,'high','s');k=0:511;fk=0:6000/512:6000;wk=2*pi*fk;[h,w]=freqs(BH,AH,wk);figureplot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|');axis([0,6,-200,0])N=5wc=1.5782e+004低通到带通的频率变换

低通与带通滤波器的幅度特性如图6.2.10所示。

归一化边界频率：图6.2.10带通与低通滤波器的幅度特性

表6.2.2η与λ的对应关系

由η与λ的对应关系，得到：由表6.2.2知λp对应ηu，代入上式中，有(6.2.43)式称为低通到带通的频率变换公式。利用该式将带通的边界频率转换成低通的边界频率。下面推导由归一化低通到带通的转换公式。由于

将(6.2.43)式代入上式，得到：将q=jη代入上式，得到：为去归一化，将q=s/B代入上式，得到：(6.2.44)(6.2.45)上式就是由归一化低通直接转换成带通的计算公式。下面总结模拟带通的设计步骤。(1)确定模拟带通滤波器的技术指标，即：带通上限频率Ωu，带通下限频率Ωl下阻带上限频率Ωs1,上阻带下限频率Ωs2

通带中心频率Ω20=ΩlΩu，通带宽度B=ΩuΩl与以上边界频率对应的归一化边界频率如下：(2)确定归一化低通技术要求：λs与-λs的绝对值可能不相等，一般取绝对值小的λs，这样保证在较大的λs处更能满足要求。

通带最大衰减仍为αp，阻带最小衰减亦为αs。(3)设计归一化低通G(p)。(4)由(6.2.45)式直接将G(p)转换成带通H(s)。

例6.2.4设计模拟带通滤波器，通带带宽B=2π×200rad/s，中心频率Ω0=2π×1000rad/s，通带内最大衰减αp=3dB，阻带Ωs1=2π×830rad/s,Ωs2=2π×1200rad/s，阻带最小衰减αs=15dB。

解

(1)模拟带通的技术要求：Ω0=2π×1000rad/s,αp=3dBΩs1=2π×830rad/s,Ωs2=2π×1200rad/s,αs=15dBB=2π×200rad/s;η0=5,ηs1=4.15,ηs2=6(2)模拟归一化低通技术要求：

取λs=1.833,αp=3dB,αs=15dB。(3)设计模拟归一化低通滤波器G(p)：

采用巴特沃斯型，有

取N=3，查表6.2.1，得(4)求模拟带通H(s)：带通滤波器设计%li626Rp=1;As=20;ws=2*pi*[2000,9000];wp=2*pi*[4000,7000];[N,wc]=buttord(wp,ws,Rp,As,'s')[BB,AA]=butter(N,wc,'s');k=0:511;fk=0:15000/512:15000;wk=2*pi*fk;[h,w]=freqs(BB,AA,wk);plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|')axis([0,15,-80,0])

3低通到带阻的频率变换

低通与带阻滤波器的幅频特性如图6.2.11所示。图6.2.11低通与带阻滤波器的幅频特性Ωl和Ωu分别是下通带截止频率和上通带截止频率，Ωs1和Ωs2分别为阻带的下限频率和上限频率，Ω0为阻带中心频率

图中，Ωl和Ωu分别是下通带截止频率和上通带截止频率，Ωs1和Ωs2分别为阻带的下限频率和上限频率，Ω0为阻带中心频率，Ω20=ΩuΩl，阻带带宽B=Ωu-Ωl，B作为归一化参考频率。相应的归一化边界频率为ηu=Ωu/B,ηl=Ωl/B,ηs1=Ωs1/B,ηs2=Ωs2/B;η20=ηuηl

表6.2.3η与λ的对应关系

根据η与λ的对应关系，可得到：

且ηu-ηl=1，λp=1，(6.2.46)式称为低通到带阻的频率变换公式。将(6.2.46)式代入p=jλ，并去归一化，可得

上式就是直接由归一化低通转换成带阻的频率变换公式。(6.2.46)(6.2.47)(6.2.48)下面总结设计带阻滤波器的步骤：(1)确定模拟带阻滤波器的技术要求，即：下通带截止频率Ωl,上通带截止频率Ωu阻带下限频率Ωs1，阻带上限频率Ωs2阻带中心频率Ω20=ΩuΩl

，阻带宽度B=Ωu-Ωl它们相应的归一化边界频率为ηl=Ωl/B,ηu=Ωu/B,ηs1=Ωs1/B;ηs2=Ωs2/B,η20=ηuηl以及通带最大衰减αp和阻带最小衰减αs。(2)确定归一化模拟低通技术要求，即：

取λs和λs的绝对值较小的λs；通带最大衰减为αp，阻带最小衰减为αs。

(3)设计归一化模拟低通G(p)。(4)按照(6.2.48)式直接将G(p)转换成带阻滤波器H(s)。

例6.2.5设计模拟带阻滤波器，其技术要求为：Ωl=2π×905rad/s,Ωs1=2π×980rad/s,Ωs2=2π×1020rad/s,Ωu=2π×1105rad/s,αp=3dB,αs=25dB。试设计巴特沃斯带阻滤波器。

解(1)模拟带阻滤波器的技术要求：Ωl=2π×905,Ωu=2π×1105;Ωs1=2π×980,Ωs2=2π×1020;Ω20=ΩlΩu=4π2×1000025B=Ωu-Ωl=2π×200;

ηl=Ωl/B=4.525,ηu=Ωu/B=5.525;ηs1=Ωs1/B=4.9,ηs2=5.1;η20=ηlηu=25(2)归一化低通的技术要求：(3)设计归一化低通滤波器G(p):αp=3dB,αs=25dB(4)带阻滤波器的H(s)为%butterworth带阻IIRRp=1;As=20;wp=2*pi*[2000,9000];ws=2*pi*[4000,7000];[N,wc]=buttord(wp,ws,Rp,As,'s')[Bs,As]=butter(N,wc,'stop','s');k=0:511;fk=0:12000/512:12000;wk=2*pi*fk;[h,w]=freqs(Bs,As,wk);subplot(121)plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|');axis([0,12,-80,0])带阻滤波器设计%ellip带阻IIRRp=1;As=20;wp=2*pi*[2000,9000];ws=2*pi*[4000,7000];[N,wc]=ellipord(wp,ws,Rp,As,'s')[Bs,As]=ellip(N,Rp,As,wc,'stop','s');k=0:511;fk=0:12000/512:12000;wk=2*pi*fk;[h,w]=freqs(Bs,As,wk);subplot(122)plot(fk/1000,20*log10(abs(h)));gridonxlabel('kHz');ylabel('|H(f)|');axis([0,12,-80,0])6.3用脉冲响应不变法设计IIR

数字低通滤波器上节我们讲到模拟滤波器设计方法，现在我们要讲如何将设计好的模拟滤波器系统函数转换成我们所需的数字滤波器系统函数。转换有两种方法：冲激响应不变法和双线性变换法。冲激响应不变法是从时域出发，要求数字滤波器的冲激响应h(n)对应于模拟滤波器ha(t)的等间隔抽样。h(n)=ha(nT)，其中T是抽样周期。因此时域逼近良好。

为了保证转换后的H(z)稳定且满足技术要求，对转换关系提出两点要求：(1)因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的。(2)数字滤波器的频率响应模仿模拟滤波器的频响，s平面的虚轴映射z平面的单位圆，相应的频率之间成线性关系。一、变换原理

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

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

式中si为Ha(s)的单阶极点。将Ha(s)进行逆拉氏变换得到ha(t)：(6.3.2)

式中u(t)是单位阶跃函数。对ha(t)进行等间隔采样，采样间隔为T，t=nT代入得到：(6.3.3)对上式进行Z变换，得到数字滤波器的系统函数H(z)：(6.3.4)二、数字滤波器与模拟滤波器的频率响应映射关系设ha(t)的采样信号用ha(t)表示，对进行拉氏变换，得到:h

式中ha(nT)是ha(t)在采样点t=nT时的幅度值，它与序列h(n)的幅度值相等，即h(n)=ha(nT)，因此得到:(6.3.5)

上式表示采样信号的拉氏变换与相应的序列的Z变换之间的映射关系可用下式表示：

(6.3.6)数字滤波器的冲激响应为对应模拟滤波器冲激响应的抽样，由抽样定理可知其频谱为模拟滤波器频谱的周期延拓。将s=jΩ代入上式，得由(6.3.5)式和(6.3.8)式得到:(6.3.7)(6.3.8)(6.3.9)模拟信号ha(t)的傅里叶变换Ha(jΩ)和其采样信号的傅里叶变换之间的关系如下：(Ωs=2π/T)

上式表明将模拟信号ha(t)的拉氏变换在s平面上沿虚轴按照周期Ωs=2π/T延拓后，再按照(6.3.6)式映射关系，映射到z平面上，就得到H(z)。(6.3.6)式可称为标准映射关系。下面进一步分析这种映射关系。设按照(6.3.6)式，得到:因此得到:(6.3.10)

那么σ=0,则r=1；σ<0,则r<1；σ>0,则r>1

另外，注意到z=esT是一个周期函数，可写成为任意整数

由于在时域抽样，导致在频域内周期延拓，数字滤波器的频率响应H(ejw)为模拟滤波器频率响应的周期延拓.

存在多对一的映射关系。图6-9z=esT,s平面与z平面之间的映射关系映射规则的要点S平面上每一条宽为的横带部分，将重叠映射到z平面的整个平面上。每一横条的左半边映射到z平面单位园内，每一横条的右半边映射到z平面单位园外。S平面的虚轴（j)轴映射到z平面单位园上，虚轴上每一段长为的线段都映射到z平面单位园上一周。数字滤波器的频响并不是简单地重现模拟滤波器的频响，而是模拟滤波器频响的周期延拓。三、脉冲响应不变法的频率混叠现象

由于在时域抽样，导致在频域内周期延拓。数字滤波器的频率响应H(ejw)为模拟滤波器频率响应的周期延拓。从模拟信号到采样信号，其拉氏变换以为周期，沿虚轴进行周期延拓。如果原模拟信号的频带不限于之间，则会在的奇数倍附近产生频率混叠。设计出的数字滤波器在附近的频率特性偏离模拟滤波器在附近的频率特性，严重时使数字滤波器不能满足给定的技术指标。希望设计的滤波器是带限滤波器，若不是带限的，如高通、带阻滤波器，则要在前面加保护滤波器，滤除高频信号，以免产生频率混叠现象。

只有模拟滤波器的频谱限带于折叠频率内时，即要满足才能避免混叠失真。而实际的滤波器并非严格限带，所以用冲激响应不变法设计的数字滤波器不可避免地会产生混叠失真。所以此法只适于设计带限滤波器。图6-10脉冲响应不变法的频率混叠现象

假设没有频率混叠现象，即满足

按照(6.3.9)式，并将关系式s=jΩ代入，ω=ΩT，代入得到：

令四、公式推导

一般Ha(s)的极点si是一个复数，且以共轭成对的形式出现，在(6.3.1)式中将一对复数共轭极点放在一起，形成一个二阶基本节。如果模拟滤波器的二阶基本节的形式为极点为

可以推导出相应的数字滤波器二阶基本节(只有实数乘法)的形式为

如果模拟滤波器二阶基本节的形式为极点为五、冲激响应不变法设计IIRDF的优缺点（1）冲激响应不变法使得数字滤波器的冲激响应完全模仿模拟滤波器的冲激响应，也就是时域逼近良好。（2）模拟频率Ω和数字频率w之间呈线性关系：w=ΩT如：一个线性相位的模拟滤波器（例贝塞尔滤波器）可以映射成一个线性相位的数字滤波器。（3）缺点：由于有频率混叠效应，所以冲激响应不变法只适用于限带的模拟滤波器,如低通、带通滤波器。六、冲激不变法应用的局限性

由于具有频率的混叠效应，所以高通和带阻滤波器不宜采用冲激不变法。因为它们高频部分不衰减，将完全混淆在低频中，从而使整个频响面目全非。若要对高通和带阻实行冲激不变法，则必须先对高通和带阻滤波器加一保护滤波器，滤掉高于折叠频率以上的频带。它会增加设计的复杂性和滤波器的阶数，因而只有在一定要追求频率线性关系或保持网络瞬态响应不变时才使用。对于带通和低通滤波器，需充分限带，若阻带衰减越大，则混叠效应越小。

例6.3.1已知模拟滤波器的传输函数Ha(s)为

用脉冲响应不变法将Ha(s)转换成数字滤波器的系统函数H(z)。

解首先将Ha(s)写成部分分式：极点为b=[10.64490.7079];roots(b)ans=-0.3225+0.7771i-0.3225-0.7771i用MATLAB求极点：按照(6.3.4)式，并经过整理，得到

设T=1s时用H1(z)表示，T=0.1s时用H2(z)表示，则

转换时，也可以直接按照(6.3.13),(6.3.14)式进行转换。首先将Ha(s)写成(6.3.13)式的形式，如极点s1,2=σ1±jΩ1,则那么H(z)的极点为再按照(6.3.14)式，H(z)为图6.3.3例6.3.1的幅度特性例2由于模拟滤波器不是充分限带，所以数字滤波器产生很大的频谱混叠失。|Ha(jΩ)|Ωw|H(ejw)|%li6.3.2T=1;Rp=1;As=10;wp=0.2*pi/T;ws=0.35*pi/T;[N,wc]=buttord(wp,ws,Rp,As,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=impinvar(B,A);[h,w]=freqs(B,A);[H,w1]=freqz(Bz,Az);printsys(b,a,'s')printsys(Bz,Az,'z')subplot(221)plot(w,20*log10(abs(h)));gridonxlabel('f/Hz');ylabel('|H(f)dB');axis([0,5,-50,0])subplot(222)plot(w1,20*log10(abs(H)));gridonxlabel('w/pi');ylabel('|H(f)dB|');axis([0,5,-50,0])%T=0.1sT=0.1;fs=10;Rp=1;As=10;wp=0.2*pi/T;ws=0.35*pi/T;[N,wc]=buttord(wp,ws,Rp,As,'s');[B,A]=butter(N,wc,'s');printsys(b,a,'s')[Bz,Az]=impinvar(B,A,fs);printsys(Bz,Az,'z')[h,w]=freqs(B,A);[H,w1]=freqz(Bz,Az);subplot(223);plot(w,20*log10(abs(h)));gridonxlabel('f/Hz');ylabel('|H(f)dB');axis([0,50,-50,0])subplot(224);plot(w1,20*log10(abs(H)));gridonxlabel('w/pi');ylabel('|H(f)dB|');axis([0,5,-50,0])6.4用双线性变换法设计IIR数字

低通滤波器

冲激不变法：是使数字滤波器在时域上模仿模拟滤波器，但它的缺点：产生频率响应的混叠失真。这是由于从S平面->Z平面是多值的映射关系所造成的。为了克服这一缺点，我们采用双线性变换法.1.定义双线性变换法：是从频域出发，使DF的频率响应与AF的频率响应相似的一种变换法。采用频率压缩方法，将整个频率轴上的频率范围压缩到之间，再用转换到z平面上。一、变换原理2、双线性变换法的映射关系实现S平面与Z平面一一对应的关系。第一次变换：频率压缩第二次变换：数字化S平面S1平面Z平面3、双线性变换法的映射规则

映射规则（1）频率压缩：把整个S平面压缩变换到某一中介的S1平面的一条横带里。（2）数字化：将S1平面通过标准变换关系变换到z平面。（1）频率压缩把整个S平面压缩变换到某一中介的S1平面的一条横带里。（2）数字化将S1平面通过标准变换关系变换到z平面。（3）变换常数C的选择调节C，可使AF与DF在不同频率点处有对应的关系。（a）使AF与DF在低频处有较确切的对应关系。看出在低频处，AF的低频特性近似等于DF的低频特性。二、性能分析1.解决了冲激不变法的混叠失真问题。2.它是一种简单的代数关系。只须将上述关系代入AF的Ha(s)中（对直接、级联、并联结构都适用）即可求出DF的H(z),设计十分方便。3.由于双线性变换中，即模拟角频率与数字角频率存在非线性关系。所以双线性变换避免了混叠失真，却又带来了非线性的频率失真。4.双线性变换法不适用于如下设计：（1）设计线性相位的DF（2）它要求AF的幅频响应是分段常数型.(即幅度变换是线性的）。（一般低通，高通，带通，带阻型滤波器的频率响应特性都是分段常数）

5.同时，看出双线性变换：（1）在零频附近，模拟角频率与数字角频率变换关系接近线性关系。（2）又要求AF的幅频响应是分段常数型，即幅度变换是线性的所以称之为双线性变换。频率升高时，非线性失真严重。6.对于分段常数型AF滤波器，经双线性变换后，仍得到幅频特性为分段常数的DF.但在各个分段边缘的临界频率点产生畸变，这种频率的畸变，可通过频率预畸变加以校正。例1

一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位。如一个模拟微分器将不能通过双线性变换成为数字微分器。模拟微分器数字图6.4.3双线性变换法幅度和相位特性的非线性映射例2

对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘临界频率点产生了畸变。这种频率的畸变，可以通过频率的预畸变加以校正，也就是临界频率事先加以畸变，然后经变换后正好映射到所需要的频率。三、设计流程1.根据要求，确定数字低通滤波器的技术指标：通带截止频率ωp、通带衰减αp、阻带截止频率ωs、阻带衰减αs。2.将各分段频率临界点预畸变。3.将数字滤波器的性能指标转换为中间模拟滤波器的性能指标。4.根据设计要求，选定双线性变换常数C。5.设计中间模拟滤波器的系统函数Ha(s).6.将代入Ha(s)中，得到数字低通滤波器DF的系统函数的H(z).

为简化设计，已将模拟滤波器各系数和经双线性变换法得到的数字滤波器的各系数之间关系，列成表格供设计时使用。

表6-2系数关系表

例3试分别用脉冲响应不变法和双线性不变法将图6.4.4所示的RC低通滤波器转换成数字滤波器。

解首先按照图6.4.4写出该滤波器的传输函数Ha(s)为

利用脉冲响应不变法转换，数字滤波器的系统函数H1(z)为

利用双线性变换法转换，数字滤波器的系统函数H2(z)为H1(z)和H2(z)的网络结构分别如图6.4.5(a),(b)所示。图6.4.5例6.4.1图——H1(z)和H2(z)的网络结构(a)H1(z);(b)H2(z)图6.4.6例6.4.1图——数字滤波器H1(z)和H2(z)的幅频特性例6.4.2设计低通数字滤波器，要求在通带内频率低于0.2πrad时，容许幅度误差在1dB以内；在频率0.3π到π之间的阻带衰减大于15dB。指定模拟滤波器采用巴特沃斯低通滤波器。试分别用脉冲响应不变法和双线性变换法设计滤波器。

解(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③设计巴特沃斯低通滤波器。先计算阶数N及3dB截止频率Ωc。

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

根据阶数N=6，查表6.2.1，得到归一化传输函数为

为去归一化，将p=s/Ωc代入Ha(p)中，得到实际的传输函数Ha(s),④用脉冲响应不变法将Ha(s)转换成H(z)。首先将Ha(s)进行部分分式，并按照(6.3.11)式、(6.3.12)式，或者(6.3.13)式和(6.3.14)式，得到：图6.4.7例6.4.2图——用脉冲响应不变法设计的数字低通滤波器的幅度特性(2)用双线性变换法设计数字低通滤波器。①数字低通技术指标仍为ωp=0.2πrad,αp=1dB;ωs=0.3πrad,αs=15dB②模拟低通的技术指标为③设计巴特沃斯低通滤波器。阶数N计算如下：

取N=6。为求Ωc,将Ωs和αs代入(6.2.18)式中，得到Ωc=0.7662rad/s。这样阻带技术指标满足要求，通带指标已经超过。

根据N=6，查表6.2.1得到的归一化传输函数Ha(p)与脉冲响应不变法得到的相同。为去归一化，将p=s/Ωc代入Ha(p)，得实际的Ha(s)，④用双线性变换法将Ha(s)转换成数字滤波器H(z)：图6.4.8例6.4.2图——用双线性变换法设计的数字低通滤波器的幅度特性Rp=1;As=15;wp=0.2*pi;ws=0.3*pi;T=1;Fs=1;Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2);ksp=((10^(0.1*Rp)-1)/(10^(0.1*As)-1))^0.5Asp=Omegas/OmegapN=-(log(ksp)/log(Asp))N=ceil(N)%不小于自变量的最小整数wc=Omegas/((10^(As/10)-1)^(1/(2*N)))[z,p,k]=buttap(N)%求零极点[b,a]=zp2tf(z,p,k)%归一化系统函数[bt,at]=lp2lp(b,a,wp)%去归一化[bb,ab]=bilinear(bt,at,Fs)%双线性变换法转换为数字滤波器ksp=0.0920Asp=1.5682N=5.3044N=6wc=0.7662z=[]p=-0.2588+0.9659i-0.2588-0.9659i-0.7071+0.7071i-0.7071-0.7071i-0.9659+0.2588i-0.9659-0.2588ik=1运行结果如下：b=0000001a=1.00003.86377.46419.14167.46413.86371.0000bt=0.0615at=1.00002.42762.94672.26761.16330.37840.0615bb=0.00030.00170.00430.00580.00430.00170.0003ab=1.0000-3.65435.8693-5.21922.6894-0.75740.0908%li6.4.2设计IIRRp=1;As=15;wp=0.2*pi;ws=0.3*pi;T=1;Fs=1;Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2);[N,wc]=buttord(wp,ws,Rp,As,'s')[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs);[h,w]=freqs(B,A);subplot(121)plot(w,20*log10(abs(h)));gridonxlabel('Hz');ylabel('|H(f)|')axis([0,1,-20,0])[H,w1]=freqz(Bz,Az)subplot(122)plot(w1,20*log10(abs(H)));gridonxlabel('w/pi');ylabel('|H(f)dB|')axis([0,1,-20,0])

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

例如高通数字滤波器等。具体设计步骤如下：(1)确定所需类型数字滤波器的技术指标。(2)将所需类型数字滤波器的技术指标转换成所需类型模拟滤波器的技术指标，转换公式为

(3)将所需类型模拟滤波器技术指标转换成模拟低通滤波器技术指标(具体转换公式参考本章6.2节)。(4)设计模拟低通滤波器。(5)将模拟低通通过频率变换，转换成所需类型的模拟滤波器。(6)采用双线性变换法，将所需类型的模拟滤波器转换成所需类型的数字滤波器。

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

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

解(1)数字高通的技术指标为ωp=0.8πrad,αp=3dB;ωs=0.44πrad,αs=15dB(2)模拟高通的技术指标计算如下：

令T=1，则有(3)模拟低通滤波器的技术指标计算如下：15

将Ωp和Ωs对3dB截止频率Ωc归一化，这里Ωc=Ωp,(4)设计归一化模拟低通滤波器G(p)。模拟低通滤波器的阶数N计算如下：

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

为去归一化，将p=s/Ωc代入上式得到：(5)将模拟低通转换成模拟高通。将上式中G(s)的变量换成1/s，得到模拟高通Ha(s)：(6)用双线性变换法将模拟高通H

(s)转换成数字高通H(z)：实际上(5)、(6)两步可合并成一步，即%li6.5.1设计IIRRp=3;As=15;wp=0.8;ws=0.44;[N,wc]=buttord(wp,ws,Rp,As)[B,A]=b