FIR數(shù)字濾波器的Matlab實(shí)現(xiàn)

1、第7章 FIR數(shù)字濾波器的Matlab實(shí)現(xiàn)7.1 實(shí)驗(yàn)?zāi)康膌 學(xué)習(xí)用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的原理及其設(shè)計(jì)步驟；l 學(xué)習(xí)編寫數(shù)字濾波器的設(shè)計(jì)程序的方法，并能進(jìn)行正確編程；l 根據(jù)給定的濾波器指標(biāo)，給出設(shè)計(jì)步驟。 7.2 實(shí)驗(yàn)原理及實(shí)例分析 FIR低通數(shù)字濾波器的設(shè)計(jì)原理如果系統(tǒng)的沖激響應(yīng)為已知，則系統(tǒng)的輸入/輸出關(guān)系為： 對于低通濾波器，只要設(shè)計(jì)出低通濾波器的沖激響應(yīng)函數(shù)，就可以由上式得到系統(tǒng)的輸出了。假設(shè)所希望的數(shù)字濾波器的頻率響應(yīng)為，它是頻域的周期函數(shù)，周期為2，那么它與相對應(yīng)的傅立葉系數(shù)為以為沖激響應(yīng)的數(shù)字濾波器將具有頻域響。但是將作為濾波器脈沖響應(yīng)有兩個(gè)問題：（1） 它是無限長的，

2、與FIP濾波器脈沖響應(yīng)有限長這一前提不一致（2） 它是非因果的，對此，要采取以下的措施，（1） 將截短 （2）將其往右平移，由此得到的實(shí)際頻域響應(yīng)，與理想頻域響應(yīng)相近，但不完全一致。理論證明上述現(xiàn)象是對進(jìn)行簡單截短處理的必然結(jié)果，一般稱為吉布斯現(xiàn)象，為盡可能的減少吉布斯現(xiàn)象，應(yīng)對進(jìn)行加窗截取，即以作為FIR濾波器的系數(shù)。常用的窗函數(shù)有矩形窗、海明窗和布萊克曼窗等。 用窗函數(shù)法設(shè)計(jì)FIR濾波器Matlab設(shè)計(jì)FIR濾波器有多種方法和對應(yīng)的函數(shù)，見表71。表71 matlab設(shè)計(jì)FIR濾波器的方法和函數(shù)方法描述函數(shù)窗方法使用窗函數(shù)和逆傅立葉變換實(shí)現(xiàn)fir1,fir2,kaiserord等多帶方法包

3、含子帶頻率域firls，remez等最小二乘法使用最小二乘法將整個(gè)頻率域上的錯(cuò)誤幾率壓縮到最小fircls ，fircls1等任意響應(yīng)法使用任意響應(yīng)，包括非線性相位以及復(fù)濾波器cremez等余弦法使用三角函數(shù)的低通響應(yīng)firrcos等窗函數(shù)方法不僅在數(shù)字濾波器的設(shè)計(jì)中占有重要的地位，同時(shí)可以用于功率譜的估計(jì)，從根本上講，使用窗函數(shù)的目的就是消除由無限序列的截短而引起的Gibbs現(xiàn)象所帶來的影響。窗函數(shù)設(shè)計(jì)線性相位FIR濾波器步驟如下：（1） 確定數(shù)字濾波器的性能要求，臨界頻率，濾波器單位脈沖響應(yīng)長度N（2） 根據(jù)性能要求，合理選擇單位脈沖響應(yīng)h(n)的奇偶對稱性，從而確定理想頻率響應(yīng)的幅頻特性

4、和相頻特性（3） 求理想單位脈沖響應(yīng)，在實(shí)際計(jì)算中，可對采樣，并對其求IDFT的，用代替（4） 選擇適當(dāng)?shù)拇昂瘮?shù)w（n），根據(jù)求所需設(shè)計(jì)的FIR濾波器單位脈沖響應(yīng)（5） 求，分析其幅頻特性，若不滿足要求，可適當(dāng)改變窗函數(shù)形式或長度N，重復(fù)上述設(shè)計(jì)過程，以得到滿意的結(jié)果FIR1.m調(diào)用格式：a) b=fir1(N,wn)b) b=fir1(N,wn,hign)N為濾波器的階次，wn是通帶截止頻率，其值在01之間，1對應(yīng)抽樣頻率的一半；b是設(shè)計(jì)好的濾波器系數(shù)h(n)。對于格式a)，若wn是標(biāo)量，則用來設(shè)計(jì)低通濾波器，若wn是1×L的向量，則用來設(shè)計(jì)帶通(帶阻)濾波器。其格式為：b=fir

5、1(N,wn,DC-1)或b=fir1(N,wn,DC-0),其中，前者保證第一個(gè)帶為通帶，后者保證第一個(gè)帶為阻帶。b）用來設(shè)計(jì)高通濾波器。用于產(chǎn)生窗函數(shù)的matlab文件有如下8個(gè)：bartlett.m(巴特利特窗)blackman.m(布萊克曼窗)rectangle.m(矩形窗)hamming.m（漢明窗）hanning.m（漢寧窗）triang.m(三角窗)chebwin.m（切比雪夫窗）Kaiser.m（凱澤窗）【實(shí)例7-1】 設(shè)計(jì)一個(gè)34階的高通濾波器，截止頻率為0.48，使用具有30dB波紋的chebyshev窗。解：源程序如下： b=fir1(34,0.48,'high&

6、#39;,chebwin(35,30);freqz(b,1,512)其響應(yīng)波形如圖7-1所示。圖71 帶通FIR濾波器【實(shí)例7-2】設(shè)計(jì)一個(gè)24階FIR帶通濾波器，通帶頻率。解：源程序如下：wn=0.35 0.65;N=24;b=fir1(2*N,wn,'DC-0');% DC-0表示第一個(gè)帶為阻帶，這樣設(shè)計(jì)出來的濾波器為帶通濾波器。freqz(b,1,512)圖72 帶通FIR濾波器【實(shí)例7-3】 設(shè)計(jì)具有下面指標(biāo)的低通FIR濾波器由于其最小阻帶衰減為50dB，因此可以選擇hamming窗來實(shí)現(xiàn)這個(gè)濾波器，因?yàn)樗哂休^小的過渡帶。解：MATLAB源程序?yàn)?數(shù)字濾波器指標(biāo)wp=

7、0.2*pi;ws=0.3*pi;tr_width=ws-wp; %通帶寬度N=ceil(6.6*pi/tr_width)+1; %濾波器長度if rem(N,2)=0 N=N+1;endNw=N; %濾波器長度為奇數(shù)n=0:1:Nw-1;wc=(ws+wp)/2; %截止頻率alpha=(Nw-1)/2;m=n-alpha+0.00001; %如果不加0.00001，當(dāng)n= alpha 時(shí)，m為0，下式分母為0，不能執(zhí)行hd=sin(wc*m)./(pi*m); %理想濾波器的沖激響應(yīng)% 生成hamming窗w_ham=(hamming(Nw)'% 頻域圖像的繪制h=hd.*w_ha

8、m; % 濾波器實(shí)際沖激響應(yīng)freqz (h,1)figure(2);subplot(2,2,1),stem(n,hd);title('理想脈沖響應(yīng)')axis(0 Nw-1 -0.3 0.3);xlabel('n');ylabel('hd(n)')subplot(2,2,2),stem(n,w_ham);title('hamming窗')axis(0 Nw-1 -0.3 1.2);xlabel('n');ylabel('w(n)')subplot(2,2,3),stem(n,h);title(&#

9、39;實(shí)際脈沖響應(yīng)')axis(0 Nw-1 -0.3 0.3);xlabel('n');ylabel('h(n)')其響應(yīng)波形如圖7-2所示。 圖73 hamming窗函數(shù)設(shè)計(jì)FIR濾波器【實(shí)例7-4】設(shè)帶通濾波器的指標(biāo)為選擇Blackman窗來實(shí)現(xiàn)這個(gè)濾波器。解：MATLAB源程序?yàn)? 數(shù)字濾波器指標(biāo)ws1=0.2*pi; wp1=0.35*pi;ws2=0.65*pi; wp2=0.8*pi;As=60;tr_width=min(wp1-ws1),(wp2-ws2); % 過渡帶帶寬N=ceil(11*pi/tr_width)+1;if rem(N

10、,2)=0 N=N+1;endNw=N;n=0:1:Nw-1;wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;alpha=(Nw-1)/2;m=n-alpha+0.00001; %如果不加0.00001，當(dāng)n= alpha 時(shí)，m為0，下式分母為0，不能執(zhí)行hd=sin(wc2*m)./(pi*m)- sin(wc1*m)./(pi*m); %理想濾波器的沖激響應(yīng)% 生成blackman窗w_bla=(blackman(M)'h=hd.*w_bla;% 頻域圖像的繪制freqz(h,1)figure(2);subplot(2,2,1),stem(n,hd);title(

11、'idael impulse response')axis(0 Nw-1 -0.4 0.5);xlabel('n');ylabel('hd(n)')subplot(2,2,2),stem(n,w_bla);title('blackman window')axis(0 Nw-1 0 1.1);xlabel('n');ylabel('w(n)')subplot(2,2,3),stem(n,h);title('actual impulse response')axis(0 Nw-1 -0.4 0.5);xlabel('n');ylabel('h(n)')其響應(yīng)波形如圖7-4所示。圖74 bl