有限脉冲响应数字滤波器设计实验报告

1、成绩：数字信号处理作业与上机实验（第二章）班 级：学 号：姓 名：任课老师：完成时间：信息与通信工程学院20142015学年第1学期第7章 有限脉冲响应数字滤波器设计1、教材 p238:19 .设信号x(t) = s(t) + v(t),其中v(t)是干扰，s与v(t)的频谱不混 叠，其幅度谱如题19图所示。要求设计数字滤波器，将干扰滤除， 指标是允许|s|在0Wfwi5 kHz频率范围中幅度失真为± 2%( 5 1 = 0.02) ;f > 20 kHz,衰减大于40 dB ( 8 2=0.01);希望分别设计 性价比最高的FIR和IIR两种滤波器进行滤除干扰。 请选择合适

2、的滤波器类型和设计方法进行设计，最后比较两种滤波器的幅频特性、相频4,叔郴数。_ y/kHz题19图(1)matlab 代码：%基于双线性变换法直接设计IIR数字滤波器Fs=80000;fp=15000;fs=20000;rs=40;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Rp=-20*log10(1-0.02);As=40;N1,wp1=ellipord(wp/pi,ws/pi,Rp,As);B,A=ellip(N1,Rp,As,wp1);Hk,wk1=freqz(B,A,1000);mag=abs(Hk);pah=angle(Hk);%窗函数法设计FIR数字滤波器Bt=ws

3、-wp;alph=0.5842*(rs-21)A0.4+0.07886*(rs-21);N=ceil(rs-8)/2.285/Bt);wc=(wp+ws)/2/pi;hn=fir1(N,wc,kaiser(N+1,alph);M=1024;Hk=fft(hn,M);k=0:M/2-1;wk=(2*pi/M)*k;%画出各种比较结果图 figure(2);plot(wk/pi,20*log10(abs(Hk(k+1), hold on,'linewidth',2.5)；plot(wk1/pi,20*log10(mag), hold off'linewidth',2

4、);legend( 'FIR 滤波器，IIR axis(0,1,-80,5);xlabel( title('损耗函数');figure(3) plot(wk/pi,angle(Hk(k+1)/pi, hold on滤波器');'w/pi');ylabel( '幅度/dB','linewidth',2.5););图1损耗函数比较图图2相频特性比较图plot(wk1/pi,pah/pi, hold off'linewidth' ,2);legend( 'FIR 滤波器，IIR 滤波器);xlab

5、el( 'w/pi' );ylabel( '相位八pi');title('相频特性曲线);1、2(2)两种数字滤波器的损耗函数和相频特性的比较分别如图所示:损耗函数0-10-20FIR滤波器一 IIR滤波器4610.80.60.4O 1,20 2 4 6s-1 o-0-0-0-0/.位相FIR滤波器IIR滤波器(3)IIR数字滤波器阶数：N=5FIR数字滤波器阶数：N=36(4)运行结果分析：由图2及阶数可见，IIR阶数低得多，但相 位特性存在非线性失真，FIR具有线性相位特性。20 .调用MATLAB工具箱函数万门设计线性相位低通FIR滤波器， 要求希

6、望逼近的理想低通滤波器通带截止频率3 c=兀/4 rad,滤波器长度N=21。分别选用矩形窗、Hanning窗、Hamming窗和Blackman窗进行设计，绘制用每种窗函数设计的单位脉冲响应h(n)及其损耗函数曲线，并进行比较,观察各种窗函数的设计性台匕 目匕。(1)matlab 代码:wc=pi/4;N=21;hn_boxcar=fir1(N-1,wc/pi,boxcar(N);hn_hanning=fir1(N-1,wc/pi,hanning(N);hn_hamming=fir1(N-1,wc/pi,hamming(N);hn_blackman=fir1(N-1,wc/pi,blackm

7、an(N);n=0:N-1; plot(n,hn_boxcar);hold on plot(n,hn_hanning, plot(n,hn_hamming, plot(n,hn_blackman, hold offxlabel( 'n' );ylabel(':','linewidth''+' , 'linewidth' 'o');'h(n)');,2);,2);legend('矩形窗，汉宁窗，哈明窗，布莱克曼窗);title( '单位冲激响应');M=1024

8、;Hk=fft(hn_boxcar,M);k=0:M/2-1;wk=(2*pi/M)*k;'linewidth',2);figure();plot(wk/pi,20*log10(abs(Hk(k+1), Hk=fft(hn_hanning,M)hold onplot(wk/pi,20*log10(abs(Hk(k+1), Hk=fft(hn_hamming,M)plot(wk/pi,20*log10(abs(Hk(k+1), Hk=fft(hn_blackman,M)plot(wk/pi,20*log10(abs(Hk(k+1), hold off,'linewidth

9、''o'););,3)；图4损耗函数比较图legend('矩形窗，汉宁窗，哈明窗，布莱克曼窗)；axis(0,1,-80,5);xlabel('w/pi' );ylabel( '幅度/dB');title('损耗函数');(2)四种窗函数设计的单位脉冲响应的比较如图3所示:单位冲激响应0.3矩形窗-0.05-0.1101214161820图3 单位脉冲响应比较图四种窗函数设计的损耗函数的比较如图4所示:损耗函数0-10-20-30-40-50-60-70-8000.10.20.30.40.50.60.70.80.9

10、1w/ 7：.运行结果分析：由图4可见，当滤波器长度N不变时，矩形窗设计的滤波器的过渡带最窄，阻带最小衰减最小；布莱克曼窗设计的滤波器的过渡带最宽，同时阻带最小衰减最大。21 .将要求改成设计线性相位高通 FIR滤波器，重作题20。(1)matlab 代码:wc=pi/4;N=21;':','linewidth''+' , 'linewidth' 'o');,2);,2);hn_boxcar=fir1(N-1,wc/pi, hn_hanning=fir1(N-1,wc/pi, hn_hamming=fir1(N-1

11、,wc/pi, hn_blackman=fir1(N-1,wc/pi, n=0:N-1; plot(n,hn_boxcar); hold on plot(n,hn_hanning, plot(n,hn_hamming, plot(n,hn_blackman, hold offhigh' ,boxcar(N);'high' ,hanning(N);'high' ,hamming(N);'high' ,blackman(N)xlabel( 'n' );ylabel( 'h(n)');legend（'矩形

12、窗，汉宁窗，哈明窗，布莱克曼窗）;title（ '单位冲激响应'）;M=1024;Hk=fft（hn_boxcar,M）;k=0:M/2-1;wk=(2*pi/M)*k;'linewidth' ,2);':','linewidth' ,3);'o');figure();plot(wk/pi,20*log10(abs(Hk(k+1), Hk=fft(hn_hanning,M)hold onplot(wk/pi,20*log10(abs(Hk(k+1), Hk=fft(hn_hamming,M)plot(wk/pi,2

13、0*log10(abs(Hk(k+1), Hk=fft(hn_blackman,M)plot(wk/pi,20*log10(abs(Hk(k+1), hold off legend（'矩形窗，汉宁窗，哈明窗，布莱克曼窗）；axis（0,1,-80,5）;xlabel（'w八pi'）;ylabel（ '幅度/dB'）;title（'损耗函数'）;(2)四种窗函数设计的单位脉冲响应的比较如图5所示:单位冲激响应0.80.60.40.20-0.2*矩形窗汉宁窗*哈明窗16182002468101214n图5 单位脉冲响应比较图四种窗函数设计的损

14、耗函数的比较如图 6所示:图6损耗函数比较图布莱克曼窗(5)运行结果分析：由图6可见，当滤波器长度N不变时，矩形窗设计的滤波器的过渡带最窄，阻带最小衰减最小；布莱克曼窗 设计的滤波器的过渡带最宽，同时阻带最小衰减最大。25 .调用MATLAB工具箱函数fi门设计线性相位高通FIR滤波器。要求通带截止频率为0.6兀rad,阻带截止频率为0.45兀，通带最大衰减为0.2 dB,阻带最小衰减为45 dBo显示所设计的单位脉冲响应h(n)的数据，弁画出损耗函数曲线(1)matlab 代码：wp=0.6*pi;ws=0.45*pi;Bt=wp-ws;N0=ceil(6.6*pi/Bt);N=N0+mod

15、(N0+1,2);wc=(wp+ws)/2/pi;hn=fir1(N-1,wc, 'high' ,hamming(N);M=1024;Hk=fft(hn,M);n=0:N-1;stem(n,hn);xlabel( 'n' );ylabel( 'h(n)');title( '单位冲激响应');k=0:M/2-1;wk=(2*pi/M)*k;figure(2);plot(wk/pi,20*log10(abs(Hk(k+1););ylabel( '幅度/dB');损耗函数如图 7、8所示:axis(0,1,-80,5);

16、xlabel('w/pi'title( '损耗函数');grid on(2)高通FIR滤波器的单位脉冲响应、0.5单位冲激响应0.40.30.20.1I I-0.1-0.2-0.3-0.4图8损耗函数图7单位脉冲响应26 .调用MATLAB工具箱函数firl设计线性相位带通FIR滤波器。要求通带截止频率为 0.55兀rad和0.7兀rad,阻带截止频率为 0.45兀rad和0.8兀rad,通带最大衰减为 0.15 dB,阻带最小衰 减为40 dBo显示所设计的单位脉冲响应 h(n)的数据，弁画出损 耗函数曲线。(1)matlab 代码:wp1=0.55*pi;w

17、p2=0.7*pi;ws1=0.45*pi;ws2=0.8*pi;Bt=wp2-wp1;N=ceil(6.2*pi/Bt);wc=(wp1+ws1)/2/pi,(ws2+wp2)/2/pi;hn=fir1(N-1,wc,hanning(N);M=1024;Hk=fft(hn,M);n=0:N-1;stem(n,hn);xlabel( 'n' );ylabel( 'h(n)');title( '单位冲激响应');k=0:M/2-1;wk=(2*pi/M)*k;figure(2);plot(wk/pi,20*log10(abs(Hk(k+1););y

18、label( '幅度/dB');axis(0,1,-80,5);xlabel('w/pi'title( '损耗函数');grid on(2)带通FIR滤波器的单位脉冲响应、损耗函数如图9、10所示:1图10损耗函数图9单位脉冲响应2、 某信号 x(t)为：x(t) =0.5cos(2n f1t) +0.7cos(2兀 f2t +0.1兀)+ 0.4cos(2n f3t), 其中fi =100Hz, f2 =130Hz, fi =600Hz.设计最低阶FIR数字滤波器， 按下图所示对x(t)进行数字滤波处理，实现：1)将f3频率分量以高于50dB的

19、衰减抑制，同时以低于 2dB的衰减通过i和f2频率分量；一、基于窗函数法设计FIR数字滤波器：(1)matlab 代码：Fs=3800;fp=130;fs=600;rs=50;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Bt=ws-wp;alph=0.5842*(rs-21)A0.4+0.07886*(rs-21);N=ceil(rs-8)/2.285/Bt);wc=(wp+ws)/2/pi;hn=fir1(N,wc,kaiser(N+1,alph);M=1024;Hk=fft(hn,M);k=0:M/2-1;wk=(2*pi/M)*k;figure(2);plot(wk/pi,2

20、0*log10(abs(Hk(k+1);axis(0,1,-80,5);xlabel('w八pi');ylabel('幅度 /dB');title( '损耗函数');grid onfigure(3)plot(wk/pi,angle(Hk(k+1)/pi);gridonxlabel('w八pi');ylabel('相位八pi'); title('相频特性曲线);(2)数字滤波器的损耗函数和相频特性分别如图11、12所示：损耗函数相频特性曲线图11损耗函数曲线图12相频特性曲线二、按直接型网络结构编程编写滤波程

21、序(1)matlab 代码：N=500;n=0:N-1;f=2800;T=1/f;t=n*T;x=0.5*cos(2*pi*100*t)+0.7*cos(2*pi*130*t+0.1*pi)+0.4*cos(2*pi*600*t)；m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;m7=0;m8=0;m9=0;m10=0;m11=0;m12=0;m1 3=0;m14=0;m15=0;m16=0;m17=0;m18=0;m19=0;m20=0;m21=0;m22=0;m23=0;m24=0;for m=1:length(x) y(m)=0.0012*x(m)+m1*0.0011-m2*

22、0.0014-m3*0.0072-m4*0.0147-m5*0.019 3-.m6*0.0145+m7*0.0055+m8*0.0423+m9*0.0910+m10*0.1410+m11*0.1786+m12*0.1926+.m13*0.1786+m14*0.1410+m15*0.0910+m16*0.0423+m17*0.0055-m18*0.0145 m19*0.0193-m20*0.0147-m21*0.0072-m22*0.0014+m23*0.0011+m24*0.0012;m24=m23;m23=m22;m22=m21;m21=m20;m20=m19;m19=m18;m18=m17

23、;m17=m16;m 16=m15;m15=m14;m14=m13;m13=m12;m12=m11;m11=m10;m10=m9;m9=m8;m8=m7;m7=m6;m6=m5;m5=m4;m4=m3;m3=m2;m2=m1;m1=x(m);endplot(n,x);title( '信号x(n)');ylabel( '幅值');xlabel( 'n');S=fft(x,N);fs=n/(N*T);figure(2)plot(fs,abs(S);axis(0,1500,0,180);title('原信号幅度频谱(采样点数为500)'

24、);xlabel('频率/Hz' );ylabel( '幅值');figure(3)plot(n,y);title('信号 y(n)');ylabel('幅值');xlabel( 'n');S=fft(y,N);fs=n/(N*T);figure(4)plot(fs,abs(S);axis(0,1500,0,160);title( '幅度频谱');xlabel( '频率/Hz' );ylabel( '幅值);(2)原信号及其幅度频谱分别如图13、14所示:信号x(n)21.5

25、-2，050100150200250300350400450500n图13信号x(n)波形图14幅度频谱(3)滤波后信号y(n)及其幅度频谱分别如图15、16所示:图15信号y(n)波形图16幅度频谱2)将fi和f2频率分量以高于50dB的衰减抑制，同时以低于2dB的衰减通过f3频率分量；一、基于频率采样法设计 FIR数字滤波器：(1)matlab 代码:T=0.48;Fs=3800;fp=600;fs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; datB=wp-ws;wc=wp; m=1;N=ceil(m+1)*2*pi/datB+1); N=N+mod(N+1,2);

26、Np=fix(wc/(2*pi/N);Ns=N-2*Np-1;Ak=zeros(1,Np+1),ones(1,Ns),zeros(1,Np); Ak(Np+2)=T;Ak(N-Np尸T; thetak=-pi*(N-1)*(0:N-1)/N;Hk=Ak.*exp(1j*thetak); hn=real(ifft(Hk); M=1024; Hk=fft(hn,M); k=0:M/2-1; wk=(2*pi/M)*k; figure(2);plot(wk/pi,20*log10(abs(Hk(k+1);axis(0,1,-80,5);xlabel('w/pi'title( '

27、;损耗函数');grid on figure(3) plot(wk/pi,angle(Hk(k+1)/pi);grid xlabel('w八pi');ylabel( '相位八pi'title('相频特性曲线););ylabel( '幅度/dB');on)；(2)数字滤波器的损耗函数和相频特性分别如图17、18所示:图17损耗函数曲线图18相频特性曲线二、按直接型网络结构编程编写滤波程序：(1)matlab 代码：N=500;n=0:N-1;f=2800;T=1/f;t=n*T;x=0.5*cos(2*pi*100*t)+0.7*c

28、os(2*pi*130*t+0.1*pi)+0.4*cos(2*pi*600*t)；m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;m7=0;m8=0;m9=0;m10=0;m11=0;m12=0;m13=0;m14=0;m15=0;m16=0;for m=1:length(x)y(m)=-0.0009*x(m)-m1*0.0169-m2*0.0128+m3*0.0282+m4*0.0627+m5*0.0198- .m6*0.1198-m7*0.2827+m8*0.6447-m9*0.2827-m10*0.1198+m11*0.0198+m12*0.0627+.m13*0.0282

29、-m14*0.0128-m15*0.0169-m16*0.0009;m16=m15;m15=m14;m14=m13;m13=m12;m12=m11;m11=m10;m10=m9;m9=m8;m8=m 7;m7=m6;m6=m5;m5=m4;m4=m3;m3=m2;m2=m1;m1=x(m);endplot(n,x);title( '信号x(n)' );ylabel( '幅值');xlabel( 'n');S=fft(x,N);fs=n/(N*T);figure(2)plot(fs,abs(S);axis(0,1500,0,180);title(&

30、#39;原信号幅度频谱(采样点数为500)');xlabel('频率/Hz' );ylabel( '幅值');figure(3)plot(n,y);title( '信号 y(n)');ylabel('幅值');xlabel( 'n');S=fft(y,N);fs=n/(N*T);figure(4)plot(fs,abs(S);axis(0,1500,0,160);title( '幅度频谱)；xlabel('频率/Hz' );ylabel( '幅值');(4)原信号及其

31、幅度频谱分别如图19、20所示：信号x(n)原信号幅度频谱(采样点数为500)2180 值幅001500图19信号x(n)波形图20幅度频谱(5)滤波后信号y(n)及其幅度频谱分别如图21、22所示:幅度频谱160信号y(n)0.8-0.61111,111050100150200250300350400450500n值幅5001000频率/Hz1500图21信号y(n)波形图22幅度频谱要求：按数字滤波器直接型结构图编写滤波程序，求得 y(n)； 1)中的FIR滤波器采用窗函数法设计；2)中的FIR滤波器采用频率采样法设计。画出所设计的滤波器频率特性图、信号时域图；给出滤波器设计的 MATLA

32、B代码与滤波器实现的代码；选择合适的信号采样周期To3)与第6章作业2的IIR滤波方法进行比较研究。一、低通滤波器部分：(1)matlab 代码：Fs=3800;fp=130;fs=600;rs=50;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Bt=ws-wp;alph=0.5842*(rs-21)A0.4+0.07886*(rs-21);N=ceil(rs-8)/2.285/Bt);wc=(wp+ws)/2/pi;hn=fir1(N,wc,kaiser(N+1,alph);M=1024;Hk=fft(hn,M);k=0:M/2-1;wk=(2*pi/M)*k;wp2=2*fp/

33、Fs;ws2=2*fs/Fs;Rp=2;As=50;N1,wp1=ellipord(wp2,ws2,Rp,As);B,A=ellip(N1,Rp,As,wp1)Hk1,wk1=freqz(B,A);mag=abs(Hk1);pah=angle(Hk1);on':','linewidth',3);gridonplot(wk1/pi,20*log10(mag);grid hold onplot(wk/pi,20*log10(abs(Hk(k+1), onhold offxlabel('w八pi');ylabel(title('损耗函数曲线&#

34、39;);legend( 'IIR' , 'FIR');figure(2)plot(wk1/pi,pah/pi);gridhold onplot(wk/pi,angle(Hk(k+1)/pi, hold offxlabel('w八pi');ylabel(title('相频特性曲线)； legend( 'IIR' , 'FIR');幅度/dB')on':','linewidth' ,3);grid相位八pi');(2)两种滤波器的损耗函数、相频特性的比较图见图2

35、3、 24:损耗函数曲线加. * itI*iIIR11. VV #i.*t T* 事1* a r .K4, * .sS . : a 'tI>14-iS” ： , n-f* « H"*I* 4ft,: ii.i胃 ll.I a t-| -ft *«> 1 1 .» W以L*i1事; n =; * * «；»i-*A 1* *1 P HM n k 1 _ , k 11. * « I* i 事. . * *1 I il* r h； .* *W 4 i- r * *B .；» ,* K:JV 1I； 41-一同,! «1 p： 1 J_L二 u1 , Iw*m v1 ;* «相频特性曲线0.80.60.40.20-0.2-0.4-0.6-0.8/二 位相00.10.20.30.40.50.60.70.80.91w/ ：-1图24相频特性比较图图23损耗函数比较图IIR滤波器的阶数：N1=3FIR滤波器的阶数：N=17(1)matlab 代码：T=0.48;Fs=3800;fp=600;fs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;datB=wp-ws;wc=wp;m=1;N=ceil(m+1)