《数字信号处理教学课件》iir滤波器.ppt

《《数字信号处理教学课件》iir滤波器.ppt》由会员分享，可在线阅读，更多相关《《数字信号处理教学课件》iir滤波器.ppt（45页珍藏版）》请在三一办公上搜索。

1、IIR滤波器的设计,一、模拟低通滤波器的设计,1、巴特沃思低通滤波器的设计,MATLAB提供了一个叫做z,p,k=buttap(N)的函数，来设计N阶归一化的巴特沃思模拟原形滤波器，它返回数组z和p（零点和极点）以及增益k。但是，我们需要具有任意c的未归一化的巴特沃思滤波器。,可以用自定义的u_buttap(N,Omegac)函数来设计未归一化的巴特沃思模拟原形滤波器。,function b,a=u_buttap(N,Omegac);%b,a=u_buttap(N,Omegac);%b=numerator polynomial coefficients of Ha(s)%a=denominat

2、or polynomial coefficients of Ha(s)%N=Order of the Butterworth Filter%Omegac=Cutoff frequency in radians/secz,p,k=buttap(N);%z,p,k为归一化的零、极点和增益p=p*Omegac;%极点去归一化k=k*OmegacN;%增益乘 B=real(poly(z);%由零点恢复多项式B(z)b0=k;b=k*B;a=real(poly(p);%由极点恢复多项式A(z),模拟低通滤波器是由参数p，Rp，s和As描述的，设计巴特沃思滤波器的基本思想是根据这些指标得到阶数N和截止频率

3、c。可以编写函数afd_butt来实现：,function b,a=afd_butt(Wp,Ws,Rp,As);if Wp 0)endif Ws Passband edge)endif(Rp 0)end,N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);fprintf(n*Butterworth Filter Order=%2.0f n,N)OmegaC=Wp/(10(Rp/10)-1)(1/(2*N);b,a=u_buttap(N,OmegaC);,例1：设计一个低通巴特沃思滤波器，以满足：通带截止频率：p=0.2，通带波动：Rp

4、=7dB 阻带截止频率：s=0.3，阻带波动：As=16dB,wp=0.2*pi;Ws=0.3*pi;Rp=7;As=16;%模拟滤波器设计:b,a=afd_butt(Wp,Ws,Rp,As)%计算频率响应:w=0:1:500*pi/500;h=freqs(b,a,w);%画图subplot(2,2,1);plot(w/pi,abs(h);title(幅度响应);grid;subplot(2,2,3);plot(w/pi,20*log10(abs(h);title(幅度响应(dB);grid;subplot(2,2,2);plot(w/pi,angle(h)*180/pi);title(相位响

5、应);grid;subplot(2,2,4);impulse(b,a);title(脉冲响应);grid;,2、切比雪夫型低通滤波器的设计,MATLAB提供了一个函数z,p,k=cheblap(N,Rp)，来设计一个阶数为N，通带波动为Rp的归一化切比雪夫型原型滤波器。,可以自编写函数u_chblap用于设计一个未归一化的切比雪夫型原型滤波器，它返回直接形式的Ha(s)。,以函数u_chblap为基础，可以编写一个afd_chb1函数，它根据给定的技术指标设计模拟切比雪夫型滤波器。,例2：设计一个低通切比雪夫型滤波器，以满足：通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0

6、.3，阻带波动：As=16dB,Wp=0.2*pi;Ws=0.3*pi;Rp=1;As=16;%模拟滤波器设计:b,a=afd_chb1(Wp,Ws,Rp,As);%计算频率响应:w=0:1:500*pi/500;h=freqs(b,a,w);%画图subplot(2,2,1);plot(w/pi,abs(h);title(幅度响应);grid;subplot(2,2,3);plot(w/pi,20*log10(abs(h);title(幅度响应(dB);grid;subplot(2,2,2);plot(w/pi,angle(h)*180/pi);title(相位响应);grid;subplo

7、t(2,2,4);impulse(b,a);title(脉冲响应);grid;,3、切比雪夫型低通滤波器的设计,MATLAB提供了一个函数z,p,k=cheb2ap(N,As)，来设计一个阶数为N，通带波动为As的归一化切比雪夫型原型滤波器。,函数u_chb2ap用于设计一个未归一化的切比雪夫型原型滤波器，它返回直接形式的Ha(s)。,以函数u_chb2ap为基础，可以编写了一个afd_chb2函数，它根据给定的技术指标设计模拟切比雪夫型滤波器。,例3：设计一个低通切比雪夫型滤波器，以满足：通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=16dB,Wp

8、=0.2*pi;Ws=0.3*pi;Rp=1;As=16;%模拟滤波器设计:b,a=afd_chb2(Wp,Ws,Rp,As);%计算频率响应:w=0:1:500*pi/500;h=freqs(b,a,w);%画图subplot(2,2,1);plot(w/pi,abs(h);title(幅度响应);grid;subplot(2,2,3);plot(w/pi,20*log10(abs(h);title(幅度响应(dB);grid;subplot(2,2,2);plot(w/pi,angle(h)/pi);title(相位响应);grid;subplot(2,2,4);impulse(b,a);

9、title(脉冲响应);grid;,4、椭圆低通滤波器的设计,MATLAB提供了一个函数z,p,k=ellipap(N,Rp,As)，来设计一个阶数为N，通带波动为Rp，阻带衰减为As的归一化椭圆原型滤波器。,函数u_elipap用于设计一个未归一化的椭圆原型滤波器，它返回直接形式的Ha(s)。,以函数u_elipap为基础，我们提供了一个afd_elip函数，它根据给定的技术指标设计模拟椭圆低通滤波器。,例4：设计一个低通椭圆滤波器，以满足：通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=16dB,wp=0.2*pi;Ws=0.3*pi;Rp=1;A

10、s=16;%模拟滤波器设计:b,a=afd_elip(Wp,Ws,Rp,As);%计算频率响应:w=0:1:500*pi/500;h=freqs(b,a,w);%画图subplot(2,2,1);plot(w/pi,abs(h);title(幅度响应);grid;subplot(2,2,3);plot(w/pi,20*log10(abs(h);title(幅度响应(dB);grid;subplot(2,2,2);plot(w/pi,angle(h)/pi);title(相位响应);grid;subplot(2,2,4);impulse(b,a);title(脉冲响应);grid;,二、模拟数字

11、滤波器变换,1、脉冲响应不变变换,脉冲响应不变变换法的步骤如下：,1)确定T并选择模拟频率：p=p/T,s=s/T 2)根据指标p，Rp，s和As，设计模拟滤波器Ha(s)，这 可由上面所讲的三种原形之一来实现。3)利用部分分式展开式，把Ha(s)展成：4)把模拟极点Pk转换成数字极点，得到数字滤波器,MATLAB工具箱提供了函数impinvar，它采用冲激响应不变法来实现模拟滤波器到数字滤波器的转换。,bz，az=impinvar(b，a，Fs)可将模拟滤波器(b，a)变换成数字滤波器(bz，az)，两者的冲激响应不变，即模拟滤波器的冲激响应按h抽样后等同于数字滤波器的冲激响应。当缺少参数F

12、s时，抽样频率Fs取缺省值1Hz。,例5：用巴特沃思滤波器原型设计一个低通数字滤波器，满足通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=15dB,wp=0.2*pi;%digital Passband freq in Hzws=0.3*pi;%digital Stopband freq in HzRp=1;%Passband ripple in dBAs=15;%Stopband attenuation in dB%模拟原型指标的频率逆映射T=1;Fs=1/T;%Set T=1OmegaP=wp/T;%Prototype Passband freqO

13、megaS=ws/T;%Prototype Stopband freq%模拟巴特沃思原型滤波器计算cs,ds=afd_butt(OmegaP,OmegaS,Rp,As);%脉冲响应不变变换b,a=impinvar(cs,ds,Fs);,%画图h,w=freqz(b,a);subplot(2,2,1);plot(w/pi,abs(h);title(幅度响应);grid;subplot(2,2,2);plot(w/pi,angle(h);title(相位响应);grid;subplot(2,2,3);plot(w/pi,20*log10(abs(h);title(幅度响应(dB);grid;n=0

14、:1:59;imp=1;zeros(59,1);y=filter(b,a,imp);subplot(2,2,4);plot(n,y);title(脉冲响应);grid;,2、双线性变换,MATLAB工具箱提供了函数bilinear()实现双线性变换。,numd，dend=bilinear(num，den，Fs)把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型，其中Fs为抽样频率。双线性变换步骤如下:,1)选择T值，设计时它是任意的，可以置为1。2)根据截止频率p和s进行预修正，即：3)设计模拟滤波器。4)最后，设 化简得到为z-1 的有理函数的H(z)。,例6：用切比雪夫I型滤波器原型

15、设计一个低通数字滤波器,满足通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=15dB,%数字滤波器指标wp=0.2*pi;%digital Passband freq in Hzws=0.3*pi;%digital Stopband freq in HzRp=1;%Passband ripple in dBAs=15;%Stopband attenuation in dB%模拟原型指标的频率逆映射T=1;Fs=1/T;%Set T=1OmegaP=(2/T)*tan(wp/2);%预畸 OmegaS=(2/T)*tan(ws/2);%预畸%模拟切比雪夫

16、I型原型滤波器的计算cs,ds=afd_chb1(OmegaP,OmegaS,Rp,As);%双线性变换b,a=bilinear(cs,ds,Fs);,%画图h,w=freqz(b,a);subplot(2,2,1);plot(w/pi,abs(h);title(幅度响应);grid;subplot(2,2,2);plot(w/pi,angle(h);title(相位响应);grid;subplot(2,2,3);plot(w/pi,20*log10(abs(h);title(幅度响应(dB);axis(0,1,-40,5);grid;n=0:1:59;imp=1;zeros(59,1);y=

17、filter(b,a,imp);subplot(2,2,4);plot(n,y);title(脉冲响应);grid;,三、选频IIR数字滤波器设计方法,（一）典型IIR选频滤波器设计方法,(1)按一定规则将给出的数字滤波器的技术指标转换成模拟低 通滤波器的技术指标。(2)根据转换后的技术指标使用滤波器阶数选择函数，确定最 小阶数N和固定频率Wn。(3)运用最小阶数N产生模拟滤波器原型。(4)运用固有频率Wn把模拟低通滤波器原型转换成模拟低通、高通、带通、带阻滤波器。(5)运用冲激响应不变法或双线性变换法把模拟滤波器转换成 数字滤波器。,模拟低通原型滤波器的创建buutap cheb1ap ch

18、ab2apellipap,频率转换lp2lp lp2hplp2bp lp2bs,离散化处理bilinearimpinvar,1、模拟低通原型滤波器的设计,MATLAB提供了buutap、cheb1ap、chab2ap和ellipap函数设计不同类型的模拟低通原型滤波器。,它们的函数形式如下：z.p,k=buttap(n)z.p,k=cheb1ap(n,rp)z.p,k=cheb2ap(n,rs)z.p,k=ellipap(n,rp,rs),buutap、cheb1ap、chab2ap和ellipap函数以零、极点增益的形式给出模拟原型低通滤波器的设计方案。其中，z、p、k分别表示零点、极点和增

19、益。,MATLAB提供函数zp2tf将零极点增益模型转换为传递函数模型。函数形式为：,b,a=zp2tf(z,p,k),由于通常给出的数字滤波器的设计指标是wp、ws、rp和rs，而不是n，所以在利用函数buutap、cheb1ap、chab2ap和ellipap函数设计模拟低通原型滤波器之前，应求出滤波器的阶数n。,利用函数buttord、cheb1ord、cheb2ord和ellipord函数可以求出模拟滤波器的阶数。n,wn=buttord(wp,ws,rp,rs,s)n,wn=cheb1ord(wp,ws,rp,rs,s)n,wn=cheb2ord(wp,ws,rp,rs,s)n,wn

20、=ellipord(wp,ws,rp,rs,s),注意：（1）buttord、cheb1ord、cheb2ord和ellipord函数中的s表示所求的滤波器为模拟滤波器。因此，指标wp和ws实际上是模拟角频率而不是数字角频率，其单位是rad/s，而不是rad。（2）对于低通滤波器：wpws。对于带通滤波器：wp和ws均为二维向量，且wp(1)ws(1)、wp(2)ws(2)。,2、频率转换,MATLAB提供函数lp2lp、lp2hp、lp2bp和lp2bs实现原型低通滤波器向低通、高通、带通和带阻滤波器的转换,（1）从低通到低通的转换,lp2lp函数可将截止频率为1rad/s的模拟低通滤波器原

21、型变换成截止频率为Wn的低通滤波器。bt，at=lp21p(b，a，Wn)其中Wn为截止频率。,（2）从低通到高通的转换,lp2hp函数可将截止频率为1rad/s的模拟低通滤波器原型变换成截止频率为Wn的高通滤波器。bt，at=lp2hP(b，a，Wn),（3）从低通到带通的转换,lp2bp函数可将截止频率为1rad/s的模拟低通滤波器原型转换成具有指定带宽Bw和中心频率Wn的带通滤波器。bt，at=1p2bp(b，a，Wn，Bw)其中心频率为Wn，带宽为Bw。如果要求的滤波器的低端截止频率为w1，高端截止频率为w2，则可计算出Wn和Bw：Wn=sqrt(w1,w2)Bw=w2-w1,（4）从

22、低通到带阻的转换,lp2bs函数可将截止频率为1rad/s的模拟低通滤波器原型转换成具有指定带宽Bw和中心频率Wn的带阻滤波器。bt，at=lp2bs(b，a，Wn，Bw)其中心频率为Wn，带宽为Bw。Wn和Bw的计算方法与（3）相同。如果要求的滤波器的低端截止频率为w1，高端截止频率为w2，则可计算出Wn和Bw：Wn=sqrt(w1,w2)Bw=w2-w1,3、滤波器离散化,（1）冲激响应不变法 MATLAB工具箱提供了函数impinvar，它采用冲激响应不变法来实现模拟滤波器到数字滤波器的转换。,（2）双线性变换法 为了克服冲激响应不变法产生的频率混叠现象，需要使s平面与z平面建立一一对应

23、的单值映射关系，可采用双线性变换法。MATLAB工具箱提供了函数bilinear()，实现双线性变换。,例7：设计一个巴特沃斯高通数字滤波器，它的抽样频率为Fs为100Hz,通带临界频率为40Hz,rp=0.5dB，rs=40dB，通带上限临界频率为40Hz，阻带下限临界频率为30Hz。,wp=2*pi*40;ws=2*pi*30;rp=0.5;rs=40;%滤波器指标 Fs=100;T=1/Fs;wp=(2/T)*tan(wp/200)%预畸 ws=(2/T)*tan(ws/200)%预畸n,wn=buttord(wp,ws,rp,rs,s)%模拟滤波器设计z,p,k=buttap(n);b

24、1,a1=zp2tf(z,p,k);b2,a2=lp2hp(b1,a1,wn);%频带转换b,a=bilinear(b2,a2,Fs)%滤波器离散化h,w=freqz(b,a)%画滤波器的频响plot(w*Fs/(2*pi),abs(h);grid;,例8：试用双线性变换法设计一个带通椭圆数字滤波器,其技术指标为wpl=l0Hz，wsl=9Hz，wp2=20Hz，ws2=21Hz，rp=0.5B，rs=50dB，抽样频率为100Hz。,wp1=2*pi*10;wp2=2*pi*20;ws1=2*pi*9;ws2=2*pi*21;Fs=100;rp=0.5;rs=50;wp1=(2/T)*tan

25、(wp1*T/2);wp2=(2/T)*tan(wp2*T/2)ws1=(2/T)*tan(ws1*T/2);ws2=(2/T)*tan(ws2*T/2)wp=wp1 wp2;ws=ws1 ws2;n,wn=ellipord(wp,ws,rp,rs,s)Bw=wn(2)-wn(1);Wo=sqrt(wn(2)*wn(1);z,p,k=ellipap(n,rp,rs);b1,a1=zp2tf(z,p,k)b2,a2=lp2bp(b1,a1,Wo,Bw)b,a=bilinear(b2,a2,Fs)H,W=freqz(b,a);plot(W*Fs/(2*pi),abs(H);grid;axis(0

26、50 0 1);,（二）利用MATLAB函数直接设计数字滤波器,1、巴特沃思滤波器,MATLAB提供了函数butter完成巴特沃思滤波器的设计。形式如下：b,a=butter(N,wn):设计N阶低通滤波器，wn为它的3dB截止频率。b,a=butter(N,wn,high):设计N阶高通滤波器，wn为它的3dB截止频率。b,a=butter(N,wn):当wn为具有两个元素的矢量wn=w1 w2时，它设计2N 阶带通滤波器，它具有单位为的3dB通带：w1ww2。b,a=butter(N,wn,stop):当wn为具有两个元素的矢量wn=w1 w2时，它设计2N阶带阻滤波器，它具有单位为的3d

27、B阻带：w1ww2。注：wn单位为。,为了设计任意的选频巴特沃思滤波器，必须知道阶数N和3dB截止频率矢量wn。在MATLAB的信号处理工具箱中，提供了一个叫做buttord的函数计算这些参数。给定指标：wp，ws，Rp和As，这个函数可以求出必要的参数：它的句法为：,N,wn=buttord(wp,ws,Rp,As)l 对于低通滤波器 wpwsl 对于带通滤波器 wp和ws是具有两个元素的矢量，wp=wp1,wp2和ws=ws1,ws2，并且ws1wp1wp2ws2l 对于带阻滤波器 wp和ws是具有两个元素的矢量，wp=wp1,wp2和ws=ws1,ws2，并且wp1ws1ws2wp2,2

28、、切比雪夫滤波器,b,a=cheby1(N,Rp,Wn);其中Rp为指定通带内波纹。b,a=cheby1(N,Rp,Wn,ftype);其中ftype=high或stop。N,wn=cheb1ord(wp,ws,Rp,As);b,a=cheby2(N,As,Wn);其中As为指定阻带内波纹。b,a=cheby2(N,As,Wn,ftype);其中ftype=high或stop。N,wn=cheb2ord(wp,ws,Rp,As);注：wn单位为。,请在MATLAB中参考以下函数的帮助信息：cheby1、cheb1ord、cheby2、cheb2ord,3、椭圆滤波器,b,a=ellip(N,R

29、p,As,Wn);其中Rp、As分别为指定通带内和阻带内的波纹。b,a=ellip(N,Rp,As,Wn,ftype);其中ftype=high或stop。N,wn=ellipord(wp,ws,Rp,As);注：wn单位为。,请在MATLAB中参考以下函数的帮助信息：ellip、ellipord,例7：设计一个巴特沃斯高通数字滤波器，它的抽样频率为Fs为100Hz,通带临界频率为40Hz,rp=0.5dB，rs=40dB，通带上限临界频率为40Hz，阻带下限临界频率为30Hz。,wp=0.8;ws=0.6;Fs=100;n,wn=buttord(wp,ws,rp,rs)b,a=butter(n,wn,high)h,w=freqz(b,a);plot(w*Fs/(2*pi),abs(h);grid;,例8：试用双线性变换法设计一个带通椭圆数字滤波器,其技术指标为wpl=l0Hz，wsl=9Hz，wp2=20Hz，ws2=21Hz，rp=0.5B，rs=50dB，抽样频率为100Hz。,w1=0.2 0.4;w2=0.18 0.42;n,wn=ellipord(w1,w2,rp,rs)b,a=ellip(n,rp,rs,wn)H,W=freqz(b,a);plot(W*Fs/(2*pi),abs(H);grid;,