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

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

2、= denominator 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

3、，s和As描述的，设计巴特沃思滤波器的基本思想是根据这些指标得到阶数N和截止频率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_but

4、tap(N,OmegaC);,例1：设计一个低通巴特沃思滤波器，以满足： 通带截止频率：p=0.2，通带波动：Rp=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

5、(abs(h);title(幅度响应(dB);grid;subplot(2,2,2); plot(w/pi,angle(h)*180/pi); title(相位响应);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函

6、数，它根据给定的技术指标设计模拟切比雪夫型滤波器。,例2：设计一个低通切比雪夫型滤波器，以满足：通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.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/

7、pi,20*log10(abs(h);title(幅度响应(dB);grid;subplot(2,2,2); plot(w/pi,angle(h)*180/pi); title(相位响应);grid;subplot(2,2,4); impulse(b,a);title(脉冲响应);grid;,3、切比雪夫型低通滤波器的设计,MATLAB提供了一个函数z,p,k=cheb2ap(N,As)，来设计一个阶数为N，通带波动为As的归一化切比雪夫型原型滤波器。,函数u_chb2ap用于设计一个未归一化的切比雪夫型原型滤波器，它返回直接形式的Ha(s)。,以函数u_chb2ap为基础，可以编写了一个af

8、d_chb2函数，它根据给定的技术指标设计模拟切比雪夫型滤波器。,例3：设计一个低通切比雪夫型滤波器，以满足：通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=16dB,Wp = 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);

9、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;,4、椭圆低通滤波器的设计,MATLAB提供了一个函数z,p,k=ellipap(N,Rp,As)，来设计一个阶数为N，通带波动为Rp，阻带衰减为As的归一化椭圆原型滤波器。,函数u_elipap用于设计一个未归一化的椭圆原型滤波器，它返回直接形式的Ha(s)。,以函数u_elipap为基础，我们提供

10、了一个afd_elip函数，它根据给定的技术指标设计模拟椭圆低通滤波器。,例4：设计一个低通椭圆滤波器，以满足： 通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=16dB,wp = 0.2*pi; Ws = 0.3*pi; Rp = 1; As = 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)

11、; 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;,二、模拟数字滤波器变换,1、脉冲响应不变变换,脉冲响应不变变换法的步骤如下：,1) 确定T并选择模拟频率： p=p/T , s=s/T 2) 根据指标p，Rp，s和As，设计模拟滤波器Ha(s)，这 可由上面所讲的三种原形之一来实现。3) 利用部分分式展开式，把Ha(s)展成：4) 把模拟极点P

12、k转换成数字极点 ，得到数字滤波器,MATLAB工具箱提供了函数impinvar，它采用冲激响应不变法来实现模拟滤波器到数字滤波器的转换。,bz，az=impinvar(b，a，Fs)可将模拟滤波器(b，a)变换成数字滤波器(bz，az)，两者的冲激响应不变，即模拟滤波器的冲激响应按h抽样后等同于数字滤波器的冲激响应。当缺少参数Fs时，抽样频率Fs取缺省值1Hz。,例5：用巴特沃思滤波器原型设计一个低通数字滤波器，满足通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=15dB,wp = 0.2*pi; % digital Passband freq in

13、 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 freqOmegaS = ws / T; % Prototype Stopband freq% 模拟巴特沃思原型滤波器计算cs,ds = afd_butt(OmegaP,OmegaS,Rp,As);% 脉冲响应不变变换b,

14、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:1:59;imp=1;zeros(59,1);y=filter(b,a,imp);subplot(2,2,4); plot(n,y);title(脉冲响应);grid;,

15、2、双线性变换,MATLAB工具箱提供了函数bilinear()实现双线性变换。,numd，dend=bilinear(num，den，Fs)把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型，其中Fs为抽样频率。双线性变换步骤如下:,1) 选择T值，设计时它是任意的，可以置为1。2) 根据截止频率p和s进行预修正，即： 3) 设计模拟滤波器。4) 最后，设 化简得到为z-1 的有理函数的H(z)。,例6：用切比雪夫I型滤波器原型设计一个低通数字滤波器,满足通带截止频率：p=0.2，通带波动：Rp=1dB阻带截止频率：s=0.3，阻带波动：As=15dB,% 数字滤波器指标wp = 0.

16、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); %预畸% 模拟切比雪夫I型原型滤波器的计算cs,ds = afd_chb1(OmegaP,OmegaS,Rp,As);

17、% 双线性变换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=filter(b,a,imp);subplot(2,2,4); pl

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

19、p2hp lp2bp lp2bs,离散化处理 bilinear impinvar,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分别表示零点、极点和增益。,MATLAB提供函数zp2tf将零极

20、点增益模型转换为传递函数模型。函数形式为：,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=ellipord(wp,ws,rp,rs

21、,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的模拟低通滤波器原型变换成截止频率为Wn的低通滤波器。

22、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）从低通到带阻的转换,lp2bs函数

23、可将截止频率为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平面建立一一对应的单值映射关系，可采用双线性

24、变换法。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

25、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

26、=50; wp1 = (2/T)*tan(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 = fre

27、qz(b,a);plot(W*Fs/(2*pi),abs(H); grid; axis(0 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 =

28、butter(N,wn,stop): 当wn为具有两个元素的矢量wn=w1 w2时， 它设计2N阶带阻滤波器，它具有单位为的3dB阻带：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，并且wp1ws1