matlab的fir高通數(shù)字濾波器的設(shè)計(jì)及分析

1、摘要無限長脈沖數(shù)字濾波器的設(shè)計(jì)方法只考慮了幅度特性，沒有考慮相位特性，所設(shè)的濾波器一般是某種確定的非線性相位特性。有限脈沖響應(yīng)(FIR)濾波器在保證了幅度特性滿足技術(shù)要求的同時(shí)，很容易做到有嚴(yán)格的線性相位特性。本課題利用MATLAB軟件實(shí)現(xiàn)。MATLAB是“矩陣實(shí)驗(yàn)室”(MATrixLABoratoy)的縮寫，是一種科學(xué)計(jì)算軟件，它使用方便，輸入簡捷，運(yùn)算高效，內(nèi)容豐富，因此利用MATLAB軟件，通過一系列較為系統(tǒng)的函數(shù)法，根據(jù)已知的技術(shù)指標(biāo)，就可以設(shè)計(jì)出滿足要求的濾波器。關(guān)鍵字：MATLAB;窗函數(shù)；FIR帶阻數(shù)字濾波器；線性相位目錄1 .FIR濾波器簡介31.1 FIR的特點(diǎn)3.2.2線性

2、相位3.2 .主要設(shè)計(jì)內(nèi)容53 .窗函數(shù)63.1 常用窗函數(shù)63.2 窗函數(shù)的指標(biāo)94應(yīng)用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的步驟104.1數(shù)字高通濾波器的設(shè)計(jì)：10總結(jié)11參考文獻(xiàn)12附錄131.FIR濾波器簡介數(shù)字濾波器是一種用來過濾時(shí)間離散信號(hào)的數(shù)字系統(tǒng)，通過對(duì)抽樣數(shù)據(jù)進(jìn)行數(shù)學(xué)處理來達(dá)到頻域?yàn)V波的目的。根據(jù)其單位沖激響應(yīng)函數(shù)的時(shí)域特性可分為兩類：無限沖激響應(yīng)（IIR）濾波器和有限沖激響應(yīng)（FIR）濾波器。1.1FIR的特點(diǎn)示為FIR濾波器的主要優(yōu)點(diǎn)為：系統(tǒng)總是穩(wěn)定的，F(xiàn)IR濾波器的系統(tǒng)函數(shù)可以表f(z)NzN1NH(z)="'h(n)z“書N)“h(n)zN"n為n

3、=0易知，H（z）在Z平面上有N-1個(gè)零點(diǎn)，z=0是N-1階極點(diǎn)，因此FIR系統(tǒng)總是穩(wěn)定的（極點(diǎn)都在單位圓內(nèi)）。FIR濾波器的優(yōu)點(diǎn)之二：容易實(shí)現(xiàn)線性相位。當(dāng)FIR系統(tǒng)的單位沖激響應(yīng)滿足h（n）=±h（N-1-n）時(shí)，該系統(tǒng)具有線性相位。中=-®（N12（N為奇數(shù)）（2-2）中（8）=n/2co（N1）/2（N為偶數(shù)）（2-3）FIR濾波器的優(yōu)點(diǎn)之三：允許設(shè)置多通帶（或多阻帶）濾波器。FIR濾波器的優(yōu)點(diǎn)之四：FIR濾波器可以采用FFT方法實(shí)現(xiàn)其功能，從而大大提高效率。FIR濾波器的缺點(diǎn)：由于FIR系統(tǒng)只有零點(diǎn)，因此這類系統(tǒng)不像FIR濾波器不像IIR濾波器那樣容易取得比較好的通

4、帶與阻帶衰減特性。要取得較好的衰減特性，一般要求H（z）的階次較高。綜合起來看，F(xiàn)IR濾波器具有IIR濾波器沒有的許多特點(diǎn)，得到了越來越廣泛的應(yīng)用。FIR濾波器的設(shè)計(jì)方法主要有三種：a.窗函數(shù)設(shè)計(jì)法；b.頻率抽樣發(fā)；c.最小平法抽樣法；這里我主要討論在MATLAB環(huán)境下通過調(diào)用信號(hào)分析與處理工具箱的幾類窗函數(shù)來設(shè)計(jì)濾波器并分析與比較其性能2.2線性相位一個(gè)單一頻率的正弦信號(hào)通過一個(gè)系統(tǒng)，假設(shè)它通過這個(gè)系統(tǒng)的時(shí)間需要t,則這個(gè)信號(hào)的輸出相位落后原來信號(hào)wt的相位。從這邊可以看出，一個(gè)正弦信號(hào)通過一個(gè)系統(tǒng)落后的相位等于它的w*t;反過來說，如果一個(gè)頻率為w的正弦信號(hào)通過系統(tǒng)后，它的相位落后delt

5、a,則該信號(hào)被延遲了delta/w的時(shí)間。在實(shí)際系統(tǒng)中，一個(gè)輸入信號(hào)可以分解為多個(gè)正弦信號(hào)的疊加，為了使得輸出信號(hào)不會(huì)產(chǎn)生相位失真，必須要求它所包含的這些正弦信號(hào)通過系統(tǒng)的時(shí)間是一樣的。因此每一個(gè)正弦信號(hào)的相位分別落后，w1*t,w2*t,w3*t。因此，落后的相位正比于頻率w,如果超前，超前相位的大小也是正比于頻率Wo從系統(tǒng)的頻率響應(yīng)來看，就是要求它的相頻特性是一條直線。在FIR濾波器的設(shè)計(jì)中，為了得到線性相位的性質(zhì)，通常利用實(shí)偶對(duì)稱序列的相頻特性為常數(shù)0和實(shí)奇對(duì)稱序列為相頻特性為常數(shù)90度的特點(diǎn)。因此得到的是對(duì)稱序列，不是因果序列，是不可實(shí)現(xiàn)系統(tǒng)，為了稱為物理可實(shí)現(xiàn)系統(tǒng)，需要將它向右移動(dòng)半

6、個(gè)周期，這就造成了相移特性隨時(shí)間的變化，同時(shí)也是線性變化。單位脈沖響應(yīng)h(n)(為實(shí)數(shù))具有偶對(duì)稱或奇對(duì)稱性，則FIR數(shù)字濾波器具有嚴(yán)格的線性相位特性。2 .主要設(shè)計(jì)內(nèi)容利用窗函數(shù)法、頻率取樣法及優(yōu)化設(shè)計(jì)方法設(shè)計(jì)FIR濾波器，繪制出濾波器的特性圖。利用所設(shè)計(jì)的濾波器對(duì)多個(gè)頻帶疊加的正弦信號(hào)進(jìn)行處理，對(duì)比濾波前后的信號(hào)時(shí)域和頻域圖，驗(yàn)證濾波器的效果?；舅悸罚簭臅r(shí)域出發(fā)設(shè)計(jì)h(n)逼近理想hd(n)0設(shè)理想濾波器的單位響應(yīng)在時(shí)域表達(dá)為hd(n),則Hd(n)一般是無限長的，且是非因果的，不能直接作為FIR濾波器的單位脈沖響應(yīng)。要想得到一個(gè)因果的有限長的濾波器單位抽樣響應(yīng)h(n),最直接的方法是先

7、將hd(n)往右平移，再邇行截?cái)?，即截取為有限長因果序列：h(n)=hd(n)w(n),并用合適的窗函數(shù)邇行加權(quán)作為FIR濾波器的單位脈沖響應(yīng)。按照線性相位濾波器的要求，線性相位FIR數(shù)字低通濾波器的單位抽樣響應(yīng)h(n)必須是偶對(duì)稱的。對(duì)稱中心必須等于濾波器的延時(shí)常數(shù)，即用矩形窗設(shè)計(jì)的FIR低通濾波器，所設(shè)計(jì)濾波器的幅度函數(shù)在通帶和阻帶都呈現(xiàn)出振蕩現(xiàn)象，且最大波紋大約為幅度的9%,返個(gè)現(xiàn)象稱為吉布斯(Gibbs)效應(yīng)。為了消除吉布斯效應(yīng)，一般采用其他類型的窗函數(shù)。MATLAB設(shè)計(jì)FIR濾波器有多種方法和對(duì)應(yīng)的函數(shù)。窗函數(shù)設(shè)計(jì)法不僅在數(shù)字濾波器的設(shè)計(jì)中占有重要的地位，同時(shí)可以用于功率譜的估計(jì)，從

8、根本上講，使用窗函數(shù)的目的就是消除由無限序列的截短而引起的Gibbs現(xiàn)象所帶來的影響。3 .窗函數(shù)加窗處理使得得濾波器的頻率響應(yīng)與理想濾波器的頻率響應(yīng)之間產(chǎn)生差異，表現(xiàn)為過渡帶和波動(dòng)的出現(xiàn)。我們希望所設(shè)計(jì)的濾波器盡量逼近理想濾波器，就要設(shè)法減少波動(dòng)的幅度，同時(shí)使過渡帶變窄。在設(shè)計(jì)FIR數(shù)字濾波器時(shí)，窗函數(shù)的頻譜應(yīng)該滿足：1主瓣寬度盡可能的窄，以使過渡帶盡量陡峭；2最大旁瓣相對(duì)于主瓣盡可能的小，使能量盡可能集中于主瓣內(nèi)，這樣能夠使得波動(dòng)減小。3.1 常用窗函數(shù)1 漢寧(Hanning)窗漢寧窗函數(shù)的時(shí)域形式表示為：(3-1)NeF'(3-2)2n二w(n)=0.5-0.5cos(),n=

9、0,1,2,N-1N-1頻域形式為W(j®)電0.5Wr+0.25.|Wr儂-瓦)+Wr儂+)漢寧窗函數(shù)的最大旁瓣值比主瓣值低31dB,但是主瓣寬度比矩形窗函數(shù)的主瓣寬度增加了1倍，為8tt/N。hanning函數(shù)：生成漢寧窗調(diào)用方式：(1) w=hanning(n)：輸入?yún)?shù)n是窗函數(shù)的長度；輸出參數(shù)w是由窗函數(shù)的值組成的n階向量。注意：此函數(shù)不返回是零點(diǎn)的窗函數(shù)的首尾兩個(gè)元素。(2) w=hanning(n,'symmetric'):與上面相類似。(3) w=hanning(n,'periodic'):此函數(shù)返回包括為零點(diǎn)的窗函數(shù)的首尾兩個(gè)元素。0

10、0.51圖3-1漢寧窗及其頻譜特性注釋：漢寧窗又被稱為升余弦窗，漢寧窗可以看做三個(gè)矩形時(shí)間窗的頻譜之和，漢寧窗主瓣加寬并降低，旁瓣顯著減小，分辨率下降2漢明(Hamming)窗：函數(shù)的時(shí)域形式可以表示為w(k)=0.540.46cosl2jk=1,2,，N(3-3),N-1頻域形式為：w(o)=0.54Wr(。)+0.23.WrA-2+W6+-2jLJ(3-4)_.N-1,N-1其中，WR(附為矩形窗函數(shù)的幅度頻率特性函數(shù)。海明窗函數(shù)的最大旁瓣值比主瓣值低41dB,但它和漢寧窗函數(shù)的主瓣寬度是一樣大的。Hamming函數(shù)：生成海明窗調(diào)用方式(1) w=hamming(n)：輸入?yún)?shù)n是窗函數(shù)的

11、長度；輸出參數(shù)w是由窗函數(shù)的值組成的n階向量。(2) w=hamming(n,sflag):參數(shù)sflag用來控制窗函數(shù)首尾的兩個(gè)元素值；其取值為symmetric或periodic；默認(rèn)值為symmetric。圖3-2漢明窗及其頻譜特性注釋：漢明窗和漢寧窗都是余弦窗，只是加權(quán)系數(shù)不同，漢明窗旁瓣更小3布萊克曼窗函數(shù)的時(shí)域形式可以表示為w(k)k=1,2,N(3-5)k-1k-1=0.42-0.5cos2兀+0.08cos4兀N-1N-1它的頻域特性為W()=0.42Wr0.25Wr一系Wr.-2-O04WR士Wr士_.N-1,N-1J,N-1N-1(3-6)其中，Wr(附為矩形窗函數(shù)的幅度頻

12、率特性函數(shù)。布萊克曼窗函數(shù)的最大旁瓣值比主瓣值低57dB，但是主瓣寬度是矩形窗函數(shù)的主瓣寬度的3倍，為12N。Blackman函數(shù)：生成海明窗調(diào)用方式(1) w=blackman(n)：輸入?yún)?shù)n是窗函數(shù)的長度；輸出參數(shù)w是由窗函數(shù)的值組成的n階向量。(2) w=blackman(n,sflag):參數(shù)sflag用來控制窗函數(shù)首尾的兩個(gè)元素值；其取值為symmetric或periodic；默認(rèn)值為symmetric。圖3-3布萊克曼窗及其頻譜特性注釋：布萊克曼窗最大旁瓣值比主瓣值滴57db,但主瓣寬度是矩形窗的3倍3.2窗函數(shù)的指標(biāo)表3-1如圖所示對(duì)于任意一個(gè)窗函數(shù)，求出它的頻域值，并求出在主

13、瓣邊笫1個(gè)零點(diǎn)的位置；2,求出主瓣在-3dB處的位置；3,把笫1個(gè)零點(diǎn)白位置-(-3dB)處的位置，就是過渡帶的精確帶寬。矩形窗最簡單，但其-21dB的阻帶最小衰減在實(shí)際應(yīng)用中遠(yuǎn)遠(yuǎn)不夠。另外，矩形窗還會(huì)造成很強(qiáng)的吉布斯效應(yīng)。三角窗的阻帶衰減性能與矩形窗相比有所改善，但代價(jià)是過度帶加寬。窗函數(shù)主瓣寬度精確過渡帶最大旁瓣幅度最小阻帶衰減矩形窗4翼NL8WN-13dB-21dB三角窗S立N6/元/N-27dB-25dB漢寧窗8jt/N02支N-31dB-44dB漢明窗8nN6.6支N-41dB-53dB布萊克曼窗12MN11kN-5Z1B-73dB4應(yīng)用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的步驟4.1數(shù)字高

14、通濾波器的設(shè)計(jì):令高通濾波器的頻率響應(yīng)為e”竺0cW0<nHd(=W0,0<o<ochd(n)二工cej(n)d，工二ej(n)dd2二-二2二%1 -cj-ccos-(n-.)dsin(n-.)d'2 蹙2二一二1二jcos(n-)dsin(n-)d-2二L2二csink.(n-)J-sinIc(n一.)1二(n-)(4-1)(4-2)個(gè)低通濾從上述結(jié)果可以看出，一個(gè)高通濾波器相當(dāng)于用一個(gè)全通濾波器減去波器總結(jié)設(shè)計(jì)帶通濾波器時(shí)首先要計(jì)算出過渡帶，然后查表得到不同窗函數(shù)所需要的階數(shù)，不同的窗函數(shù)所設(shè)計(jì)的濾波器的形狀各有差異，尤其在主瓣寬度、旁瓣的形狀以及主瓣與旁瓣的高

15、度差上有比較明顯得差別，實(shí)際應(yīng)用中應(yīng)根據(jù)實(shí)際情況，折衷處理，兼顧各項(xiàng)指標(biāo)。為了這次課程設(shè)計(jì)，自己自學(xué)了數(shù)字信號(hào)處理領(lǐng)域中窗函數(shù)的有關(guān)知識(shí)。實(shí)際中遇到的離散時(shí)間信號(hào)總是有限長的，因此不可避免地要遇到數(shù)據(jù)截?cái)鄦栴}。而在信號(hào)處理中，對(duì)離散序列的數(shù)據(jù)截?cái)嗍峭ㄟ^序列與窗函數(shù)相乘來實(shí)現(xiàn)的。而且，有關(guān)濾波器的設(shè)計(jì)、功率譜估計(jì)等基本概念也要用到窗函數(shù)。本次課程設(shè)計(jì)對(duì)經(jīng)常用到的下面6窗函數(shù)：矩形窗函數(shù)、三角窗函數(shù)、漢寧窗函數(shù)、哈明窗函函數(shù)、布萊克曼窗函數(shù)、凱塞窗函窗，先是做了基本概念上的闡釋，然后對(duì)其MATLAB實(shí)現(xiàn)函數(shù)做出了說明，最后又結(jié)合具體的實(shí)例，對(duì)這些窗函數(shù)的頻域特性等進(jìn)行了介紹。通過這次學(xué)習(xí)，我不但掌

16、握了FIR數(shù)字濾波器窗函數(shù)的基本知識(shí)及其實(shí)際應(yīng)用的技巧了，還提高了自己的編程和寫報(bào)告的能力，收獲頗多。參考文獻(xiàn)1數(shù)字信號(hào)處理（第三版），丁玉美，高西全.西安電子科技大學(xué)出版社，2000.2 «MATLAB及在電子信息課程中的應(yīng)用，陳懷琛，吳大正，高西全.電子工業(yè)出版社，2006.3 «MATLAB7.0從入門到精通，求是科技.人民郵電出版社，2006.4數(shù)字信號(hào)處理（第三版）學(xué)習(xí)指導(dǎo)，高西全，丁玉美.西安科技大學(xué)出版社，2001.附錄FIR高通數(shù)字濾波器程序：clear;clc;closeall;S,Fs,Bit=wavreac('CUsersAdministrat

17、orDesktopFIR和IIRFIRSHE.wav');%讀取音頻信號(hào)sounc(S,Fs)；fnoisel=10000;%加入頻率為fnoise的正弦噪聲信號(hào)fnoise2=50;%加入頻率為fnoise的正弦噪聲信號(hào)N=length(S);T=N/Fs;t=T/N:T/N:T;forj=1:NNoise1(j)=sin(2*pi*fnoise1*t(j);Noise2(j)=sin(2*pi*fnoise2*t(j);end%對(duì)于低通濾波器：%51=S+0.1*Noise1'%fp=3000;%s=4000;%寸于高通濾波器：S1=S+Noise2'fp=3200

18、;fs=3000;%對(duì)于帶通濾波器：%S1=S+0.1*Noise1'+Noise2'%fp1=500;%fs1=100;%fp2=3000;%fs2=4000;sounc(S1,Fs);ap=1;as=100;Fs=44100;%wp=(fp*2*pi)/Fs5%ws=(fs*2*pi)/Fs5%DB=wp-ws;%Fs抽樣頻率關(guān)于冗的歸一化通帶截止頻率關(guān)于冗的歸一化阻帶截止頻率過渡帶寬N0=ceil(6.2*pi/DB);%計(jì)算所需h(n)長度N0,ceil(x)取大于等于x的最小整數(shù)N=N0+mod(N0+1,2);%確保h(n)長度N是奇數(shù)wc=(wp+ws)/2;%計(jì)

19、算理想高通濾波器通帶截止頻率(關(guān)于冗歸一化)h1=fir1(N-1,wc/pi,'high',hamming(N);%調(diào)用fir1計(jì)算高通FIRDF的h(n)H,f=freqz(h1,1,N,Fs);%求濾波器幅度響應(yīng)，設(shè)置最大幅度為1figure(1);plot(f,abs(H);figure(2);freqz(h1);%濾波器幅度響應(yīng)(db)和相位響應(yīng)title('hamming濾波器頻率響應(yīng)(幅頻上/相頻下)');%h2=fir1(N-1,wc/pi,hannin.N);subplot(2,1,1);stem(h1);%濾波器單位沖擊響應(yīng)序列及幅頻響應(yīng)F=abs(fft(h1);xlabel=0:Fs/(N-1):Fs/2;subplot(2,1,2);plot(xlabel,F(1:(N+1)/2);X1=abs(fft(S);Z1=(0:150000)/150000;subplot(3,1,1);plot(Z1,X1(1:150001);X2=abs(fft(S1);Z2=(0:150000)/150000;subplot(3,1,2);plot(Z2,X2(1:150001);S2=conv(S1,h1);sound(S2,Fs);X3=abs