用MATLAB設(shè)計(jì)FIR數(shù)字濾波器

1、實(shí)驗(yàn)八 用MATLAB設(shè)計(jì)FIR數(shù)字濾波器（二）一、實(shí)驗(yàn)?zāi)康模?、加深對窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的基本原理的理解。2、學(xué)習(xí)用MATLAB語言的窗函數(shù)法編寫設(shè)計(jì)FIR數(shù)字濾波器的程序。3、了解MATLAB語言有關(guān)窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的常用函數(shù)用法。二、實(shí)驗(yàn)原理：1、用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器FIR數(shù)字濾波器的系統(tǒng)函數(shù)為這個(gè)公式也可以看成是離散LSI系統(tǒng)的系統(tǒng)函數(shù)分母a0為1，其余ak全都為0時(shí)的一個(gè)特例。由于極點(diǎn)全部集中在零點(diǎn)，穩(wěn)定和線性相位特性是FIR濾波器的突出優(yōu)點(diǎn)，因此在實(shí)際中廣泛使用。FIR濾波器的設(shè)計(jì)任務(wù)是選擇有限長度的h(n)，使傳輸函數(shù)H(ej)滿足技術(shù)要求。主要設(shè)

2、計(jì)方法有窗函數(shù)法、頻率采樣法和切比雪夫等波紋逼近法等。本實(shí)驗(yàn)主要介紹窗函數(shù)法。用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的基本步驟如下：（1）根據(jù)過渡帶和阻帶衰減指標(biāo)選擇窗函數(shù)的類型，估算濾波器的階數(shù)N。（2）由數(shù)字濾波器的理想頻率響應(yīng)H(ej)求出其單位脈沖響應(yīng)hd(n)?？捎米远x函數(shù)ideal_lp實(shí)現(xiàn)理想數(shù)字低通濾波器單位脈沖響應(yīng)的求解。程序清單如下：function hd=ideal_lp(wc,N) %點(diǎn)0到N-1之間的理想脈沖響應(yīng)%wc=截止頻率（弧度）%N=理想濾波器的長度tao=(N-1)/2;n=0:(N-1);m=n-tao+eps; %加一個(gè)小數(shù)以避免0作除數(shù)hd=sin(wc*m

3、)./(pi*m);其它選頻濾波器可以由低通頻響特性合成。如一個(gè)通帶在c1c2之間的帶通濾波器在給定N值的條件下，可以用下列程序?qū)崿F(xiàn)：Hd=ideal_lp(wc2,N)-ideal_lp(wc1,N)（3）計(jì)算數(shù)字濾波器的單位沖激響應(yīng)h(n)=w(n)hd(n)。（4）檢查設(shè)計(jì)的濾波器是否滿足技術(shù)指標(biāo)。如果設(shè)計(jì)的濾波器不滿足技術(shù)指標(biāo)，則需要重新選擇或調(diào)整窗函數(shù)的類型，估算濾波器的階數(shù)N。再重復(fù)前面的四個(gè)步驟，直到滿足指標(biāo)。常用的窗函數(shù)有矩形窗、三角形窗、漢寧窗、哈明窗、切比雪夫窗、布萊克曼窗、凱塞窗等，MATLAB均有相應(yīng)的函數(shù)可以調(diào)用。另外，MATLAB信號處理工具箱還提供了fir1函數(shù)，

4、可以用于窗函數(shù)法設(shè)計(jì)FIR濾波器。由于第一類線性相位濾波器（類型）能進(jìn)行低通、高通、帶通、帶阻濾波器的設(shè)計(jì)，因此，本實(shí)驗(yàn)所有濾波器均采用第一類線性相位濾波器。三、實(shí)驗(yàn)內(nèi)容：1、閱讀并輸入實(shí)驗(yàn)原理中介紹的例題程序，觀察輸出的數(shù)據(jù)和圖形，結(jié)合基本原理理解每一條語句的含義。2、選擇合適的窗函數(shù)設(shè)計(jì)FIR數(shù)字低通濾波器，要求：p=0.2，Rp=0.05dB；s=0.3，As=40dB。描繪該濾波器的脈沖響應(yīng)、窗函數(shù)及濾波器的幅頻響應(yīng)曲線和相頻響應(yīng)曲線。解：程序清單如下：wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp;N0=ceil(6.6*pi/deltaw);N=N0+mod(N0

5、+1,2) %為實(shí)現(xiàn)FIR類型1偶對稱濾波器，應(yīng)確保N為奇數(shù)windows=(hamming(N);wc=(ws+wp)/2;%截止頻率hd=ideal_lp(wc,N);b=hd.*windows;db,mag,pha,grd,w=freqz_m(b,1);%dbn=0:N-1;dw=2*pi/1000;Rp=-(min(db(1:wp/dw+1) %檢驗(yàn)通帶波動(dòng)As=-round(max(db(ws/dw+1:501) %檢驗(yàn)最小阻帶衰減subplot(2,2,1);stem(n,b);axis(0,N,1.1*min(b),1.1*max(b);title(實(shí)際脈沖響應(yīng));xlabel(

6、n);ylabel(h(n);subplot(2,2,2);stem(n,windows);axis(0,N,0,1.1);title(窗函數(shù)特性);xlabel(n);ylabel(wd(n);subplot(2,2,3);plot(w/pi,db);axis(0,1,-80,10);title(幅度頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(H(ejomega);set(gca,XTickMode,manual,XTick,0,wp/pi,ws/pi,1);set(gca,YTickMode,manual,YTick,-50,-20,-3,0);gridsubplot(2

7、,2,4);plot(w/pi,pha);axis(0,1,-4,4);title(相位頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(phi(omega);set(gca,XTickMode,manual,XTick,0,wp/pi,ws/pi,1);set(gca,YTickMode,manual,YTick,-3.1416,0,3.1416,4);grid程序運(yùn)行結(jié)果： N =67 Rp =0.0394 As =52波形如下：3、 用凱塞窗設(shè)計(jì)一個(gè)FIR數(shù)字高通濾波器，要求：p=0.3，Rp=0.1dB；s=0.2，As=50dB。描繪該濾波器的脈沖響應(yīng)、窗函數(shù)及濾波器的幅

8、頻響應(yīng)曲線和相頻響應(yīng)曲線。解：程序清單如下：wp=0.3*pi;ws=0.2*pi;deltaw=wp-ws;N0=ceil(6.1*pi/deltaw);N=N0+mod(N0+1,2) %為實(shí)現(xiàn)FIR類型1偶對稱濾波器，應(yīng)確保N為奇數(shù)windows=(kaiser(N);wc=(ws+wp)/2;hd=ideal_lp(pi,N)-ideal_lp(wc,N);b=hd.*windows;db,mag,pha,grd,w=freqz_m(b,1);n=0:N-1;dw=2*pi/1000;Rp=-(min(db(wp/dw+1:501) %檢驗(yàn)通帶波動(dòng)As=-round(max(db(1:

9、ws/dw+1) %檢驗(yàn)最小阻帶衰減subplot(2,2,1);stem(n,b);axis(0,N,1.1*min(b),1.1*max(b);title(實(shí)際脈沖響應(yīng));xlabel(n);ylabel(h(n);subplot(2,2,2);stem(n,windows);axis(0,N,0,1.1);title(窗函數(shù)特性);xlabel(n);ylabel(wd(n);subplot(2,2,3);plot(w/pi,db);axis(0,1,-40,2);title(幅度頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(H(ejomega);set(gca,XTic

10、kMode,manual,XTick,0,ws/pi,wp/pi,1);set(gca,YTickMode,manual,YTick,-20,-3,0);gridsubplot(2,2,4);plot(w/pi,pha);axis(0,1,-4,4);title(相位頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(phi(omega);set(gca,XTickMode,manual,XTick,0,ws/pi,wp/pi,1);set(gca,YTickMode,manual,YTick,-pi,0,pi);grid 運(yùn)行結(jié)果： N =61 Rp = 1.1756 As =29

11、波形如下：4、 選擇合適的窗函數(shù)設(shè)計(jì)一個(gè)FIR數(shù)字帶通濾波器，要求：fp1=3.5kHz，fp2=6.5kHz，Rp=0.05dB；fs1=2.5kHz，fs2=7.5kHz，As=60dB。濾波器采樣頻率Fs=20kHz。描繪該濾波器的脈沖響應(yīng)、窗函數(shù)及濾波器的幅頻響應(yīng)曲線和相頻響應(yīng)曲線。解：程序清單如下：fp1=3500;fp2=6500;fs1=2500;fs2=7500;Fs=20000;ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi;wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi;deltaw=wp1-ws1;N0=ceil(11*pi

12、/deltaw);N=N0+mod(N0+1,2)%為實(shí)現(xiàn)FIR類型1偶對稱濾波器，應(yīng)確保N為奇數(shù)windows=blackman(N);wc1=(ws1+wp1)/2/pi;wc2=(ws2+wp2)/2/pi;b=fir1(N-1,wc1,wc2,windows);db,mag,pha,grd,w=freqz_m(b,1);n=0:N-1;dw=2*pi/1000;Rp=-(min(db(wp1/dw+1:wp2/dw+1)%檢驗(yàn)通帶波動(dòng)ws0=1:ws1/dw+1,ws2/dw+1:501;%建立阻帶頻率樣點(diǎn)數(shù)組As=-round(max(db(ws0)%檢驗(yàn)最小阻帶衰減subplot(

13、2,2,1);stem(n,b);axis(0,N,1.1*min(b),1.1*max(b);title(實(shí)際脈沖響應(yīng));xlabel(n);ylabel(h(n);subplot(2,2,2);stem(n,windows);axis(0,N,0,1.1);title(窗函數(shù)特性);xlabel(n);ylabel(wd(n);subplot(2,2,3);plot(w/pi,db);axis(0,1,-150,10);title(幅度頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(H(ejomega);set(gca,XTickMode,manual,XTick,0,fs1

14、,fp1,fp2,fs2,8000);set(gca,YTickMode,manual,YTick,-150,-40,-3,0);gridsubplot(2,2,4);plot(w/pi,pha);axis(0,1,-4,4);title(相位頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(phi(omega);set(gca,XTickMode,manual,XTick,0,fs1,fp1,fp2,fs2,500);set(gca,YTickMode,manual,YTick,-pi,0,pi);grid運(yùn)行結(jié)果： N = 111 Rp =0.0034 As =74波形圖如下：

15、5、 選擇合適的窗函數(shù)設(shè)計(jì)一個(gè)FIR數(shù)字帶阻濾波器，要求：fp1=1kHz，fp2=4.5kHz，Rp=0.1dB；fs1=2kHz，fs2=3.5kHz，As=40dB。濾波器采樣頻率Fs=10kHz。描繪該濾波器的脈沖響應(yīng)、窗函數(shù)及濾波器的幅頻響應(yīng)曲線和相頻響應(yīng)曲線。解：程序清單如下：fs1=2000;fs2=3500;fp1=1000;fp2=4500;Fs=10000;ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi;wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi;deltaw=ws1-wp1;N0=ceil(6.1*pi/deltaw);

16、N=N0+mod(N0+1,2) windows=(hanning(N);wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(wc2,N);%建立理想帶阻b=hd.*windows;db,mag,pha,grd,w=freqz_m(b,1);n=0:N-1;dw=2*pi/1000;wp0=1:wp1/dw+1,wp2/dw+1:501;%建立通帶頻率樣點(diǎn)數(shù)組As=-round(max(db(ws1/dw+1:ws2/dw+1) %檢驗(yàn)最小阻帶衰減Rp=-(min(db(wp0) %檢驗(yàn)通帶波動(dòng)s

17、ubplot(2,2,1);stem(n,b);axis(0,N,1.1*min(b),1.1*max(b);title(實(shí)際脈沖響應(yīng));xlabel(n);ylabel(h(n);subplot(2,2,2);stem(n,windows);axis(0,N,0,1.1);title(窗函數(shù)特性);xlabel(n);ylabel(wd(n);subplot(2,2,3);plot(w/pi,db);axis(0,1,-150,10);title(幅度頻率響應(yīng));xlabel(頻率（單位：pi）);ylabel(H(ejomega);set(gca,XTickMode,manual,XTick,0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1);set(gca,YTickMode,