用MATLAB设计FIR数字的滤波器某实验八

1、精选优质文档-倾情为你奉上课程名称：数字信号处理 实验成绩： 指导教师： 实 验 报 告院系： 信息工程学院 班级： 学号： 姓名： 日期： 2011. 11.26 实验八：用MATLAB设计FIR数字滤波器1、 选择合适的窗函数设计FIR数字低通滤波器，要求： wp=0.2，Rp=0.05dB； ws=0.3，As=40dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。 程序清单如下：wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp; N0=ceil(6.6*pi/deltaw); N=N0+mod(N0+1,2) %为实现FIR类型1偶对称滤波器，应

2、确保N为奇数 windows=(hamming(N)' wc=(ws+wp)/2; hd=ideal_lp(wc,N); b=hd.*windows; db,mag,pha,grd,w=freqz_m(b,1); n=0:N-1;dw=2*pi/1000; Rp=-(min(db(1:wp/dw+1) %检验通带波动 As=-round(max(db(ws/dw+1:501) %检验最小阻带衰减 subplot(2,2,1);stem(n,b); axis(0,N,1.1*min(b),1.1*max(b);title('实际脉冲响应'); xlabel('n&

3、#39;);ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis(0,N,0,1.1);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db); axis(0,1,-80,10);title('幅度频率响应'); xlabel('频率（单位：pi）');ylabel('H(ejomega)'); set(gca,'XTickMode&

4、#39;,'manual','XTick',0,wp/pi,ws/pi,1); set(gca,'YTickMode','manual','YTick',-50,-20,-3,0);grid subplot(2,2,4);plot(w/pi,pha); axis(0,1,-4,4);title('相位频率响应'); xlabel('频率（单位：pi）');ylabel('phi(omega)');set(gca,'XTickMode','man

5、ual','XTick',0,wp/pi,ws/pi,1); set(gca,'YTickMode','manual','YTick',-3.1416,0,3.1416,4);grid函数 ideal_lp 调用部分如下 function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应 %wc=截止频率（弧度） %N=理想滤波器的长度 tao=(N-1)/2; n=0:(N-1); m=n-tao+eps; %加一个小数以避免0作除数hd=sin(wc*m)./(pi*m);函数 freqz_m 调用

6、部分如下functiondb,mag,pha,grd,w=freqz_m(b,a) H,w=freqz(b,a,1000,'whole'); H=(H(1:501)'w=(w(1:501)' mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle(H);grd=grpdelay(b,a,w);程序运行结果如下：N = 67Rp = 0.0394As = 522、 用凯塞窗设计一个FIR数字高通滤波器，要求： wp=0.3，Rp=0.1dB；ws=0.2，As=50dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频

7、响应曲线和相频响应曲线。程序清单如下： N=61;As=50; wp=0.3*pi;ws=0.2*pi; beta=0.112*(As-8.7) windows=kaiser(N,beta); wc=(ws+wp)/2/pi b=fir1(N-1,wc,'high',windows); db,mag,pha,grd,w=freqz_m(b,1); n=0:N-1;dw=2*pi/1000; As=-round(max(db(1:ws/dw+1) %检验最小阻带衰减 Rp=-(min(db(wp/dw+1:501) %检验通带波动 subplot(2,2,1);stem(n,b)

8、; axis(0,N,1.1*min(b),1.1*max(b);title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis(0,N,0,1.1);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db); axis(0,1,-80,10);title('幅度频率响应'); xlabel(&#

9、39;频率（单位：pi）');ylabel('H(ejomega)'); set(gca,'XTickMode','manual','XTick',0,ws/pi,wp/pi,1); set(gca,'YTickMode','manual','YTick',-50,-20,-3,0);grid subplot(2,2,4);plot(w/pi,pha); axis(0,1,-4,4);title('相位频率响应'); xlabel('频率（单位：pi）

10、');ylabel('phi(omega)'); set(gca,'XTickMode','manual','XTick',0,ws/pi,wp/pi,1); set(gca,'YTickMode','manual','YTick',-3.1416,0,3.1416,4);grid函数 freqz_m 调用部分如下functiondb,mag,pha,grd,w=freqz_m(b,a) H,w=freqz(b,a,1000,'whole'); H=(H(1:

11、501)'w=(w(1:501)' mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle(H);grd=grpdelay(b,a,w);程序运行结果如下：beta = 4.6256wc = 0.2500As = 51Rp = 0.04603、 选择合适的窗函数设计一个FIR数字带通滤波器，要求：fp1=3.5kHz，fp2=6.5kHz，Rp=0.05dB；fs1=2.5kHz，fs2=7.5kHz，As=60dB。滤波器采样频率Fs=20kHz。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。程序清单如下

12、： fp1=3500;fp2=6500; fs1=2500;fs2=7500; Fs=20000; ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi; wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi; wp=wp1,wp2;ws=ws1,ws2; deltaw=wp1-ws1; N0=ceil(11*pi/deltaw); N=N0+mod(N0+1,2) %为实现FIR类型1偶对称滤波器，应确保N为奇数 windows=(blackman(N)' wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2; hd=ideal_lp(w

13、c2,N)-ideal_lp(wc1,N); b=hd.*windows; db,mag,pha,grd,w=freqz_m(b,1); n=0:N-1;dw=2*pi/1000; Rp=-(min(db(wp1/dw+1:wp2/dw+1) %检验通带波动 ws0=1:ws1/dw+1,ws2/dw+1:501;%建立阻带频率样点数组 As=-round(max(db(ws0) %检验最小阻带衰减 subplot(2,2,1);stem(n,b); axis(0,N,1.1*min(b),1.1*max(b);title('实际脉冲响应'); xlabel('n

14、9;);ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis(0,N,0,1.1);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db); axis(0,1,-150,10);title('幅度频率响应'); xlabel('频率（单位：pi）');ylabel('H(ejomega)'); set(gca,'XTickMode&#

15、39;,'manual','XTick',0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1); set(gca,'YTickMode','manual','YTick',-100,-65,-20,-3,0);grid subplot(2,2,4);plot(w/pi,pha); axis(0,1,-4,4);title('相位频率响应'); xlabel('频率（单位：pi）');ylabel('phi(omega)'); set(gca,'XT

16、ickMode','manual','XTick',0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1);set(gca,'YTickMode','manual','YTick',-pi,0,pi);grid函数 ideal_lp 调用部分如下 function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应 %wc=截止频率（弧度） %N=理想滤波器的长度 tao=(N-1)/2; n=0:(N-1); m=n-tao+eps; %加一个小数以避免0作除数hd=sin(w

17、c*m)./(pi*m);函数 freqz_m 调用部分如下functiondb,mag,pha,grd,w=freqz_m(b,a) H,w=freqz(b,a,1000,'whole'); H=(H(1:501)'w=(w(1:501)' mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle(H);grd=grpdelay(b,a,w);程序运行结果如下：N = 111Rp = 0.0034As = 744、 选择合适的窗函数设计一个FIR数字带阻滤波器，要求：fp1=1kHz，fp2=4.5kHz，Rp=

18、0.1dB；fs1=2kHz，fs2=3.5kHz，As=40dB。滤波器采样频率Fs=10kHz。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。程序清单如下： fp1=1000;fp2=4500; fs1=2000;fs2=3500; Fs=10000; ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi; wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi; wp=wp1,wp2;ws=ws1,ws2; deltaw=ws1-wp1; N0=ceil(6.2*pi/deltaw); N=N0+mod(N0+1,2) %为实现FI

19、R类型1偶对称滤波器，应确保N为奇数 windows=(hanning(N)' wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2; hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(wc2,N);%建立理想带阻 b=hd.*windows; db,mag,pha,grd,w=freqz_m(b,1); n=0:N-1;dw=2*pi/1000; wp0=1:wp1/dw+1,wp2/dw+1:501;%建立通带频率样点数组 As=-round(max(db(ws1/dw+1:ws2/dw+1) %检验最小阻带衰减 Rp=-(min(db

20、(wp0) %检验通带波动 subplot(2,2,1);stem(n,b); axis(0,N,1.1*min(b),1.1*max(b);title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis(0,N,0,1.1);title('窗函数特性'); xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db); axis(0,1,-150

21、,10);title('幅度频率响应'); xlabel('频率（单位：pi）');ylabel('H(ejomega)'); set(gca,'XTickMode','manual','XTick',0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1); set(gca,'YTickMode','manual','YTick',-150,-40,-3,0);grid subplot(2,2,4);plot(w/pi,pha); axis(

22、0,1,-4,4);title('相位频率响应'); xlabel('频率（单位：pi）');ylabel('phi(omega)'); set(gca,'XTickMode','manual','XTick',0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1);set(gca,'YTickMode','manual','YTick',-pi,0,pi);grid函数 ideal_lp 调用部分如下 function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应 %wc=截止频率（弧度） %N=理想滤波器的长度 tao=(N-1)/2; n=0:(N-1); m=n-tao+eps; %加一个小数以避免0作除数hd=sin(wc*m)./(pi*m);函数 freqz_m 调用部分如下functiondb,mag,pha,grd,w=freqz_m(b,a) H,w=freqz(b,a,1000,'whole'); H=(H(1:501)'w=(w(1:501)' mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle