語音信號濾波去噪——使用脈沖響應(yīng)不變法設(shè)計的巴特沃斯濾波器.

1、 語音信號濾波去噪使用脈沖響應(yīng)不變法設(shè)計的巴特沃斯濾波器 第 頁 共 18 頁 語音信號濾波去噪使用脈沖響應(yīng)不變法設(shè)計的巴特沃斯濾波器 摘 要 本課程設(shè)計是利用雙線性變換法設(shè)計的切比雪夫I型濾波器對語音信號濾波去噪。用windows工具中的錄音機采集一段自己說的話（語音信號），然后在MATLAB7.0軟件上，用wavread函數(shù)求出語音信號的三個參數(shù)，對錄制的話（語音信號）進行讀取，并繪制其時域和頻域圖；然后給語音信號加噪聲，再對加了噪聲之后的信號做傅立葉變換，繪制出時域和頻域波形。然后給定相應(yīng)技術(shù)指標，用雙線性變換法設(shè)計的一個滿足指標的切比雪夫I型IIR濾波器，對該語音信號進行濾波去噪處理，

2、比較濾波前后的波形和頻譜并進行分析。通過分析濾波前后語音信號的時域圖和頻域圖發(fā)現(xiàn)，噪聲基本被濾掉，然后播放濾波后的語音文件，發(fā)現(xiàn)噪聲基本被消除。 關(guān)鍵詞 濾波去噪；切比雪夫I型IIR濾波器；雙線性變換法；MATLAB1 引 言本課程設(shè)計主要利用麥克風采集一段8000Hz的單聲道語音信號，并繪制波形觀察其頻譜，再用MATLAB利用雙線性變換法設(shè)計的切比雪夫I型濾波器，將該語音信號進行濾波去噪處理。1.1 課程設(shè)計目的課程設(shè)計是教學的最后一個步驟，課程設(shè)計有利于基礎(chǔ)知識的理解，我們掌握了基礎(chǔ)知識和基本技能，但是要真正接觸才能真正理解課程的深入部分；還有利于邏輯思維的鍛煉，在許多常規(guī)學科的日常教學中

3、，我們不難發(fā)現(xiàn)這樣一個現(xiàn)象，不少學生的思維常常處于混亂的狀態(tài)，寫起作文來前言不搭后語，解起數(shù)學題來步驟混亂，這些都是缺乏思維訓練的結(jié)果，所以我們可以通過實踐來分析問題、解決問題、預(yù)測目標等目的；同時也有利于與其他學科的整合，例如我們這次的課程設(shè)計就要運用MATLAB軟件的幫助才能實現(xiàn)；最重要的有利于治學態(tài)度的培養(yǎng)，在課程設(shè)計中，我們可能經(jīng)常犯很多小錯誤，可能要通過好幾次的反復(fù)修改、調(diào)試才能成功，但這種現(xiàn)象會隨著學校的深入而慢慢改觀。這當中就有一個嚴謹治學、一絲不茍的科學精神的培養(yǎng)，又有一個不怕失敗、百折不饒品格的鍛煉。數(shù)字信號處理課程設(shè)計是在學生完成數(shù)字信號處理和MATLAB的結(jié)合后的基本實驗

4、以后開設(shè)的。本課程設(shè)計的目的是為了讓學生綜合數(shù)字信號處理和MATLAB并實現(xiàn)一個較為完整的小型濾波系統(tǒng)。這一點與驗證性的基本實驗有本質(zhì)性的區(qū)別。開設(shè)課程設(shè)計環(huán)節(jié)的主要目的是通過系統(tǒng)設(shè)計、軟件仿真、程序安排與調(diào)試、寫實習報告等步驟，使學生初步掌握工程設(shè)計的具體步驟和方法，提高分析問題和解決問題的能力，提高實際應(yīng)用水平。1.2 課程設(shè)計的要求（1）濾波器指標必須符合工程設(shè)計。（2）設(shè)計完后應(yīng)檢查其頻率響應(yīng)曲線是否滿足指標。（3）處理結(jié)果和分析結(jié)論應(yīng)該一致，而且應(yīng)符合理論。（4）獨立完成課程設(shè)計并按要求編寫課程設(shè)計報告。1.3 設(shè)計平臺本課程設(shè)計的平臺為MATLAB7.0。MATLAB是美國Math

5、Works公司生產(chǎn)的一個為科學和工程計算專門設(shè)計的交互式大型軟件，是一個可以完成各種精確計算和數(shù)據(jù)處理的、可視化的、強大的計算工具。它集圖示和精確計算于一身，在應(yīng)用數(shù)學、物理、化工、機電工程、醫(yī)學、金融和其他需要進行復(fù)雜數(shù)值計算的領(lǐng)域得到了廣泛應(yīng)用。它不僅是一個在各類工程設(shè)計中便于使用的計算工具，而且也是一個在數(shù)學、數(shù)值分析和工程計算等課程教學中的優(yōu)秀的數(shù)學工具，在世界各地的高和大型計算機上運行，適用于Windows、UNIX等多種系統(tǒng)平臺。MATLAB作為一種科學計算的高級語言之所以受歡迎，就是因為它有豐富的函數(shù)資源和工具箱資源，編程人員可以根據(jù)自己的需要選擇函數(shù)，而無需再去編寫大量繁瑣的程

6、序代碼，從而減輕了編程人員的工作負擔，被稱為第四代編程語言的MATLAB最大的特點就是簡潔開放的程序代碼和直觀實用的開發(fā)環(huán)境1。2 設(shè)計原理用麥克風采集一段語音信號，繪制波形并觀察其頻譜，給定相應(yīng)技術(shù)指標，用雙線性變換法設(shè)計的一個滿足指標的切比雪夫I型濾波器，對該語音信號進行濾波去噪處理，比較濾波前后的波形和頻譜并進行分析。2.1 IIR濾波器從離散時間來看，若系統(tǒng)的單位抽樣(沖激)響應(yīng)延伸到無窮長，稱之為“無限長單位沖激響應(yīng)系統(tǒng)”，簡稱為IIR系統(tǒng)。（1）無限長單位沖激響應(yīng)（IIR）濾波器有以下幾個特點:（2）系統(tǒng)的單位沖激響應(yīng)h(n)是無限長；（3）系統(tǒng)函數(shù)H(z)在有限z平面（0<

7、<）；結(jié)構(gòu)上存在著輸出到輸入的反饋，也就是結(jié)構(gòu)上是遞歸型的。IIR濾波器采用遞歸型結(jié)構(gòu)，即結(jié)構(gòu)上帶有反饋環(huán)路。同一種系統(tǒng)函數(shù)H（z）可以有多種不同的結(jié)構(gòu)，基本網(wǎng)絡(luò)結(jié)構(gòu)有直接型、直接型、級聯(lián)型、并聯(lián)型四種，都具有反饋回路。同時，IIR數(shù)字濾波器在設(shè)計上可以借助成熟的模擬濾波器的成果，巴特沃斯(Butterworth)濾波器、切比雪夫(Chebyshev)濾波器、橢圓(Cauer)濾波器、貝塞爾(Bessel)濾波器等，這些典型的濾波器各有特點。有現(xiàn)成的設(shè)計數(shù)據(jù)或圖表可查，在設(shè)計一個IIR數(shù)字濾波器時，我們根據(jù)指標先寫出模擬濾波器的公式，然后通過一定的變換，將模擬濾波器的公式轉(zhuǎn)換成數(shù)字濾波器

8、的公式。2.2 切比雪夫I型濾器切比雪夫濾波器（又譯車比雪夫濾波器）是在通帶或阻帶上頻率響應(yīng)幅度等波紋波動的濾波器。在通帶波動的為“I型切比雪夫濾波器”，在阻帶波動的為“II型切比雪夫濾波器”。切比雪夫濾波器在過渡帶比巴特沃斯濾波器的衰減快，但頻率響應(yīng)的幅頻特性不如后者平坦。切比雪夫濾波器和理想濾波器的頻率響應(yīng)曲線之間的誤差最小，但是在通頻帶內(nèi)存在幅度波動。 這種濾波器來自切比雪夫多項式，因此得名，用以紀念俄羅斯數(shù)學家巴夫尼提·列波維其·切比雪夫切比雪夫濾波器的 在通帶范圍內(nèi)是等幅起伏的，所以在同樣的通常內(nèi)衰減要求下，其階數(shù)較巴特沃茲濾波器要小。 切比雪夫濾波器的振幅平方函

9、數(shù)為 (2-1) 式中   c有效通帶截止頻率    與通帶波紋有關(guān)的參量， 大,波紋大 0< <1    VN（x）N階切比雪夫多項式 （2-2）|x|1時,|VN(x)|1 |x|>1時, |x|, VN(x) 切比雪夫濾波器的振幅平方特性如圖所示,通帶內(nèi)， 的變化范圍為 1(max) (min)時，|x|>1,隨 ， 0 (迅速趨于零)當 =0時， （2-3） N為偶數(shù)，cos2( )=1，得到min， （2-4）N為奇數(shù)，cos2( ，得到max， （2-5）2.3 雙線性變換法 雙

10、線性變換法是使數(shù)字信號濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)相似的一種變換方法。為了客服多值映射這一缺點，我們首先把整個s平面壓縮變換到某一中介的s1平面的一條橫帶里（寬度為,即從到),其次再通過上面討論過的標準變換關(guān)系將此橫帶變換到這個z平面上去，這樣就使s平面與z平面式一一對應(yīng)的關(guān)系，消除了多值變換性，也就消除了頻譜混疊現(xiàn)象。將s平面整個軸壓縮變換到s1平面軸上的到一段，可以采用以下變換關(guān)系： （2-6） 這樣，變換到，變到可將（2-6）式寫成 （2-7）解析延拓到整個s平面和s1平面，令，則得 （2-8）再將s1平面通過以下標準變換關(guān)系映射到z平面： （2-9）從而得到s平面和z平面的單

11、值映射的關(guān)系為 （2-10） （2-11）一般來說，為了使模擬濾波器的某一頻率與數(shù)字濾波器的任一頻率有對應(yīng)的關(guān)系，可以引入待定常數(shù)c，使（2-6）式和（2-7）式變換成 （2-12） （2-13）仍將代入（2-13）式，可得 （2-14） （2-15） （2-14）式和（2-15）式是s平面與z平面之間的單值映射關(guān)系，這種變換就稱為雙線性變換。3.設(shè)計步驟3.1設(shè)計流程圖語音信號濾波去噪使用脈沖不變響應(yīng)法設(shè)計的巴特沃斯濾波器的設(shè)計流程如圖3.1所示：開始語音信號的采集（wavread函數(shù)），畫時域圖快速傅里葉變換，并且畫頻譜圖，噪聲頻率fn=2500設(shè)定濾波器性能指標，通帶截止頻率fp=230

12、0，阻帶截止頻率fc=2450，通帶波紋Rp=3.3，阻帶波紋As=15雙線性變換法設(shè)計切比雪夫I型濾波器驗證并進行頻譜分析設(shè)計好的濾波器進行濾波處理比較濾波前后語音信號的波形及頻譜回放語音信號結(jié)束圖3.1 脈沖響應(yīng)不變法巴特沃斯濾波器對語音信號去噪流程圖3.2語音信號的采集 利用Windows下的錄音機，錄制語音信號，時間在1s左右，要求為8000HZ，8位單聲道的音頻格式。然后在Matlab軟件平臺下，利用函數(shù)wavread對語音信號進行采樣，函數(shù)為x,fs,bits= wavread('yuyin.wav'),記住采樣頻率和采樣點數(shù)，如圖3.2：圖3.2 語音信號設(shè)置3.

13、3語音信號的頻譜分析首先畫出語音信號的時域波形，再對語音信號進行快速傅里葉變換，得到信號的頻譜特性。程序如下：>> x,fs,bits= wavread('yuyin.wav');>> sound(x,fs,bits);>> N=length(x);>> fn=2500;>> t=0:1/fs:(N-1)/fs;>> x=x' y=x+0.03*sin(fn*2*pi*t);>> sound(y,fs,bits);>> X=abs(fft(x);Y=abs(fft(y);&g

14、t;> X=X(1:N/2); Y=Y(1:N/2);>> deltaf=fs/N;>> f=0:deltaf:fs/2-deltaf; 得到原始語音信號時域圖形，原始語音信號幅度譜，加入噪聲之后的語音信號時域圖形，加入噪聲之后的語音信號幅度譜如圖3.3所示：圖3.3 原始語音信號和加入噪聲之后的信號的時域波形和幅度譜3.4濾波器設(shè)計將數(shù)字濾波器的設(shè)計指標設(shè)為通帶截止頻率fp=2300HZ,阻帶頻率fc=2450HZ,通帶波紋Rp=3.5dB，阻帶波紋As=15dB,要求確定H(z)。fp=fn-200;fc=fn-50; % 低通濾波器設(shè)計指標wp=fp/fs*

15、2*pi;ws=fc/fs*2*pi; % 將Hz為單位的模擬頻率換算為rad為單位的數(shù)字頻率T=1;OmegaP=wp/T; OmegaS=ws/T; %定義采樣間隔，截止頻率線性變換Rp=25;As=35; % 定義通帶波紋和阻帶衰減T=1;Fs=1/T; OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);% Analog chebyshev-1Prototype Filter Calculation;cs,ds=afd_chb1(OmegaP, OmegaS,Rp,As);% Bilinear transformation:b,a=bilinea

16、r (cs,ds,Fs);db,mag,pha,grd,w = freqz_m(b,a); %驗證濾波器是否達到指定性能n=0:100;h=impz(b,a,n);所以得到圖形如圖3.4：圖3.4 利用脈沖響應(yīng)不變法設(shè)計的數(shù)字巴特沃斯濾波器（w(單位：)）3.5 信號濾波信號前面已經(jīng)用雙線性變換法設(shè)計好了我們要的切比雪夫I型濾波器，接著就對語音信號進行濾波處理，看自己設(shè)計的切比雪夫I型濾波器有沒有對我們的語音信號進行處理。所以就用filter函數(shù)進行濾波，即y_fil=filter(h,1,y)。我們將濾波前后的時域波形進行比較，并對其進行快速傅里葉變換，即Y_fil=abs(fft(y_fi

17、l)，目的是對比前后的頻域頻譜，具體分析設(shè)計的濾波器是否達到設(shè)計要求。如圖3.5濾波前后的語音信號時域?qū)Ρ龋?gt;> y_fil=filter(h,1,y);>> Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2);圖3.5 濾波前后的時域波形和頻譜圖對比3.6 結(jié)果分析 我們先采集語音信號，再按照步驟用雙線性變換法設(shè)計切比雪夫I型濾波器，得到圖3.4。并且由圖3.5可知，縱坐標差不多剛好在As處，所以設(shè)計的濾波器達到要求。我們觀察到圖3.5濾波前后語音信號的波形對比圖，發(fā)現(xiàn)時域波形中的變化不明顯，可能是因為我們采集的語音信號噪聲不是很大，但是還

18、是有濾去噪聲的；但是可以看到在頻域波形中，很明顯地反應(yīng)出設(shè)計的濾波器濾去了我們采集的語音信號中的噪聲。所以，運用雙線性變換法設(shè)計的切比雪夫I型濾波器達到了設(shè)計要求。4.出現(xiàn)的問題及解決方法在這次的課程設(shè)計中，由于理論知識的不踏實以及其他各種原因，我們遇到了不少問題。（1）在進行語音信號提取時，進過多次錄取才得到理想的語音信號，在得到理想的波形時，通過多次嘗試，和查找書籍及同學討論，最后猜得到理想的語音信號的時域圖和頻譜圖（2）在運用Matlab設(shè)計濾波器時，當編輯完前面兩條程序時無法放出聲音，后來發(fā)現(xiàn)我們應(yīng)當把采集的語音信號wav文件放到Matlab的work文件夾中。（3）所有的時間波形橫坐

19、標都要化為時間，濾波前后頻譜的橫坐標應(yīng)是頻率，這樣在觀察通帶截止頻率和阻帶截止頻率時更加精確，誤差較小。5.結(jié)束語 兩周的課程設(shè)計即將結(jié)束，這次做的濾波器要濾去語音信號中的噪聲，覺得很有成就，做了2次的課程設(shè)計了，當自己親自來做的時候就會發(fā)現(xiàn)自己還存在許多缺陷。其實用MATLAB軟件做實驗是要細心的，因為很多的語法和常量變量的定義我們都要仔細，一個不小心看錯了或者輸入不認真是容易出錯誤。也許一個小小的語法錯誤和常量變量的定義的錯誤就造成整個程序出現(xiàn)問題，得不到所需的波形，導致實驗結(jié)果不正確。 經(jīng)過這次課程設(shè)計，讓我有機會將自己學到的理論知識運用到實際中，提高了自己的動手能力和思維能力。在課程設(shè)

20、計中發(fā)現(xiàn)自己的不足，所以在今后的學習和生活中我們要更加努力，學習好我們的專業(yè)知識并要能運用到實際。 這次的設(shè)計在老師的指導和同學的幫助下通過自己的努力和思考成功的完成了。能將自己平時學到的東西能運用到實際中，讓理論和實際得以結(jié)合還是很不錯的。也讓我在課程設(shè)計中找到了動手的樂趣和思考的快樂，很有成就感。希望這次的經(jīng)歷能讓我們在以后的學習生活中不斷成長。最后，在此衷心地感謝老師和同學對我的幫助，也感謝學校給我們的機會，讓我們能夠?qū)⒆约簩W到的知識運用到實際中！參考文獻1 張圣勤.MATLAB7.0實用教程M .機械工程出版社.2006.3 2 維納·K·英格爾，約翰·G

21、·普羅克斯（著），劉樹棠（譯）.數(shù)字信號處理（MATLAB版）M .西安交通大學出版社.2008.13 吳鎮(zhèn)揚.數(shù)字信號處理M .高等教育出版社.2004.94 陳懷琛.數(shù)字信號處理教程：MATLAB釋義與實現(xiàn)M .電子工業(yè)出版社.2004.125 羅軍輝.MALAB7.0在數(shù)字信號處理中的應(yīng)用M .機械工業(yè)出版社.2005.5附錄1：源程序清單% 程序名稱：語音信號濾波去噪% 程序功能：用雙線性變換法設(shè)計的切比雪夫I型濾波器，并對加了噪聲后的語音信號進行濾波去噪。% 程序作者：姜成林% 程序修改日期：2011-3-4>> x,fs,bits= wavread('

22、e:yuyin.wav'); % 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放N=length(x); % 計算信號x的長度fn=2500; % 單頻噪聲頻率t=0:1/fs:(N-1)/fs; % 計算時間范圍，樣本數(shù)除以采樣頻率x=x' y=x+0.03*sin(fn*2*pi*t);%將原語音信號中加入噪聲sound(y,fs,bits); %可以明顯聽出有尖銳的單頻嘯叫聲X=abs(fft(x);Y=abs(ff

23、t(y); % 對原始信號和加噪信號進行fft變換,取幅度譜X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率范圍subplot(2,2,1);plot(t,x);xlabel('時間（單位：s）');ylabel('幅度');title('原始語音信號')subplot(2,2,2);plot(f,X);xlabel('頻率（單位：Hz）');ylabel('幅度譜');title(

24、9;原始語音信號幅度譜') ;subplot(2,2,3);plot(t,y);xlabel('時間（單位：s）');ylabel('幅度');title('加入單頻干擾后的語音信號')subplot(2,2,4);plot(f,Y);xlabel('頻率（單位：Hz）');ylabel('幅度譜');title('加入干擾后的語音信號幅度譜');fp=fn-200;fc=fn-50; % 低通濾波器設(shè)計指標wp=fp/fs*2*pi;ws=fc/fs*2*pi; % 將Hz為單位的模擬頻率

25、換算為rad為單位的數(shù)字頻率T=1;OmegaP=wp/T; OmegaS=ws/T; %定義采樣間隔，截止頻率線性變換Rp=1;As=15; % 定義通帶波紋和阻帶衰減T=1;Fs=1/T; OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);% Analog chebyshev-1Prototype Filter Calculation;cs,ds=afd_butt(OmegaP, OmegaS,Rp,As);% Bilinear transformation:b,a=imp_invr(cs,ds,Fs);db,mag,pha,grd,w = fre

26、qz_m(b,a); %驗證濾波器是否達到指定性能n=0:100;h=impz(b,a,n);subplot(2,2,1);plot(w/pi,mag);xlabel('w');ylabel('|H|');grid on; title('濾波器幅度圖')subplot(2,2,2);plot(w/pi,pha);xlabel('w');ylabel('相位');grid on; title('濾波器相位響應(yīng)圖')subplot(2,2,3);plot(w/pi,db);xlabel('w&#

27、39;);ylabel('db');grid on; title('濾波器幅度響應(yīng)圖')hold on;line(0,1,-35,-35,'linestyle',':','color','r');line(0.18,0.18,-120,20,'linestyle',':','color','r');line(0.35,0.35,-120,20,'linestyle',':','color'

28、,'r'); hold off;subplot(2,2,4);stem(n,h);xlabel('n');ylabel('h(n)');grid on; title('濾波器脈沖響應(yīng)圖')y_fil=filter(h,1,y); % 用設(shè)計好的濾波器對y進行濾波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 計算頻譜取前一半subplot(3,2,1);plot(t,x);xlabel('時間（單位：s）');ylabel('幅度');title('原始語音信號')subplot(3,2,2);plot(f,X);xlabel('頻率（單位：Hz）');ylabel('幅度譜');title('原始語音信號幅度譜') ;subplot(3,2,3);plot(t,y);xlabel('時間（單位：s）');ylabel('幅度');title('加入單頻干擾后的語音信號')subplot(3