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

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

2、波器类型和设计方法进行设计，最后比较两种滤波器的幅频特性、相频特性和阶数。 题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-wp; alph

3、=0.5842*(rs-21)0.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),':','linewidth',2.5);hold onplot(wk1/pi,20*log10(mag),'linewidth',2);hol

4、d offlegend('FIR滤波器','IIR滤波器');axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数');figure(3)plot(wk/pi,angle(Hk(k+1)/pi,':','linewidth',2.5);hold onplot(wk1/pi,pah/pi,'linewidth',2);hold offlegend('FIR滤波器','IIR滤波器&#

5、39;);xlabel('w/pi');ylabel('相位/pi');title('相频特性曲线');(2)两种数字滤波器的损耗函数和相频特性的比较分别如图1、2所示： 图1 损耗函数比较图图2 相频特性比较图(3) IIR数字滤波器阶数：N=5 FIR数字滤波器阶数：N=36(4) 运行结果分析：由图2及阶数可见，IIR阶数低得多，但相位特性存在非线性失真，FIR具有线性相位特性。 20. 调用MATLAB工具箱函数fir1设计线性相位低通FIR滤波器， 要求希望逼近的理想低通滤波器通带截止频率c=/4 rad， 滤波器长度N=21。分别选用

6、矩形窗、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,blackman(N);n=0:N-1;plot(n,hn_boxcar);hold onplot(n,h

7、n_hanning,':','linewidth',2);plot(n,hn_hamming,'+','linewidth',2);plot(n,hn_blackman,'o');hold offxlabel('n');ylabel('h(n)');legend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');title('单位冲激响应');M=1024;Hk=fft(hn_boxcar,M);k=

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

9、39;*');hold offlegend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数'); (2)四种窗函数设计的单位脉冲响应的比较如图3所示：图3 单位脉冲响应比较图(3)四种窗函数设计的损耗函数的比较如图4所示：图4 损耗函数比较图 (4) 运行结果分析：由图4可见，当滤波器长度N不变时，矩形窗设计的滤波器的过渡带最窄，阻带最小衰减最小；布莱克

10、曼窗设计的滤波器的过渡带最宽，同时阻带最小衰减最大。21.将要求改成设计线性相位高通FIR滤波器，重作题20。(1)matlab代码：wc=pi/4;N=21;hn_boxcar=fir1(N-1,wc/pi,'high',boxcar(N);hn_hanning=fir1(N-1,wc/pi,'high',hanning(N);hn_hamming=fir1(N-1,wc/pi,'high',hamming(N);hn_blackman=fir1(N-1,wc/pi,'high',blackman(N);n=0:N-1;plot

11、(n,hn_boxcar);hold onplot(n,hn_hanning,':','linewidth',2);plot(n,hn_hamming,'+','linewidth',2);plot(n,hn_blackman,'o');hold offxlabel('n');ylabel('h(n)');legend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');title('单位冲激响应');

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

13、wk/pi,20*log10(abs(Hk(k+1),'*');hold offlegend('矩形窗','汉宁窗','哈明窗','布莱克曼窗');axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数'); (2)四种窗函数设计的单位脉冲响应的比较如图5所示：图5 单位脉冲响应比较图(3)四种窗函数设计的损耗函数的比较如图6所示：图6 损耗函数比较图(5) 运行结果分析：由图6可见，当滤波器长度N不变时，

14、矩形窗设计的滤波器的过渡带最窄，阻带最小衰减最小；布莱克曼窗设计的滤波器的过渡带最宽，同时阻带最小衰减最大。25. 调用MATLAB工具箱函数fir1设计线性相位高通FIR滤波器。 要求通带截止频率为0.6 rad，阻带截止频率为0.45，通带最大衰减为0.2 dB，阻带最小衰减为45 dB。显示所设计的单位脉冲响应h(n)的数据，并画出损耗函数曲线。(1)matlab代码：wp=0.6*pi;ws=0.45*pi;Bt=wp-ws; N0=ceil(6.6*pi/Bt); N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi; hn=fir1(N-1,wc,'high&#

15、39;,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);axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数');grid on (2)高通FIR滤波器的单位脉冲响应、

16、损耗函数如图7、8所示： 图7 单位脉冲响应图8 损耗函数26. 调用MATLAB工具箱函数fir1设计线性相位带通FIR滤波器。 要求通带截止频率为0.55 rad和0.7 rad，阻带截止频率为0.45 rad和0.8 rad，通带最大衰减为0.15 dB，阻带最小衰减为40 dB。显示所设计的单位脉冲响应h(n)的数据，并画出损耗函数曲线。 (1)matlab代码：wp1=0.55*pi;wp2=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

17、=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);axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数');grid on (2)带通FIR滤波

18、器的单位脉冲响应、损耗函数如图9、10所示： 图9 单位脉冲响应 图10 损耗函数2、某信号为：，其中设计最低阶FIR数字滤波器，按下图所示对进行数字滤波处理，实现：1） 将频率分量以高于50dB的衰减抑制，同时以低于2dB的衰减通过和频率分量；一、基于窗函数法设计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)0.4+0.07886*(rs-21); N=ceil(rs-8)/2.285/Bt); wc=(wp+ws)/2/pi;

19、 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);axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数');grid onfigure(3)plot(wk/pi,angle(Hk(k+1)/pi);grid onxlabel('w/pi');ylabel('相位/pi'

20、);title('相频特性曲线');(2)数字滤波器的损耗函数和相频特性分别如图11、12所示：图11 损耗函数曲线图12 相频特性曲线二、按直接型网络结构编程编写滤波程序:(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;m13=0;m14=0;m15=0;m16=0;m17=

21、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*0.0014-m3*0.0072-m4*0.0147-m5*0.0193-. 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*

22、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;m17=m16;m16=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=ff

23、t(x,N);fs=n/(N*T);figure(2)plot(fs,abs(S);axis(0,1500,0,180);title('原信号幅度频谱(采样点数为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('

24、;幅度频谱');xlabel('频率/Hz');ylabel('幅值'); (2) 原信号及其幅度频谱分别如图13、14所示： 图13 信号x(n)波形图14 幅度频谱(3) 滤波后信号y(n)及其幅度频谱分别如图15、16所示： 图15 信号y(n)波形图16 幅度频谱2） 将和频率分量以高于50dB的衰减抑制，同时以低于2dB的衰减通过频率分量；一、基于频率采样法设计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

25、=wp; m=1;N=ceil(m+1)*2*pi/datB+1); N=N+mod(N+1,2); 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)

26、;axis(0,1,-80,5);xlabel('w/pi');ylabel('幅度/dB');title('损耗函数');grid onfigure(3)plot(wk/pi,angle(Hk(k+1)/pi);grid onxlabel('w/pi');ylabel('相位/pi');title('相频特性曲线');(2)数字滤波器的损耗函数和相频特性分别如图17、18所示：图17 损耗函数曲线图18 相频特性曲线二、按直接型网络结构编程编写滤波程序:(1)matlab代码：N=500;n=0:

27、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;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*

28、0.2827-m10*0.1198+m11*0.0198+m12*0.0627+. m13*0.0282-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=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);f

29、igure(2)plot(fs,abs(S);axis(0,1500,0,180);title('原信号幅度频谱(采样点数为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(

30、'频率/Hz');ylabel('幅值');(4) 原信号及其幅度频谱分别如图19、20所示： 图19 信号x(n)波形 图20 幅度频谱(5) 滤波后信号y(n)及其幅度频谱分别如图21、22所示： 图21 信号y(n)波形图22 幅度频谱要求：按数字滤波器直接型结构图编写滤波程序，求得；1)中的FIR滤波器采用窗函数法设计；2)中的FIR滤波器采用频率采样法设计。画出所设计的滤波器频率特性图、信号时域图；给出滤波器设计的MATLAB代码与滤波器实现的代码；选择合适的信号采样周期T。3） 与第6章作业2的IIR滤波方法进行比较研究。一、低通滤波器部分：(1)m

31、atlab代码：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)0.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/Fs;ws2=2*fs/Fs;Rp=2;As=50;N1,wp1=ellipord(wp2,ws2,Rp,As);B,A

32、=ellip(N1,Rp,As,wp1)Hk1,wk1=freqz(B,A);mag=abs(Hk1);pah=angle(Hk1); plot(wk1/pi,20*log10(mag);grid onhold onplot(wk/pi,20*log10(abs(Hk(k+1),':','linewidth',3);grid onhold offxlabel('w/pi');ylabel('幅度/dB')title('损耗函数曲线');legend('IIR','FIR');figu

33、re(2)plot(wk1/pi,pah/pi);grid onhold onplot(wk/pi,angle(Hk(k+1)/pi,':','linewidth',3);grid onhold offxlabel('w/pi');ylabel('相位/pi');title('相频特性曲线');legend('IIR','FIR');(2) 两种滤波器的损耗函数、相频特性的比较图见图23、24：图23 损耗函数比较图图24 相频特性比较图(3) IIR滤波器的阶数：N1=3 FIR滤波器的阶数：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)*2*pi/datB+1); N=N+mod(N+1,2); Np=fix(wc/(2*pi/N);Ns=N-2*Np-1;Ak=zeros(1,Np+1),ones(1,Ns),zer