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

1、有限脉冲响应数字滤波器设计实验报告DOC成 绩： 数字信号处理作业与上机实验（第二章）班 级： 学 号： 姓 名： 任课老师： 完成时间： 信息与通信工程学院 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

3、-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;%画出各种比较结果图figure(2);plot(wk/pi,20*log10(abs(Hk(k+1),:,linewidth,2.5);hold onplot(wk1/pi,20*log10(mag),linewidth,2);hold offlegend(FIR滤波器,II

4、R滤波器);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滤波器);xlabel(w/pi);ylabel(相位/pi);title(相频特性曲线);(2)两种数字滤波器的损耗函数和相频特性的比较分别如图1、2所示： 图1 损耗函数比较图 图2 相频特性比较图(3)IIR数字滤波器阶数：N=5 FI

5、R数字滤波器阶数：N=36(4)运行结果分析：由图2及阶数可见，IIR阶数低得多，但相位特性存在非线性失真，FIR具有线性相位特性。 20. 调用MATLAB工具箱函数fir1设计线性相位低通FIR滤波器， 要求希望逼近的理想低通滤波器通带截止频率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_han

6、ning=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,hn_hanning,:,linewidth,2);plot(n,hn_hamming,+,linewidth,2);plot(n,hn_blackman,o);hold offxlabel(n);ylabel(h(n);legend(矩形窗,汉宁窗,哈明窗,布莱克曼窗);title(单位冲激响应);M=

7、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(wk/pi,20*log10(abs(Hk(k+1),*);hold offlege

8、nd(矩形窗,汉宁窗,哈明窗,布莱克曼窗);axis(0,1,-80,5);xlabel(w/pi);ylabel(幅度/dB);title(损耗函数); (2)四种窗函数设计的单位脉冲响应的比较如图3所示：图3 单位脉冲响应比较图(3)四种窗函数设计的损耗函数的比较如图4所示：图4 损耗函数比较图 (4)运行结果分析：由图4可见，当滤波器长度N不变时，矩形窗设计的滤波器的过渡带最窄，阻带最小衰减最小；布莱克曼窗设计的滤波器的过渡带最宽，同时阻带最小衰减最大。21.将要求改成设计线性相位高通FIR滤波器，重作题20。(1)matlab代码：wc=pi/4;N=21;hn_boxcar=fir1

9、(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(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(

10、h(n);legend(矩形窗,汉宁窗,哈明窗,布莱克曼窗);title(单位冲激响应);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)

11、plot(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不变时，矩形窗设计的滤波器的过渡带最窄，阻带最小衰减最小；布莱克曼窗设计的滤波器的过渡带最宽，同时阻带最小衰减最大。25.调用MATLAB工具箱函数fir1设计

12、线性相位高通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,hamming(N); M=1024;Hk=fft(hn,M);n=0:N-1;stem(n,hn);xlabel(n);ylabel(h(n);title(单位冲激

13、响应);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滤波器的单位脉冲响应、损耗函数如图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。显示所设

14、计的单位脉冲响应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=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*

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

16、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;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

17、+1)/pi);grid onxlabel(w/pi);ylabel(相位/pi);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

18、=0;m13=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*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.01

19、45-. 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;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

20、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(幅度频谱);xlabel(频率/Hz);ylabel(幅值); (2)原信号及其幅度频谱分别如图13、14所示： 图13 信号x(n

21、)波形 图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=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,

22、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);ylabel(幅度/dB);title(损耗函数);grid onfigure(3)plot(wk/pi,angle(Hk(k

23、+1)/pi);grid onxlabel(w/pi);ylabel(相位/pi);title(相频特性曲线);(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*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

24、=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-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

25、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(原信号幅度频谱(采样点数为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,

26、abs(S);axis(0,1500,0,160);title(幅度频谱);xlabel(频率/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滤波

27、方法进行比较研究。一、低通滤波器部分：(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; 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=elli

28、pord(wp2,ws2,Rp,As);B,A=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);figure(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/