有限脈沖響應(yīng)數(shù)字濾波器設(shè)計(jì)實(shí)驗(yàn)報(bào)告

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

2、的濾波器類型和設(shè)計(jì)方法進(jìn)行設(shè)計(jì)，最后比較兩種濾波器的幅頻特性、相頻4,叔郴數(shù)。_ y/kHz題19圖(1)matlab 代碼：%基于雙線性變換法直接設(shè)計(jì)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è)計(jì)FIR數(shù)字濾波器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;%畫出各種比較結(jié)果圖 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('損耗函數(shù)');figure(3) plot(wk/pi,angle(Hk(k+1)/pi, hold on濾波器');'w/pi');ylabel( '幅度/dB','linewidth',2.5););圖1損耗函數(shù)比較圖圖2相頻特性比較圖plot(wk1/pi,pah/pi, hold off'linewidth' ,2);legend( 'FIR 濾波器，IIR 濾波器);xlab

5、el( 'w/pi' );ylabel( '相位八pi');title('相頻特性曲線);1、2(2)兩種數(shù)字濾波器的損耗函數(shù)和相頻特性的比較分別如圖FIR濾波器IIR濾波器8664202 o o O-04 6 80 0。一/,位相FIR濾波器一 IIR濾波器665OWB夕度幅(3)IIR數(shù)字濾波器階數(shù)：N=5FIR數(shù)字濾波器階數(shù)：N=36(4)運(yùn)行結(jié)果分析：由圖2及階數(shù)可見，IIR階數(shù)低得多，但相 位特性存在非線性失真，F(xiàn)IR具有線性相位特性。20 .調(diào)用MATLAB工具箱函數(shù)萬門設(shè)計(jì)線性相位低通FIR濾波器， 要求希望逼近的理想低通濾波器通帶截止頻率

6、3 c=兀/4 rad,濾波器長(zhǎng)度N=21。分別選用矩形窗、Hanning窗、Hamming窗和Blackman窗進(jìn)行設(shè)計(jì)，繪制用每種窗函數(shù)設(shè)計(jì)的單位脈沖響應(yīng)h(n)及其損耗函數(shù)曲線，并進(jìn)行比較,觀察各種窗函數(shù)的設(shè)計(jì)性臺(tái)匕 目匕。(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; pl

7、ot(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( '單位沖激響應(yīng)');M=1024;Hk=fft(hn_boxcar

8、,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''o')

9、；,3)；圖4損耗函數(shù)比較圖legend('矩形窗，漢寧窗，哈明窗，布萊克曼窗)；axis(0,1,-80,5);xlabel('w/pi' );ylabel( '幅度/dB');title('損耗函數(shù)');(2)四種窗函數(shù)設(shè)計(jì)的單位脈沖響應(yīng)的比較如圖3所示:0.3單位沖激響應(yīng)0.250.20.150.10.050-0.05-0.1矩形窗漢寧窗 哈明窗布萊克曼窗n圖3 單位脈沖響應(yīng)比較圖四種窗函數(shù)設(shè)計(jì)的損耗函數(shù)的比較如圖4所示:損耗函數(shù)0-10-20-30-40-50-60-70-8000.10.20.30.40.50.60.70.80.

10、91w/ 7：.運(yùn)行結(jié)果分析：由圖4可見，當(dāng)濾波器長(zhǎng)度N不變時(shí)，矩形窗設(shè)計(jì)的濾波器的過渡帶最窄，阻帶最小衰減最小；布萊克曼窗設(shè)計(jì)的濾波器的過渡帶最寬，同時(shí)阻帶最小衰減最大。21 .將要求改成設(shè)計(jì)線性相位高通 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-

11、1,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（ '單位沖激響應(yīng)'）;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,

13、20*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（'損耗函數(shù)'）;(2)四種窗函數(shù)設(shè)計(jì)的單位脈沖響應(yīng)的比較如圖5所示:單位沖激響應(yīng)0.80.60.40.20-0.2*矩形窗漢寧窗*哈明窗16182002468101214n圖5 單位脈沖響應(yīng)比較圖四種窗函數(shù)設(shè)計(jì)的

14、損耗函數(shù)的比較如圖 6所示:圖6損耗函數(shù)比較圖布萊克曼窗(5)運(yùn)行結(jié)果分析：由圖6可見，當(dāng)濾波器長(zhǎng)度N不變時(shí)，矩形窗設(shè)計(jì)的濾波器的過渡帶最窄，阻帶最小衰減最??；布萊克曼窗 設(shè)計(jì)的濾波器的過渡帶最寬，同時(shí)阻帶最小衰減最大。25 .調(diào)用MATLAB工具箱函數(shù)fi門設(shè)計(jì)線性相位高通FIR濾波器。要求通帶截止頻率為0.6兀rad,阻帶截止頻率為0.45兀，通帶最大衰減為0.2 dB,阻帶最小衰減為45 dBo顯示所設(shè)計(jì)的單位脈沖響應(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+mo

15、d(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( '單位沖激響應(yīng)');k=0:M/2-1;wk=(2*pi/M)*k;figure(2);plot(wk/pi,20*log10(abs(Hk(k+1););ylabel( '幅度/dB');損耗函數(shù)如圖 7、8所示:axis(0,1,-80,5)

16、;xlabel('w/pi'title( '損耗函數(shù)');grid on(2)高通FIR濾波器的單位脈沖響應(yīng)、0.5單位沖激響應(yīng)0.40.30.20.1I I-0.1-0.2-0.3-0.4圖8損耗函數(shù)圖7單位脈沖響應(yīng)26 .調(diào)用MATLAB工具箱函數(shù)firl設(shè)計(jì)線性相位帶通FIR濾波器。要求通帶截止頻率為 0.55兀rad和0.7兀rad,阻帶截止頻率為 0.45兀rad和0.8兀rad,通帶最大衰減為 0.15 dB,阻帶最小衰 減為40 dBo顯示所設(shè)計(jì)的單位脈沖響應(yīng) h(n)的數(shù)據(jù)，弁畫出損 耗函數(shù)曲線。(1)matlab 代碼:wp1=0.55*pi;

17、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( '單位沖激響應(yīng)');k=0:M/2-1;wk=(2*pi/M)*k;figure(2);plot(wk/pi,20*log10(abs(Hk(k+1););

18、ylabel( '幅度/dB');axis(0,1,-80,5);xlabel('w/pi'title( '損耗函數(shù)');grid on(2)帶通FIR濾波器的單位脈沖響應(yīng)、損耗函數(shù)如圖9、10所示:1圖10損耗函數(shù)圖9單位脈沖響應(yīng)2、 某信號(hào) x(t)為：x(t) =0.5cos(2n f1t)+0.7cos(2兀 f2t+0.1兀) + 0.4cos(2n f3t), 其中fi =100Hz, f2 =130Hz, fi =600Hz.設(shè)計(jì)最低階FIR數(shù)字濾波器， 按下圖所示對(duì)x(t)進(jìn)行數(shù)字濾波處理，實(shí)現(xiàn)：1)將f3頻率分量以高于50dB的

19、衰減抑制，同時(shí)以低于 2dB的衰減通過i和f2頻率分量；一、基于窗函數(shù)法設(shè)計(jì)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)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( '損耗函數(shù)');grid onfigure(3)plot(wk/pi,angle(Hk(k+1)/pi);gridonxlabel('w八pi');ylabel('相位八pi'); title('相頻特性曲線);(2)數(shù)字濾波器的損耗函數(shù)和相頻特性分別如圖11、12所示：損耗函數(shù)相頻特性曲線w/ ：,圖11損耗函數(shù)曲線圖12相頻特性曲線二、按直接型網(wǎng)絡(luò)結(jié)構(gòu)編程

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.001

22、1-m2*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;m1

23、8=m17;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( '信號(hào)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('原信號(hào)幅度頻譜(采樣點(diǎn)數(shù)為500)

24、');xlabel('頻率/Hz' );ylabel( '幅值);figure(3)plot(n,y);title('信號(hào) 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)原信號(hào)及其幅度頻譜分別如圖13、14所示:信號(hào)x(n)21.5

25、-2，050100150200250300350400450500n圖13信號(hào)x(n)波形圖14幅度頻譜(3)濾波后信號(hào)y(n)及其幅度頻譜分別如圖15、16所示:圖15信號(hào)y(n)波形圖16幅度頻譜2)將fi和f2頻率分量以高于50dB的衰減抑制，同時(shí)以低于2dB的衰減通過f3頻率分量；一、基于頻率采樣法設(shè)計(jì) 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=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、;損耗函數(shù)');grid on figure(3) plot(wk/pi,angle(Hk(k+1)/pi);grid xlabel('w八pi');ylabel( '相位八pi'title('相頻特性曲線););ylabel( '幅度/dB');on)；(2)數(shù)字濾波器的損耗函數(shù)和相頻特性分別如圖17、18所示:圖17損耗函數(shù)曲線圖18相頻特性曲線二、按直接型網(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*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( '信號(hào)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;原信號(hào)幅度頻譜(采樣點(diǎn)數(shù)為500)');xlabel('頻率/Hz' );ylabel( '幅值);figure(3)plot(n,y);title( '信號(hào) 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)原信號(hào)及其幅度頻譜分別如圖19

31、、20所示：信號(hào)x(n)原信號(hào)幅度頻譜(采樣點(diǎn)數(shù)為500)2180 值幅001500圖19信號(hào)x(n)波形圖20幅度頻譜(5)濾波后信號(hào)y(n)及其幅度頻譜分別如圖21、22所示:幅度頻譜160信號(hào)y(n)0.8-0.61111,111050100150200250300350400450500n值幅5001000頻率/Hz1500圖21信號(hào)y(n)波形圖22幅度頻譜要求：按數(shù)字濾波器直接型結(jié)構(gòu)圖編寫濾波程序，求得 y(n)； 1)中的FIR濾波器采用窗函數(shù)法設(shè)計(jì)；2)中的FIR濾波器采用頻率采樣法設(shè)計(jì)。畫出所設(shè)計(jì)的濾波器頻率特性圖、信號(hào)時(shí)域圖；給出濾波器設(shè)計(jì)的 MATLAB代碼與濾波器實(shí)現(xiàn)的

32、代碼；選擇合適的信號(hào)采樣周期To3)與第6章作業(yè)2的IIR濾波方法進(jìn)行比較研究。一、低通濾波器部分：(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/Fs;ws2=2*f

33、s/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('損耗函數(shù)曲線)； legend( &

34、#39;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)兩種濾波器的損耗函數(shù)、相頻特性的比較圖見圖23、 24:損耗函數(shù)曲線加.

35、 * 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損耗函數(shù)比較圖IIR濾波器的階數(shù)：N1=3FIR濾波器的階數(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)*