基于窗函數法的FIR帶通濾波器的軟件實現(xiàn)

1、課程設計任務書20102011學年第一學期專業(yè)： 學號： 姓名： 課程設計名稱： 數字信號處理課程設計 設計題目： 基于窗函數法的FIR帶通濾波器的軟件實現(xiàn) 完成期限：自 年 月 日至 年 月 日共 周一設計目的1鞏固所學的理論知識。2提高綜合運用所學理論知識獨立分析和解決問題的能力。3更好地將理論與實踐相結合。4掌握信號分析與處理的基本方法與實現(xiàn)。5熟練使用MATLAB語言進行編程實現(xiàn)。二設計內容 編寫MATLAB程序實現(xiàn)FIR帶通濾波器的設計。指標如下：下阻帶截至頻率；上阻帶截止頻率；通帶下限頻率；通帶上限頻率；通帶最大衰減；阻帶最小衰減。三設計要求1、根據指標要求選擇合適的窗函數進行設計

2、；2、繪出及其幅頻響應并分析設計結果，驗證所設計的濾波器是否滿足要求。四設計條件計算機、MATLAB語言環(huán)境五、參考資料1數字信號處理（第三版），丁玉美，高西全.西安電子科技大學出版社，2000.2MATLAB及在電子信息課程中的應用，陳懷堔，吳大正，高西全.電子工業(yè)出版社，2006.3MATLAB 7.0從入門到精通，求是科技.人民郵電出版社，2006.4數字信號處理（第三版）學習指導，高西全，丁玉美.西安科技大學出版社，2001.指導教師（簽字）： 教研室主任（簽字）： 批準日期： 年 月 日摘 要無限長脈沖數字濾波器的設計方法只考慮了幅度特性，沒有考慮相位特性，所設的濾波器一般是某種確定

3、的非線性相位特性。有限脈沖響應（FIR）濾波器在保證了幅度特性滿足技術要求的同時，很容易做到有嚴格的線性相位特性。本課題利用MATLAB軟件實現(xiàn)。MATLAB是“矩陣實驗室”（MATrix LABoratoy）的縮寫，是一種科學計算軟件，它使用方便，輸入簡捷，運算高效，內容豐富，因此利用MATLAB軟件，通過一系列較為系統(tǒng)的函數法，根據已知的技術指標，就可以設計出滿足要求的濾波器。關鍵詞：MATLAB；窗函數；FIR帶阻數字濾波器；線性相位目 錄1課題描述12 MATLAB簡介13 窗函數設計法原理23.1線性相位43.2 基本窗函數53.2.1 矩形窗函數53.2.2 三角窗函數53.2.3

4、漢寧窗函數63.2.4 哈明窗函數73.2.5 布萊克曼窗函數73.2.6 凱塞窗84 方案設計與分析94.1方案設計程序94.2分析13總結145參考文獻151課題描述現(xiàn)代圖像、語聲、數據通信對線性相位的要求是普遍的。正是此原因，使得具有線性相位的FIR數字濾波器得到大力發(fā)展和廣泛應用。在實際進行數字信號處理時，往往需要把信號的觀察時間限制在一定的時間間隔內，只需要選擇一段時間信號對其進行分析。這樣，取用有限個數據，即將信號數據截斷的過程，就等于將信號進行加窗函數操作。而這樣操作以后，常常會發(fā)生頻譜分量從其正常頻譜擴展開來的現(xiàn)象，即所謂的“頻譜泄漏”。當進行離散傅立葉變換時，時域中的截斷是必

5、需的，因此泄漏效應也是離散傅立葉變換所固有的，必須進行抑制。而要對頻譜泄漏進行抑制，可以通過窗函數加權抑制DFT的等效濾波器的振幅特性的副瓣，或用窗函數加權使有限長度的輸入信號周期延拓后在邊界上盡量減少不連續(xù)程度的方法實現(xiàn)。而在后面的FIR濾波器的設計中，為獲得有限長單位取樣響應，需要用窗函數截斷無限長單位取樣響應序列。另外，在功率譜估計中也要遇到窗函數加權問題。由此可見，窗函數加權技術在數字信號處理中的重要地位。2 MATLAB簡介MATLAB是“矩陣實驗室”（MATrix LABoratoy）的縮寫，是一種科學計算軟件，主要適用于矩陣運算及控制和信息處理領域的分析設計，它使用方便，輸入簡捷

6、，運算高效，內容豐富，因此很多專家在自己擅長的領域用它編寫了許多專門的MATLAB工具包，由于MATLAB功能的不斷擴展，所以是科學研究中最常用必不可少的工具。MATLAB由一系列工具組成。這些工具方便用戶使用MATLAB的函數和文件，其中許多工具采用的是圖形用戶界面。包括MATLAB桌面和命令窗口、歷史命令窗口、編輯器和調試器、路徑搜索和用于用戶瀏覽幫助、工作空間、文件的瀏覽器。隨著MATLAB的商業(yè)化以及軟件本身的不斷升級，MATLAB的用戶界面也越來越精致，更加接近Windows的標準界面，人機交互性更強，操作更簡單。而且新版本的MATLAB提供了完整的聯(lián)機查詢、幫助系統(tǒng)，極大的方便了用

7、戶的使用。簡單的編程環(huán)境提供了比較完備的調試系統(tǒng)，程序不必經過編譯就可以直接運行，而且能夠及時地報告出現(xiàn)的錯誤及進行出錯原因分析。MATLAB一個高級的矩陣/陣列語言，它包含控制語句、函數、數據結構、輸入和輸出和面向對象編程特點。用戶可以在命令窗口中將輸入語句與執(zhí)行命令同步，也可以先編寫好一個較大的復雜的應用程序（M文件）后再一起運行。新版本的MATLAB語言是基于最為流行的C語言基礎上的，因此語法特征與C語言極為相似，而且更加簡單，更加符合科技人員對數學表達式的書寫格式。使之更利于非計算機專業(yè)的科技人員使用。而且這種語言可移植性好、可拓展性極強，這也是MATLAB能夠深入到科學研究及工程計算

8、各個領域的重要原因。3 窗函數設計法原理數字濾波器可以理解為是一個計算程序或算法，將代表輸入信號的數字時間序列轉化為代表輸出信號的數字時間序列，并在轉化過程中，使信號按預定的形式變化。數字濾波器有多種分類，根據數字濾波器沖激響應的時域特征，可將數字濾波器分為兩種，即無限長沖激響應（iir）濾波器和有限長沖激響應（fir）濾波器。iir數字濾波器具有無限寬的沖激響應，與模擬濾波器相匹配。所以iir濾波器的設計可以采取在模擬濾波器設計的基礎上進一步變換的方法。fir數字濾波器的單位脈沖響應是有限長序列。它的設計問題實質上是確定能滿足所要求的轉移序列或脈沖響應的常數問題，設計方法主要有窗函數法、頻率

9、采樣法和等波紋最佳逼近法等。因此設計FIR濾波器的方法之一可以從時域出發(fā)，截取有限長的一段沖擊響應作為H(z)的系數，沖擊響應長度N就是系統(tǒng)函數H(z)的階數。只要N足夠長，截取的方法合理，總能滿足頻域的要求。一般這種時域設計、頻域檢驗的方法要反復幾個回合才能成功。要設計一個線性相位的FIR數字濾波器，首先要求理想頻率響應。是w的周期函數，周期為，可以展開成傅氏級數： （3-1）其中是與理想頻響對應的理想單位抽樣響應序列。但不能用來作為設計FIR DF用的h(n)，因為一般都是無限長、非因果的，物理上無法實現(xiàn)。為了設計出頻響類似于理想頻響的濾波器，可以考慮用h(n)來近似。窗函數的基本思想：先

10、選取一個理想濾波器（它的單位抽樣響應是非因果、無限長的），再截?。ɑ蚣哟埃┧膯挝怀闃禹憫玫骄€性相位因果FIR濾波器。這種方法的重點是選擇一個合適的窗函數和理想濾波器。設x(n)是一個長序列，w(n)是長度為N的窗函數，用w(n)截斷x(n)，得到N點序列xn(n)，即 xn(n) = x(n) w(n) （3-2）在頻域上則有 （3-3）由此可見，窗函數w(n)不僅僅會影響原信號x(n)在時域上的波形，而且也會影響到頻域內的形狀。MATLAB信號工具箱主要提供了以下幾種窗函數，如表下所示：表1 MATLAB窗函數窗窗 函 數矩形窗Boxcar三角窗Barlett布萊克曼窗Blackman哈

11、明窗Hamming漢寧窗Hanning凱塞窗Kaiser切比雪夫窗Chebwin加矩形窗后的頻譜和理想頻譜可得到以下結論：加窗使過渡帶變寬，過渡帶的帶寬取決于窗譜的主瓣寬度。矩形窗情況下的過渡帶寬是。N越大，過渡帶越窄、越陡；過渡帶兩旁產生肩峰，肩峰的兩側形成起伏振蕩。肩峰幅度取決于窗譜主瓣和旁瓣面積之比。矩形窗情況下是8.95，與N無關。工程上習慣用相對衰耗來描述濾波器，相對衰耗定義為 （3-4）這樣兩個肩峰點的相對衰耗分別是0.74dB和-21dB。其中（-0.0895）對應的點的值定義為阻帶最小衰耗。以上的分析可見，濾波器的各種重要指標都是由窗函數決定，因此改進濾波器的關鍵在于改進窗函數

12、。窗函數譜的兩個最重要的指標是：主瓣寬度和旁瓣峰值衰耗。旁瓣峰值衰耗定義為： 旁瓣峰值衰耗20lg(第一旁瓣峰值主瓣峰值) （3-5）為了改善濾波器的性能，需使窗函數譜滿足：主瓣盡可能窄，以使設計出來的濾波器有較陡的過渡帶；第一副瓣面積相對主瓣面積盡可能小，即能量盡可能集中在主瓣，外泄少，使設計出來的濾波器的肩峰和余振小。但上面兩個條件是相互矛盾的，實際應用中，折衷處理，兼顧各項指標。3.1線性相位一個單一頻率的正弦信號通過一個系統(tǒng)，假設它通過這個系統(tǒng)的時間需要t，則這個信號的輸出相位落后原來信號wt的相位。從這邊可以看出，一個正弦信號通過一個系統(tǒng)落后的相位等于它的w*t；反過來說，如果一個頻

13、率為w的正弦信號通過系統(tǒng)后，它的相位落后delta，則該信號被延遲了delta/w的時間。在實際系統(tǒng)中，一個輸入信號可以分解為多個正弦信號的疊加，為了使得輸出信號不會產生相位失真，必須要求它所包含的這些正弦信號通過系統(tǒng)的時間是一樣的。因此每一個正弦信號的相位分別落后，w1*t，w2*t，w3*t。因此，落后的相位正比于頻率w，如果超前，超前相位的大小也是正比于頻率w。從系統(tǒng)的頻率響應來看，就是要求它的相頻特性是一條直線。在FIR濾波器的設計中，為了得到線性相位的性質，通常利用實偶對稱序列的相頻特性為常數0和實奇對稱序列為相頻特性為常數90度的特點。因此得到的是對稱序列，不是因果序列，是不可實現(xiàn)

14、系統(tǒng)，為了稱為物理可實現(xiàn)系統(tǒng)，需要將它向右移動半個周期，這就造成了相移特性隨時間的變化，同時也是線性變化。單位脈沖響應h(n)（為實數）具有偶對稱或奇對稱性，則FIR數字濾波器具有嚴格的線性相位特性。 數字濾波器中，IIR數字濾波器方便簡單，但它相位的非線性，要求采用全通網絡進行相位校正，且穩(wěn)定性難以保障。FIR濾波器具有很好的線性相位特性，使得它越來越受到廣泛的重視。3.2 基本窗函數數字信號處理領域中所用到的基本窗函數主要有：矩形窗函數、三角窗函數和漢寧窗函數，哈明窗函數布萊克窗函數，凱塞窗函數等。下面就對這些窗函數展開介紹。3.2.1 矩形窗函數矩形窗(Rectangular Windo

15、w)函數的時域形式可以表示為： （3-2-1）它的頻域特性為 （3-2-2）Boxcar函數：生成矩形窗調用方式w = boxcar (N)：輸入參數N是窗函數的長度；輸出參數w是由窗函數的值組成的n階向量。從功能上講，該函數又等價于w = ones(n,1)。3.2.2 三角窗函數三角窗是最簡單的頻譜函數為非負的一種窗函數。三角窗函數的時域形式可以表示為：當n為奇數時 （3-2-3）當n為偶數時 （3-2-4）Bartlett函數：生成巴特利特窗調用方式w = bartlett(n)：(1) 輸入參數n是窗函數的長度；(2) 輸出參數w是由窗函數的值組成的n階向量。(3) 三角窗也是兩個矩形

16、窗的卷積。三角窗函數的首尾兩個數值通常是不為零的。當n是偶數時，三角窗的傅立葉變換總是非負數。3.2.3漢寧窗函數漢寧窗函數的時域形式可以表示為： （3-2-5）它的頻域特性為： （3-2-6）其中，為矩形窗函數的幅度頻率特性函數。漢寧窗函數的最大旁瓣值比主瓣值低31dB，但是主瓣寬度比矩形窗函數的主瓣寬度增加了1倍，為8/N。hanning函數：生成漢寧窗調用方式(1) w = hanning(n)：輸入參數n是窗函數的長度；輸出參數w是由窗函數的值組成的n階向量。注意：此函數不返回是零點的窗函數的首尾兩個元素。 (2) w = hanning(n,'symmetric')：

17、與上面相類似。(3) w = hanning(n,'periodic')：此函數返回包括為零點的窗函數的首尾兩個元素。3.2.4 哈明窗函數海明窗函數的時域形式可以表示為 （3-2-7）它的頻域特性為 （3-2-8）其中，為矩形窗函數的幅度頻率特性函數。海明窗函數的最大旁瓣值比主瓣值低41dB，但它和漢寧窗函數的主瓣寬度是一樣大的。Hamming函數：生成海明窗調用方式(1) w = hamming(n)：輸入參數n是窗函數的長度；輸出參數w是由窗函數的值組成的n階向量。(2) w = hamming(n,sflag)：參數sflag用來控制窗函數首尾的兩個元素值；其取值為sy

18、mmetric或periodic；默認值為symmetric。3.2.5 布萊克曼窗函數布萊克曼窗函數的時域形式可以表示為 (3-2-9)它的頻域特性為 （3-2-10）其中，為矩形窗函數的幅度頻率特性函數。布萊克曼窗函數的最大旁瓣值比主瓣值低57dB，但是主瓣寬度是矩形窗函數的主瓣寬度的3倍，為12/N。Blackman函數：生成海明窗調用方式(1) w = blackman (n)：輸入參數n是窗函數的長度；輸出參數w是由窗函數的值組成的n階向量。(2) w = blackman (n,sflag)：參數sflag用來控制窗函數首尾的兩個元素值；其取值為symmetric或periodic

19、；默認值為symmetric。3.2.6 凱塞窗上面所討論的幾種窗函數，在獲得旁瓣抑制的同時卻增加了主瓣的寬度。而凱塞窗定義了一組可調的窗函數，它是由零階貝塞爾函數構成的，其主瓣能量和旁瓣能量的比例是近乎最大的。而且，這種窗函數可以在主瓣寬度和旁瓣高度之間自由選擇它們的比重，使用戶的設計變得非常靈活。凱塞窗函數的時域形式可表示為 (3-2-11)其中，是第1類變形零階貝塞爾函數，是窗函數的形狀參數，由下式確定： （3-2-12）其中，為凱塞窗函數的主瓣值和旁瓣值之間的差值(dB)。改變的取值，可以對主瓣寬度和旁瓣衰減進行自由選擇。的值越大，窗函數頻譜的旁瓣值就越小，而其主瓣寬度就越寬。Kais

20、er函數：生成凱塞窗調用方式w = kaiser(n,beta)：輸入參數n是窗函數的長度；輸入參數beta用于控制旁瓣的高度；輸出參數w是由窗函數的值組成的n階向量。n一定時，beta越大，其頻譜的旁瓣就越小，但主瓣寬度相應的增加；當beta一定時，n發(fā)生變化，其旁瓣高度不會發(fā)生變化。4 方案設計與分析用窗函數法設計一個FIR帶阻濾波器。指示如下：下阻帶截止頻率 ；上阻帶截止頻率 ；通帶下限頻率 ；通帶上限頻率 通帶最大衰減 阻帶最小衰減 =4.1方案設計程序6中窗函數的基本參數如下窗函數類型旁瓣峰值/dB過渡帶寬度阻帶最小衰減/dB近似值精確值矩形窗-134/N1.8/N-12三角窗-25

21、8/N6.1/N-25漢寧窗-318/N6.2/N-44哈明窗-418/N6.6/N-53布萊克曼窗-5712/N11/N-74凱塞窗-5710/N-80因為阻帶最小衰減=,所以選擇布萊克曼窗或凱塞窗都可以設計，先以凱塞窗進行設計。程序步驟如下：wls = 0.2*pi; wlp = 0.35*pi; whp = 0.65*pi; wc = wlp/pi,whp/pi; %計算理想帶通濾波器通帶截止頻率 B = wlp-wls; N = ceil(7.24/0.15); n=0:N-1; window=kaiser(N); h1,w=freqz(window,1); subplot(2,2,1

22、) stem(window,'.'); xlabel('n'); title(' kaiser窗函數'); subplot(2,2,2) plot(w/pi,20*log(abs(h1)/abs(h1(1); %繪制損耗函數 grid; xlabel('w/pi'); ylabel('幅度(dB)'); title(' kaiser窗函數的頻譜'); hn = fir1(N-1,wc, kaiser (N); %fir1基于窗函數FIR數字濾波器的設計h2,w=freqz(hn,1,512); %f

23、reqz數字濾波器的頻率響應 subplot(2,2,3) stem(n,hn,'.'); xlabel('n'); ylabel('h(n)'); title(' kaiser窗函數的單位脈沖響應'); subplot(2,2,4) plot(w/pi,20*log(abs(h2)/abs(h2(1); %繪制損耗函數 grid; xlabel('w/pi'); ylabel('幅度(dB)'); title(' kaiser帶通濾波器的幅度特性');運行圖形如下：圖4.1.1bl

24、ackman窗函數程序清單如下：wls = 0.2*pi; wlp = 0.35*pi; whp = 0.65*pi; wc = wlp/pi,whp/pi; %計算理想帶通濾波器通帶截止頻率 B = wlp-wls; N = ceil(12/0.15); n=0:N-1; window=blackman(N); h1,w=freqz(window,1); subplot(2,2,1) stem(window,'.'); xlabel('n'); title('blackman窗函數'); subplot(2,2,2) plot(w/pi,20*

25、log(abs(h1)/abs(h1(1); %繪制損耗函數grid; xlabel('w/pi'); ylabel('幅度(dB)'); title(' blackman窗函數的頻譜'); hn = fir1(N-1,wc, kaiser (N); %fir1基于窗函數FIR數字濾波器的設計 h2,w=freqz(hn,1,512); %freqz數字濾波器的頻率響應 subplot(2,2,3) stem(n,hn,'.'); xlabel('n'); ylabel('h(n)'); title(' blackman窗函數的單位脈沖響應'); subplot(2,2,4) plot(w/pi,20*log(abs(h2)/abs(h2(1); %繪制損耗函數grid; xlabel('w/pi'); ylabel('幅度(d