1、摘要 本次課程設計的題目是基于語音信號去噪處理的FIR濾波器設計，主要是利用MATLAB設計一個FIR數字帶通濾波器，對一段噪聲環境下的語音信號進行濾波。在設計過程中，首先要錄制一段語音信號，并對錄制好的信號進行時域和頻域分析；然后對原始的語音信號進行加噪處理，對加噪后的信號進行分析，繪出時域和頻域的分析圖；最后利用設計出的FIR帶通濾波器，針對語音信號的性質選取一種適合的窗函數設計濾波器進行濾波，最后對仿真結果進行分析。關鍵詞: FIR帶通濾波器；語音信號；MATLAB仿真；加噪；濾波目錄前言1第一章 基本原理21.1 語音處理中的采樣原理21.2數字濾波器的設計21.3窗函數法4第二章 F

2、IR帶通濾波器設計62.1 FIR濾波器簡介62.2 FIR帶通濾波器設計要求72.3 設計方法72.4 設計步驟9第三章 FIR帶通濾波器的軟件仿真103.1 程序流程圖103.2 仿真結果及分析11總結15參考文獻16附錄17致謝22一、 前言在信號處理過程中，所處理的信號往往混有噪音，從接收到的信號中消除或減弱噪音是信號傳輸和處理中十分重要的問題。根據有用信號和噪音信號的不同特性，提取有用信號的過程稱為濾波。實現濾波功能的系統被稱為濾波器。在近代電信設備和各類控制系統中，數字濾波器應用極為廣泛。語音處理是最早應用數字濾波器的領域之一，也是最早推動數字信號處理理論發展的領域之一。該領域主要

3、包括5個方面的內容：第一，語音信號分析。即對語音信號的波形特性，統計特性，模型參數等進行分析計算。第二，語音合成。即利用專用硬件或在通用計算機上運行軟件來產生語音。第三，語音識別。即利用專用硬件或計算機識別人的講話，或者識別說話的人。第四，語音增強。即從噪音或者干擾中提取被掩蓋的語音信號。第五，語音編碼。主要用于語音數據的壓縮，目前已經建立了一系列語音編碼國際標準，大量用于通信和音頻處理。MATLAB是一種面向科學和工程計算的語言，它集數值分析、矩陣運算、信號處理和圖形顯示于一體，具有編程效率高、調試手段豐富、擴充能力強等特點。MATLAB的信號處理工具箱具有強大的函數功能，它不僅可以用來設計

4、數字濾波器，還可以使設計達到最優化，是數字濾波器設計的強有力工具。第一章 基本原理1.1 語音處理中的采樣原理因為錄制的語音信號是模擬信號，要想使用數字濾波器對疊加了噪聲的信號進行濾波，則在設計數字濾波器之前首先要進行模數轉換，將模擬信號轉換為數字信號。在進行模數轉換的過程中，當最高采樣頻率fs大于信號中最高頻率f的2倍時，即：fsmax2fmax,采樣之后的數字信號可以完整地保留原始信號中的信息，一般實際應用中保證采樣頻率為信號最高頻率的510倍；采樣定理又稱奈奎斯特定理。頻帶為F的連續信號 f(t)可用一系列離散的采樣值f(t1),f(t1±t)，f(t1±2t)，.來

5、表示,只要這些采樣點的時間間隔t1/2F，便可根據各采樣值完全恢復原來的信號f(t)。這是時域采樣定理的一種表述方式。時域采樣定理的另一種表述方式是：當時間信號函數f(t)的最高頻率分量為fm時,f(t)的值可由一系列采樣間隔小于或等于1/2fm的采樣值來確定,即采樣點的重復頻率f2fm。 采樣頻率，也稱為采樣速度或者采樣率，定義了每秒從連續信號中提取并組成離散信號的采樣個數，它用赫茲（Hz）來表示。采樣頻率的倒數是采樣周期或者叫采樣時間，它是采樣之間的時間間隔。通俗的講采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質、音調，衡量聲卡、聲音文件的質量標準。采樣位數和采樣率對于音

6、頻接口來說是最為重要的兩個指標，也是選擇音頻接口的兩個重要標準。無論采樣頻率如何，理論上來說采樣的位數決定了音頻數據最大的力度范圍。每增加一個采樣位數相當于力度范圍增加了6dB。采樣位數越多則捕捉到的信號越精確1。1.2數字濾波器的設計 數字濾波器(Digital Filter，簡稱為DF)是指用來對輸入信號進行濾波的硬件和軟件。所謂數字濾波器，是指輸入、輸出均為數字信號，通過一定運算關系改變輸入信號所含頻率成分的相對比例或者濾除某些頻率成分的器件。數字濾波器和模擬濾波器相比，因為信號的形式和實現濾波的方式不同，數字濾波器具有比模擬濾波器精度高、穩定、體積小、重量輕、靈活、不要求阻抗

7、匹配等優點。一般用兩種方法來實現數字濾波器:一是采用通用計算機，把濾波器所要完成的運算編成程序通過計算機來執行，也就是采用計算機軟件來實現;二是采用實際專用的數字處理硬件。數字濾波器根據其沖擊響應函數的時域特性，可分為兩種，即無限沖擊響應IIR濾波器和有限沖擊響應FIR濾波器。IIR濾波器的特性是，具有無限持續時間沖擊響應。這種濾波器一般需要用遞歸模型來實現，因而有時也稱之為遞歸濾波器。FIR濾波器的沖擊響應只能延續一定的時間，在工程實際中可以采用遞歸的方法來實現，也可以采用非遞歸的方式來實現。數字濾波器的設計方法有多種，如雙線性變換法，窗函數法，插值逼近法和切比雪夫逼近法等等。隨著MATLA

8、B軟件尤其是MATLAB信號處理工作的不斷完善，不僅數字濾波器的計算機輔助設計有了可能，而且還可以使設計達到最優化2。數字濾波器的設計基本步驟如下：（1）確定指標在設計一個濾波器之前，必須首先根據工程實際的需要確定濾波器的技術指標。在很多實際應用中，數字濾波器常常被用來實現選頻操作。因此，指標的形式一般在頻域中給出幅度和相位響應。幅度指標主要用以下兩種方式給出。第一種是絕對指標。它提供對幅度響應函數的要求，一般應用于FIR濾波器的設計。第二種指標是相對指標。它以分貝值的形式給出要求。運用線性相位響應指標進行濾波器設計具有如下優勢：只包含實數算法，不涉及復數運算。不存在延時失真，只有固定數量的延

9、遲。長度為N的濾波器，計算量為N/2數量級。（2）逼近確定了技術指標后，就可以建立一個目標的數字濾波器模型。通常采用理想的數字濾波器模型。之后，利用數字濾波器的設計方法，設計出一個實際濾波器模型來逼近給定的目標。（3）性能分析和計算機仿真上兩步的結果是得到以差分或者系統函數或者沖擊響應描述的濾波器。根據這個描述就可以分析器頻率特性和相位特性，以驗證設計結果是否滿足指標要求：或者利用計算機仿真實現設計的濾波器，再分析濾波器結果來判斷3。1.3窗函數法數字濾波器可以理解為是一個計算程序或算法，將代表輸入信號的數字時間序列轉化為代表輸出信號的數字時間序列，并在轉化過程中，使信號按預定的形式變化。數字

10、濾波器有多種分類，根據數字濾波器沖激響應的時域特征，可將數字濾波器分為兩種，即無限長沖激響應（IIR）濾波器和有限長沖激響應（FIR）濾波器。FIR數字濾波器的單位脈沖響應是有限長序列。它的設計問題實質上是確定能滿足所要求的轉移序列或脈沖響應的常數問題，設計方法主要有窗函數法、頻率采樣法和等波紋最佳逼近法等2。用窗函數設計濾波器首先要對濾波器提出性能指標。一般是給定一個理想的頻率響應，使所設計的FIR濾波器的頻率響應去逼近所要求的理想的濾波器的響應。窗函數法設計的任務在于尋找一個可實現（有限長單位脈沖響應）的傳遞函數。用一定寬度窗函數截取無限脈沖響應序列獲得有限長的脈沖響應序列4。主要設計步驟

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

12、取一個理想濾波器（它的單位抽樣響應是非因果、無限長的），再截取（或加窗）它的單位抽樣響應得到線性相位因果FIR濾波器，這種方法的重點是選擇一個合適的窗函數和理想濾波器6。（2）由性能指標確定窗函數W(n)和窗口長度N。設x(n)是一個長序列，w(n)是長度為N的窗函數，用w(n)截斷x(n)，得到N點序列xn(n)，即 （1-2）在頻域上則有 （1-3） (3) 求得實際濾波器的單位脈沖響應h(n)， h(n)即為所設計FIR濾波器系 數向量。 （1-4）由此可見，窗函數w(n)不僅僅會影響原信號x(n)在時域上的波形，而且也會影響到頻域內的形狀4。第二章 FIR帶通濾波器設計2.1 FIR濾

13、波器簡介數字濾波器包括FIR（有限單位脈沖響應）濾波器與IIR（無限單位脈沖響應）濾波器兩種。在現代信號處理技術中，例如數據傳輸、雷達接收以及一些要求較高的電子系統，都越來越多地要求信道具有線性的相位特性。在這方面，FIR濾波器具有獨到的優點，它可以在幅度特性隨意設計的同時，保證精確、嚴格的線性相位特性5。FIR濾波器的單位脈沖響應h（n）是有限長的（0nN-1），其z變換為的（N-1）階多項式： （2-1）可得FIR濾波器的系統差分方程為： （2-2）因此，FIR濾波器又稱為卷積濾波器。FIR濾波器的頻率響應表達式為： （2-3）信號通過FIR濾波器不失真條件是在通帶內具有恒定的幅頻特性和線

14、性相位特性。理論上可以證明：當FIR濾波器的系數滿足下列中心對稱條件： （2-4）或者 （2-5）時，濾波器設計在逼近平直幅頻特性的同時，還能獲得嚴格的線性相位特性。線性相位FIR濾波器的相位滯后和群延遲在整個頻帶上是相等且不變的。對于一個 N 階的線性相位FIR濾波器，群延遲為常數，即濾波后的信號簡單地延遲常數個時間步長。這一特性使通帶頻率內信號通過濾波器后仍保持原有波形形狀而無相位失真7。FIR濾波器的設計任務是選擇有限長度的h(n)，使傳輸函數滿足技術要求。FIR濾波器的設計方法有多種，如窗函數法、頻率采樣法及其它各種優化設計方法，本次設計使用窗函數法設計FIR帶通濾波器。2.2 FIR

15、帶通濾波器設計要求利用MATLAB仿真軟件系統結合窗函數法設計一個數字帶通FIR濾波器。設計要求：（1）用所設計的濾波器對受噪聲影響的信號進行濾波，畫出濾波后語音信號的時域波形圖和頻譜圖；（2）對濾波前后的信號進行對比，分析信號的變化；回放語音信號，并與原始語音信號對比。技術指標：低端阻帶截止頻率 fc11000 Hz 低端通帶截止頻率 fb11200 Hz 高端通帶截止頻率 fb23000 Hz高端阻帶截止頻率 fc23200 Hz通帶衰減系數 ap1dB 阻帶衰減系數 as100 dB2.3 設計方法設計FIR數字濾波器的方法通常有三種：窗函數法，頻率抽樣法，等紋波逼近法。本次課程設計討論

16、的是第一種窗函數法。這種方法也叫傅里葉級數法。一般是先給定所要求的理想濾波器頻率響應，導出，我們知道理想濾波器的沖擊響應是無限長的非因果序列，而我們要設計的是h(n)是有限長的FIR濾波器，所以要用有限長序列h(n)來逼近無限長序列，設： （2-6）常用的方法是使用有限長的窗函數w(n)來截取即： （2-7）根據在時域是相乘關系，在頻域則是卷積關系： （2-8）其中為矩形窗譜，是FIR濾波器頻率響應8。在設計過程中，將無限長序列變為有限長序列是通過時域加矩形窗乘積來實現的截斷。常見的窗函數有：矩形窗，漢寧窗，海明窗，布拉克曼窗，凱澤窗等，下面分別來討論著幾種常見的窗函數：矩形窗：窗譜： 幅度函

17、數： 海明窗：布拉克曼窗：凱澤窗：其中(x)是第一類變形修正零階貝塞爾函數 這些窗函數的基本參數如表2-1所示：表2-1 窗函數的基本參數名稱旁瓣峰值/dB近似過渡帶寬精確過渡帶寬最小阻帶衰減/dB矩形窗-134/N1.8/N-21海明窗-418/N6.6/N-53布拉克曼窗-5712/N11/N-74凱澤窗(=7.865)-5710/N-802.4 設計步驟 窗函數法設計FIR濾波器的主要步驟如下： a.給出希望設計的濾波器的頻率響應的函數； b.根據允許的過渡帶寬度及阻帶衰減，初步選定窗函數及其長度N； c.根據技術要求確定待求濾波器的單位取樣響應， d.將與窗函數相乘得FIR數字濾波器的

18、單位取樣響應， e.按如下方法計算FIR數字濾波器的頻率響應，并驗證是否達到所要求的技術指標： （2-9） 或 由計算幅度響應和相位響應。計算式（15）時可用FFT算法。如果或不滿足要求，可根據具體情況重復b,c,d,e步驟，直到滿足技術要求9。第三章 FIR帶通濾波器的軟件仿真3.1 程序流程圖開始在Window下錄制語音將語音格式改為wav對語音信號進行頻譜分析，畫出時域和頻域波形圖加入噪聲畫出其頻率響應用FIR濾波器對語音信號進行濾波畫出語音信號濾波前后的波形并進行比較結束用Blackman窗設計FIR帶通濾波器圖3.1 程序設計流程圖首先要錄制一段語音信號，并對錄制好的信號進行時域和頻

19、域分析；然后對原始的語音信號進行加噪處理，對加噪后的信號進行分析，繪出時域和頻域的分析圖；最后利用設計出的FIR帶通濾波器，針對語音信號的性質選取一種適合的窗函數設計濾波器進行濾波，最后對仿真結果進行分析。3.2 仿真結果及分析對錄制好的語音信號進行分析，編寫程序對其進行仿真，原始信號時域和頻域仿真圖如圖3.2所示。（程序見附錄）圖3.2 原始信號的時域與頻域分析 對原始信號加入8.5khz的余弦噪聲，并將噪聲信號及其加噪后的信號進行時域和頻域的分析，仿真結果如圖3.3所示。（程序見附錄） 圖3.3 噪聲信號及加噪信號的時域與頻域仿真圖仿真分析：通過對圖3.2和3.3兩張圖進行對比，可以看出加

20、噪后的語音信號時域波形比原始語音信號渾濁了許多，加噪后的信號的頻帶寬度也明顯比原始語音信號的窄；再通過對原始語音信號的回放效果與加噪后的語音信號回放的效果的對比，人耳可以明顯辨別出兩種語音信號不一樣了，加噪后的語音信號在聽覺上比原始語音信號要渾濁很多，而且還有混雜音。 利用給定的參數設計FIR帶通濾波器，對其在Matlab上進行仿真，結果如圖3.4所示。（程序見附錄）圖3.4 帶通濾波器的幅頻與相頻特性 將加噪后的信號通過設計出的帶通濾波器進行濾波，時域和頻域的仿真結果如圖3.5所示。（程序見附錄）圖3.5 加噪信號經過濾波后的時域與頻域仿真仿真分析：分析濾波前后信號的變化，波形圖變窄。頻譜變

21、化表明:信號的低頻率段和高頻率段被濾除，中間頻率段被保留；分析濾波前后的聲音變化，濾波后聲音一定程度上變的尖銳，與高通濾波器濾波后的聲音相比較低。這說明濾波器設計基本符合指標要求。總結 在信號處理過程中，所處理的信號往往混有噪音，從接收到的信號中消除或減弱噪音是信號傳輸和處理中十分重要的問題。根據有用信號和噪音信號的不同特性，提取有用信號的過程稱為濾波。本次課程設計主要是利用MATLAB設計一個FIR數字帶通濾波器，對一段錄制好的語音信號加入噪聲，在噪聲環境下對語音信號進行濾波。 設計主要分為以下幾個步驟：首先要錄制一段語音信號，對錄制好的信號進行時域和頻域分析；然后對原始的語音信號進行加噪處

22、理，對加噪后的信號進行分析，繪出時域和頻域的分析圖；最后利用窗函數法設計出的FIR帶通濾波器，對加噪后的信號進行濾波，對比濾波前和濾波后的時域和頻域的仿真圖對仿真結果進行分析。 經過兩周的時間，在老師和同學的幫助下雖然完成了設計任務，但是還存在許多不足的地方，比如濾波后的語音信號還是有少許的雜音，沒有做更多的濾波效果比較等。在以后的工作和學習中會更加努力來完善設計任務。 參考文獻1 丁玉美. 數字信號處理M. 西安：西安電子科技大學出版社，2003，3.2 朱冰蓮. 數字信號處理M. 北京：電子工業出版社，2003，7. 3 程佩青.數字信號處理教程（第3版）.北京：清華大學出版社，2007.

23、4 陳亞勇等.MATLAB信號處理詳解.北京：人民郵電出版社，2001.5 萬永革.數字信號處理的MATLAB實現.北京：科學出版社，2007.6 王力寧.MATLAB與通信仿真.北京：人民郵電出版社，1999.7 宋壽鵬.數字濾波器設計及工程應用.江蘇：江蘇大學出版社,2007.8 普埃克著,方艷梅譯.數字信號處理(第四版).北京：電子工業出版社,2007.9 劉波MATLAB信號處理M北京：電子工業出版社，2006.附錄%原始語音信號的時域與頻域分析y,fs,nbits=wavread ('H:信號處理課設新建文件夾標準狼叫.wav');y=y(:,1); N=length

24、 (y); %求出語音信號的長度Y=fft(y,N); %傅里葉變換subplot(2,1,1);t=(0:N-1)/fs; figure(1);plot(t,y);axis(0 12 -1.5 1.5);title('原始信號波形');subplot(2,1,2);plot(abs(Y);title('原始信號頻譜');%對噪聲信號及加噪信號進行時域與頻域的譜分析f=fs*(0:1023)/2048;Au=0.5;d=Au.*cos(2*pi*8500*t)' %噪聲為8.5kHz的余弦信號d=d(:,1); dd=fft(d,N); %傅里葉變換fi

25、gure(2);subplot(2,1,1);plot(t,d);axis(0 12 -0.6 0.6)title('噪聲信號波形');subplot(2,1,2);plot(abs(dd);title('噪聲信號頻譜');x2=y+d;S=fft(x2); figure(3);subplot(2,1,1);plot(t,x2);axis(0 12 -2 2)title('加噪信號波形');subplot(2,1,2);plot(abs(S);title('加噪信號頻譜');%帶通濾波器的設計fp1=1200;fp2=3000;f

26、s1=1000;fs2=3200;Fs=8000;%輸入設計指標wp1=fp1/(Fs/2);%計算歸一化角頻率wp2=fp2/(Fs/2);ws1=fs1/(Fs/2);ws2=fs2/(Fs/2);deltaw=wp1-ws1;%計算過渡帶的寬度N0=ceil(11/deltaw);%按布萊克曼窗設計算濾波器長度N=N0+mod(N0+1,2)%為實現FIR類型I偶對稱濾波器，應確保N為奇數windows=blackman(N);%使用布萊克曼窗wc1=(ws1+wp1)/2;%截止頻率取通阻帶頻率的平均值wc2=(ws2+wp2)/2;b=fir1(N-1,wc1,wc2,windows

27、);%用fir1子函數求系統函數系數db,mag,pha,grd,w=freqz_m(b,1);n=0:N-1;dw=2/1000;Rp=-(min(db(wp1/dw+1:wp2/dw+1)ws0=1:ws1/dw+1,ws2/dw+1:501;As=-round(max(db(ws0)figure(4);subplot(2,2,1)stem(0:N-1,b);axis(0,N,1.1*min(b),1.1*max(b);title('實際脈沖響應');xlabel('n');ylabel('h(n)');subplot(2,2,2)stem(

28、n,windows);axis(0,N,0,1.1);title('窗函數特性');xlabel('n');ylabel('wh(n)');subplot(2,2,3)plot(w/pi,db);axis(0,1,-100,10);title('幅頻響應');xlabel('頻率(×pi)');ylabel('H(e(jomega)');set(gca,'XTickMode','manual','XTick',0,wc1,wc2,1);set(gca,'XTickMode',