数字信号处理设计IIR滤波器北华大学

1、 北华大学开放实验报告 实验项目名称：设计IIR滤波器 所属课程名称：数字信号处理 实验类型:设计型 实验日期: 班级:通信 学号: 姓名一实验目的1) 加深对脉冲响应不变法以及双线性变换法设计IIR数字滤波器基本方法的了解。2) 掌握使用模拟滤波器原型进行脉冲响应变换的方法以及掌握用双线性变换法设计数字低通，高通，带通，带阻滤波器的方法。3) 了解MATLAB有关脉冲响应变换及双线性变换法的子函数。二实验步骤脉冲响应不变法设计步骤l 确定数字滤波器性能指标l 将数字滤波器频率指标转换成相应的模拟滤波器频率指标 l 根据指标设计模拟滤波器将展成部分分式形式l 把模拟极点转换成数字极点，得到数字

2、滤波器可见至间的变换关系为在MATLAB中有两种方法可以实现上述变换。方法：利用residue函数和residuez函数实现脉冲响应不变法，这两个函数的使用方法如下：r,p,k=residue(b,a)b,a=residue(r,p,k)实现多项式形式和部分分式形式之间的转换。r,p,k=residuez(b,a)b,a=residuez(r,p,k)实现多项式形式和部分分式形式之间的转换双线性变换法的设计步骤 输入给定的数字滤波器设计指标 将数字滤波器指标转换成模拟滤波器设计指标 确定模拟滤波器的最小阶数和截止频率 设计模拟低通原型滤波器的系统传递函数 利用模拟域频率变换法，求解实际模拟滤波

3、器的系统传递函数 用双线性变换法将模拟滤波器转换成数字滤波器三实验内容1、采用脉冲响应不变法设计一个巴特沃斯数字低通滤波器，要求：Wp=0.25，Rp=1dB ;Ws=0.4，As=15dB ,滤波器采样频率Fs=2000Hz。编写的程序如下：%数字滤波器指标 wp=0.25*pi; %滤波器的通带截至频率ws=0.4*pi; %滤波器的阻带截至频率Rp=1;As=15; %输入滤波器的通阻带衰减指标ripple=10(-Rp/20); %计算通带衰减对应的幅度值Attn=10(-As/20); %计算阻带衰减对应的幅度值%转换为模拟滤波器指标Fs=2000;T=1/Fs;Omgp=wp*Fs

4、;Omgs=ws*Fs;%模拟原型滤波器计算n,Omgc=buttord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率z0,p0,k0=buttap(n); %设计归一化的巴特沃斯模拟原型滤波器ba1=k0*real(poly(z0); %求原型滤波器系数baa1=real(poly(p0); %求原型滤波器系数aba,aa=lp2lp(ba1,aa1,Omgc); %变换为模拟低通滤波器%用脉冲响应不变法计算数字滤波器系数bd,ad=impinvar(ba,aa,Fs)%C,B,A=dir2par(bd,ad)H,w=freqz(bd,ad);dbH=20

5、*log10(abs(H)+eps)/max(abs(H);subplot(2,2,1),plot(w/pi,abs(H);ylabel('H');title('幅度响应');axis(0,1,0,1.1);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','XTick',0,Attn,ripple,1);gridsubplot(2,2,2),plot(w/p

6、i,angle(H)/pi);ylabel('phi');title('相位响应');axis(0,1,-1,1);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','XTick',-1,0,1);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)')ylabel('dB');xl

7、abel('频率(pi)');axis(0,1,-40,5);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','XTick',-50.-15.-1.0);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);title('零极图');频率特性图20-1：由频率特性曲线可知，该设计结果在通阻带截止频率

8、处能满足Rp1dB,As15dB的设计指标要求，且系统的极点全部在单位圆内，是一个稳定的系统。2、 采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求:Wp=0.25，Rp=1dB ;Ws=0.4，As=15dB ,滤波器采样频率Fs=1000Hz。编写的程序如下：%数字滤波器指标 wp=0.25*pi; %滤波器的通带截至频率ws=0.4*pi; %滤波器的阻带截至频率Rp=1;As=15; %输入滤波器的通阻带衰减指标ripple=10(-Rp/20);Attn=10(-As/20);Fs=100;T=1/Fs;Omgp=(2/T)*tan(wp/2);Omgs=(2/T)*tan(ws

9、/2);n,Omgc=buttord(Omgp,Omgs,Rp,As,'s')z0,p0,k0=buttap(n);ba=k0*real(poly(z0);aa=real(poly(p0);ba1,aa1=lp2lp(ba,aa,Omgc);bd,ad=bilinear(ba1,aa1,Fs)sos,g=tf2sos(bd,ad)H,w=freqz(bd,ad);dbH=20*log10(abs(H)+eps)/max(abs(H);subplot(2,2,1),plot(w/pi,abs(H);ylabel('H');title('幅度响应')

10、;axis(0,1,0,1.1);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','XTick',0,Attn,ripple,1);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi);ylabel('phi');title('相位响应');axis(0,1,-1,1);set(gca,'XTickMode',&

11、#39;manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','XTick',-1,0,1);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)')ylabel('dB');xlabel('频率(pi)');axis(0,1,-40,5);set(gca,'XTickMode','manual','XTick',0,0

12、.25,0.4,1);set(gca,'YTickMode','manual','XTick',-50.-15.-1.0);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);title('零极图');频率特性图如下21-1：由频率特性曲线21-1可知，该设计结果在通组带截止处能满足Rp，1dB、As20dB的设计指标要求。由n=3可知，由3阶的模拟低通原型用双线性变换法设计出来的切比雪夫I型数字带通滤波器是一个6阶的系统，极点全部在Z平面的单位圆内，是一个稳定的系统。

13、3、采用脉冲响应不变法法设计一个巴特沃斯数字带通滤波器，要求:通带Wp1=0.3，Wp2=0.7，Rp=1dB ;阻带Ws1=0.1，Ws2=0.9，As=15dB ,滤波器采样频率Fs=,2000Hz。试显示数字滤波器的幅频特性和零极点分布图，并写出该系统的传递函数。编写的程序如下：wp1=0.3*pi;wp2=0.7*pi;ws1=0.1*pi;ws2=0.9*pi;Rp=1;As=15;Fs=2000;T=1/Fs;Omgp1=wp1*Fs;Omgp2=wp2*Fs;Omgp=Omgp1,Omgp2;Omgs1=ws1*Fs;Omgs2=ws2*Fs;Omgs=Omgs1,Omgs2;b

14、w=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2);n,Omgn=buttord(Omgp,Omgs,Rp,As,'s')z0,p0,k0=buttap(n);ba1=k0*real(poly(z0);aa1=real(poly(p0);ba,aa=lp2bp(bal,aal,w0,bw)bd,ad=impinvar(ba,aa,Fs)H,w=freqz(bd,ad);dbH=20*log10(abs(H)+eps)/max(abs(H);subplot(2,2,1),plot(w/pi,abs(H);subplot(2,2,2),plot(w/pi,angl

15、e(H)/pi);subplot(2,2,3),plot(w/pi,dbH);subplot(2,2,4),zplane(bd,ad); 频率特性图如下20-24、用双线型变换法设计切比雪夫II型数字滤波器，列出传递函数并描绘模拟和数字滤波器的幅频响应曲线。 (1)、设计一个数字低通，要求：通带wp=0.2，Rp=1dB;阻带ws=0.35，As=15dB,滤波器采样频率Fs=10Hz. (2)设计一个数字高通，要求：通带wp=0.35，Rp=1dB;阻带 ws=0.2，As=15dB,滤波器采样频率Fs=10Hz.编写的程序如下： (1) wp=0.2*pi;ws=0.35*pi;Rp=1;

16、As=15;ripple=10(-Rp/20);Attn=10(-As/20);Fs=10;T=1/Fs;Omgp=(2/T)*tan(wp/2);Omgs=(2/T)*tan(ws/2);n,Omgc=cheb2ord(Omgp,Omgs,Rp,As,'s')z0,p0,k0=cheb2ap(n,As);ba=k0*real(poly(z0);aa=real(poly(p0);ba1,aa1=lp2lp(ba,aa,Omgc);bd,ad=bilinear(ba1,aa1,Fs)sos,g=tf2sos(bd,ad)H,w=freqz(bd,ad);dbH=20*log10(

17、abs(H)+eps)/max(abs(H);subplot(2,2,1),plot(w/pi,abs(H);ylabel('|H|');title('幅度响应');axis(0,1,0,1.1);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','YTick',0,Attn,ripple,1);gridsubplot(2,2,2),plot(w/pi,ang

18、le(H)/pi);ylabel('phi');title('相位响应');axis(0,1,-1,1);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','YTick',-1,0,1);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)');ylabel('dB');xlabel

19、('频率');axis(0,1,-40,5);set(gca,'XTickMode','manual','XTick',0,0.25,0.4,1);set(gca,'YTickMode','manual','YTick',-50,-15,-1,0);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);title('零极图'); n = 3 Omgc =10.2282 bd =0.1624 0.0056 0.

20、0056 0.1624 ad =1.0000 -1.4073 0.9056 -0.1623 sos = 1.0000 1.0000 0 1.0000 -0.2720 0 1.0000 -0.9658 1.0000 1.0000 -1.1353 0.5968 g =0.1624 频率特性图如下,21-2-1(2) Fs=10;T=1/Fs; wp=0.35*pi; ws=0.2*pi;fp=wp/(2*pi)*Fs; fs=ws/(2*pi)*Fs; Rp=1;As=15;ripple=10(-Rp/20);Attn=10(-As/20);Omgp=(2/T)*tan(wp/2);Omgs=(2

21、/T)*tan(ws/2);n,Omgc=cheb2ord(Omgp,Omgs,Rp,As,'s') z0,p0,k0=cheb2ap(n,As); ba=k0*real(poly(z0); aa=real(poly(p0);ba1,aa1=lp2hp(ba,aa,Omgc); bd,ad=bilinear(ba1,aa1,Fs) H,w=freqz(bd,ad);dbH=20*log10(abs(H)+eps)/max(abs(H); subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');ylabel('|H|')

22、;title('幅度响应');axis(0,Fs/2,0,1.1); set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2);set(gca,'YTickMode','manual','YTick',0,Attn,ripple,1);grid subplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');ylabel('phi');title('相位响应');axis(0,Fs/2,-180,180); set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2); set(gca,'YTickMode','manual','YTick',-180,0,180);gri