有限脈沖響應(yīng)數(shù)字濾波器設(shè)計實驗報告

1、成 績： 數(shù)字信號處理作業(yè)與上機實驗（第二章）班 級： 學(xué) 號： 姓 名： 任課老師： 完成時間： 信息與通信工程學(xué)院 20142015學(xué)年第 1 學(xué)期第7章 有限脈沖響應(yīng)數(shù)字濾波器設(shè)計1、教材p238：19.設(shè)信號x(t) = s(t) + v(t),其中v(t)是干擾，s(t)與v(t)的頻譜不混疊，其幅度譜如題19圖所示。要求設(shè)計數(shù)字濾波器,將干擾濾除，指標(biāo)是允許|s(f)|在0f15 kHz頻率范圍中幅度失真為±2%（1 = 0.02）;f > 20 kHz,衰減大于40 dB（2=0.01）;希望分別設(shè)計性價比最高的FIR和IIR兩種濾波器進(jìn)行濾除干擾。請選擇合適的濾

2、波器類型和設(shè)計方法進(jìn)行設(shè)計，最后比較兩種濾波器的幅頻特性、相頻特性和階數(shù)。 題19圖(1)matlab代碼：%基于雙線性變換法直接設(shè)計IIR數(shù)字濾波器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);%窗函數(shù)法設(shè)計FIR數(shù)字濾波器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;%畫出各種比較結(jié)果圖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('損耗函數(shù)');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)兩種數(shù)字濾波器的損耗函數(shù)和相頻特性的比較分別如圖1、2所示： 圖1 損耗函數(shù)比較圖圖2 相頻特性比較圖(3) IIR數(shù)字濾波器階數(shù)：N=5 FIR數(shù)字濾波器階數(shù)：N=36(4) 運行結(jié)果分析：由圖2及階數(shù)可見，IIR階數(shù)低得多，但相位特性存在非線性失真，F(xiàn)IR具有線性相位特性。 20. 調(diào)用MATLAB工具箱函數(shù)fir1設(shè)計線性相位低通FIR濾波器， 要求希望逼近的理想低通濾波器通帶截止頻率c=/4 rad， 濾波器長度N=21。分別選用

6、矩形窗、Hanning窗、Hamming窗和Blackman窗進(jìn)行設(shè)計，繪制用每種窗函數(shù)設(shè)計的單位脈沖響應(yīng)h(n)及其損耗函數(shù)曲線，并進(jìn)行比較，觀察各種窗函數(shù)的設(shè)計性能。 (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('單位沖激響應(yīng)');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('損耗函數(shù)'); (2)四種窗函數(shù)設(shè)計的單位脈沖響應(yīng)的比較如圖3所示：圖3 單位脈沖響應(yīng)比較圖(3)四種窗函數(shù)設(shè)計的損耗函數(shù)的比較如圖4所示：圖4 損耗函數(shù)比較圖 (4) 運行結(jié)果分析：由圖4可見，當(dāng)濾波器長度N不變時，矩形窗設(shè)計的濾波器的過渡帶最窄，阻帶最小衰減最??；布萊克

10、曼窗設(shè)計的濾波器的過渡帶最寬，同時阻帶最小衰減最大。21.將要求改成設(shè)計線性相位高通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('單位沖激響應(yīng)');

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('損耗函數(shù)'); (2)四種窗函數(shù)設(shè)計的單位脈沖響應(yīng)的比較如圖5所示：圖5 單位脈沖響應(yīng)比較圖(3)四種窗函數(shù)設(shè)計的損耗函數(shù)的比較如圖6所示：圖6 損耗函數(shù)比較圖(5) 運行結(jié)果分析：由圖6可見，當(dāng)濾波器長度N不變時，

14、矩形窗設(shè)計的濾波器的過渡帶最窄，阻帶最小衰減最?。徊既R克曼窗設(shè)計的濾波器的過渡帶最寬，同時阻帶最小衰減最大。25. 調(diào)用MATLAB工具箱函數(shù)fir1設(shè)計線性相位高通FIR濾波器。 要求通帶截止頻率為0.6 rad，阻帶截止頻率為0.45，通帶最大衰減為0.2 dB，阻帶最小衰減為45 dB。顯示所設(shè)計的單位脈沖響應(yīng)h(n)的數(shù)據(jù)，并畫出損耗函數(shù)曲線。(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('單位沖激響應(yīng)');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('損耗函數(shù)');grid on (2)高通FIR濾波器的單位脈沖響應(yīng)、

16、損耗函數(shù)如圖7、8所示： 圖7 單位脈沖響應(yīng)圖8 損耗函數(shù)26. 調(diào)用MATLAB工具箱函數(shù)fir1設(shè)計線性相位帶通FIR濾波器。 要求通帶截止頻率為0.55 rad和0.7 rad，阻帶截止頻率為0.45 rad和0.8 rad，通帶最大衰減為0.15 dB，阻帶最小衰減為40 dB。顯示所設(shè)計的單位脈沖響應(yīng)h(n)的數(shù)據(jù)，并畫出損耗函數(shù)曲線。 (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('單位沖激響應(yīng)');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('損耗函數(shù)');grid on (2)帶通FIR濾波

18、器的單位脈沖響應(yīng)、損耗函數(shù)如圖9、10所示： 圖9 單位脈沖響應(yīng) 圖10 損耗函數(shù)2、某信號為：，其中設(shè)計最低階FIR數(shù)字濾波器，按下圖所示對進(jìn)行數(shù)字濾波處理，實現(xiàn)：1） 將頻率分量以高于50dB的衰減抑制，同時以低于2dB的衰減通過和頻率分量；一、基于窗函數(shù)法設(shè)計FIR數(shù)字濾波器：(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('損耗函數(shù)');grid onfigure(3)plot(wk/pi,angle(Hk(k+1)/pi);grid onxlabel('w/pi');ylabel('相位/pi'

20、);title('相頻特性曲線');(2)數(shù)字濾波器的損耗函數(shù)和相頻特性分別如圖11、12所示：圖11 損耗函數(shù)曲線圖12 相頻特性曲線二、按直接型網(wǎng)絡(luò)結(jié)構(gòu)編程編寫濾波程序:(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('原信號幅度頻譜(采樣點數(shù)為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的衰減通過頻率分量；一、基于頻率采樣法設(shè)計FIR數(shù)字濾波器：(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('損耗函數(shù)');grid onfigure(3)plot(wk/pi,angle(Hk(k+1)/pi);grid onxlabel('w/pi');ylabel('相位/pi');title('相頻特性曲線');(2)數(shù)字濾波器的損耗函數(shù)和相頻特性分別如圖17、18所示：圖17 損耗函數(shù)曲線圖18 相頻特性曲線二、按直接型網(wǎng)絡(luò)結(jié)構(gòu)編程編寫濾波程序:(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('原信號幅度頻譜(采樣點數(shù)為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 幅度頻譜要求：按數(shù)字濾波器直接型結(jié)構(gòu)圖編寫濾波程序，求得；1)中的FIR濾波器采用窗函數(shù)法設(shè)計；2)中的FIR濾波器采用頻率采樣法設(shè)計。畫出所設(shè)計的濾波器頻率特性圖、信號時域圖；給出濾波器設(shè)計的MATLAB代碼與濾波器實現(xiàn)的代碼；選擇合適的信號采樣周期T。3） 與第6章作業(yè)2的IIR濾波方法進(jìn)行比較研究。一、低通濾波器部分：(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('損耗函數(shù)曲線');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) 兩種濾波器的損耗函數(shù)、相頻特性的比較圖見圖23、24：圖23 損耗函數(shù)比較圖圖24 相頻特性比較圖(3) IIR濾波器的階數(shù)：N1=3 FIR濾波器的階數(shù)：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