語音信號濾波去噪——使用脈沖響應不變法設計的巴特沃斯濾波器

1、 王雄杰語音信號濾波去噪-使用脈沖響應不變法設計的巴特沃斯濾波器第20頁 共20頁 語音信號濾波去噪使用脈沖響應不變法設計的巴特沃斯濾波器學生姓名：王雄杰 指導老師：喻小虎摘 要 本課程設計主要內(nèi)容是設計利用脈沖響應不變法設計的巴特沃斯濾波器，對一段含噪語音信號進行濾波去噪處理并根據(jù)濾波前后的波形和頻譜分析濾波性能。本課程設計仿真平臺為MATLAB7.0，開發(fā)工具是M語言編程。錄制一段語音信號，并人為加入一單頻噪聲，然后對信號進行頻譜分析以確定所加噪聲頻率，即 fn=3000Hz，并設計濾波器進行濾波去噪處理，最后比較濾波前后的波形和頻譜，得出結(jié)果為，濾波器后的語音信號與原始信號基本一致，即設

2、計的巴特沃斯濾波器能夠去除信號中所加單頻噪聲，達到了設計目的。關鍵詞 課程設計；濾波去噪；巴特沃斯濾波器；脈沖響應不變法；MATLAB1 引 言本課程設計需要錄制一段頻率為8000Hz，8位的單聲道語音信號，并繪制波形觀察其頻譜，再對其進行加噪處理，同樣要繪制加噪后的頻譜圖，再用MATLAB利用脈沖響應不變法設計一個巴特沃斯濾波器，將該語音信號進行濾波去噪處理。1.1 課程設計目的課程設計是教學的最后一個步驟，課程設計有利于基礎知識的理解，我們掌握了基礎知識和基本技能，但是要真正接觸才能真正理解課程的深入部分；還有利于邏輯思維的鍛煉，在許多常規(guī)學科的日常教學中，我們不難發(fā)現(xiàn)這樣一個現(xiàn)象，不少學

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

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

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

6、最大的特點就是簡潔開放的程序代碼和直觀實用的開發(fā)環(huán)境1。2 設計原理本課程設計需要錄制的一段頻率為8000Hz，8位的單聲道語音信號，并繪制波形觀察其頻譜，再對其進行加噪處理，同樣要繪制加噪后的頻譜圖，再用MATLAB利用脈沖響應不變法設計一個巴特沃斯濾波器，將該語音信號進行濾波去噪處理。2.1 IIR濾波器IIR濾波器具有無限長脈沖響應，因此能夠與模擬濾波器相匹敵；一般來說，所有的模擬濾波器都有無限長脈沖響應。因此，IIR濾波器設計的基本方法是利用復值映射將大家熟知的模擬濾波器變換為數(shù)字濾波器。這一方法的優(yōu)勢在于各種模擬濾波器設計（AFD）表格和映射在文獻中普遍都能獲得。這種IIR濾波器設計

7、的基本方法存在兩種途徑：設計模擬低通濾波器途徑1：實行濾波器變換SZ實行頻帶變換 S S期望的IIR濾波器途徑2：設計模擬低通濾波器實行頻帶變換ZZ實行濾波器變換SZ 期望的IIR濾波器我們將研究途徑二的設計基本方法，在這種IIR濾波器設計方法將按下列步驟進行：（1）設計模擬低通濾波器。（2）研究并實行濾波器變換以得到數(shù)字低通濾波器。（3）研究并實行頻帶變換以便從數(shù)字低通濾波器得到其他數(shù)字濾波器2。2.2 巴特沃斯濾波器巴特沃斯濾波器是用上述IIR濾波器的性質(zhì)表征的：它的幅度響應在通帶和阻帶都是平坦的。一個N階低通濾波器的幅度平方響應給出為： （2-1）式中，N為整數(shù)，是濾波器的階，是截止頻率

8、。稱為濾波器的階數(shù)，N越大，通帶和阻帶的近似性越好，過渡帶也越陡。幅度平方響應的圖如下所示。圖2.1 幅度平方響應圖從這張圖可看出下面幾個性質(zhì)：（1）=1，對全部N。（2），對全部N，這意味著在有3dB衰減。（3）是的單調(diào)下降的函數(shù)。（4）隨N向一個理想低通濾波器趨近。（5）在世最大平坦，因為在這里所有階的導數(shù)存在且等于零。 為了確定系統(tǒng)函數(shù)現(xiàn)將得到（2-2）由（2-2）式分母多項式的根（或）的極點）給出為 （2-3）（2-3）式的解釋是：（1）2N個極點均勻分布在半徑為c的圓上，相隔 pi/N （2）N 為奇數(shù)時， pk= c ej2pik/N（3）N 為偶數(shù)時， pk= c ej(pi/2

9、N+kpi/N)（4）極點對虛軸(j)對稱分布（5）極點永遠不會落在虛軸上，且僅當N為奇數(shù)時才會落在實軸上（6）選取極點全部落在左半平面可得到因果穩(wěn)定的Ha(s)作為一個例子，三階和四階巴特沃斯濾波器的極點分布如圖所示：圖2.2巴特沃斯濾波器的極點圖通過選取在左半面的極點就能給出一個穩(wěn)定和因果的，并且能將寫成：（2-4）2.3 脈沖響應不變法脈沖響應不變法是從濾波器的脈沖響應出發(fā)，使數(shù)字濾波器的單位脈沖響應序h(n)模仿模擬濾波去的沖擊響應ha(t)，使h(n)正好等于ha(t)的采樣值，即h(n)=ha(nT) （2-5）T為采樣周期。如以Ha(s)及H(z)分別表示ha(t)的拉式變換及h

10、(n)的z變換。變換前： （2-6）變換后： （2-7）由于z=是在單位圓上，s=j是在虛軸上，所以有下面從s平面到z平面的變換 z= （2-8） 系統(tǒng)函數(shù)H（z）和Ha（s）是經(jīng)過頻域混疊公式聯(lián)系的： （2-9） 在（2-8）式的映射關系下，復平面的變換如圖2.3 圖2.3脈沖響應不變法中的復平面映射映射特征：（1）記 = Re(s): < 0, 映射到 |z|<1 (單位圓內(nèi))， = 0,映射到 |z|=1 (單位圓上)。 >0,映射到 |z|>1 (單位圓外)。（2）多個 s 映射到一個 z : 多點對一點的映射，每個寬度為2/T 的全部半無限帶都映射到單位圓內(nèi)。

11、（3）因果穩(wěn)定的AF仍然映射為因果穩(wěn)定的DF。（4）如果濾波器不是真正帶限會產(chǎn)生混疊誤差(Aliasing )。3.設計步驟3.1設計流程圖語音信號濾波去噪使用脈沖響應不變法設計的巴特沃斯濾波器，設計流程圖如下3.1所示：開始1、語音信號信號的采集（wavread函數(shù)），畫時域圖快速傅里葉變換，并且畫頻譜圖2、加入3000Hz的噪聲畫出加噪后頻譜圖3、設定濾波器性能指標，通帶截止頻率fp=2600，阻帶截止頻率fs=2950，通帶波紋Rp=3，阻帶波紋As=16脈沖響應不變法法設計巴特沃斯濾波器驗證并進行頻譜分析4、設計好的濾波器進行濾波處理比較濾波前后語音信號的波形及頻譜回放語音信號結(jié)束圖3

12、.1 脈沖響應不變法巴特沃斯濾波器對語音信號去噪流程圖3.2語言信號的采集人工錄制一段為.wav的語音信號，時間在8-9s左右，要求為8000Hz，8位單聲道的音頻格式。然后在Matlab軟件平臺下，利用函數(shù)wavread對語音信號進行采樣，函數(shù)為x,fs,bits=wavread('music.wav'),記住采樣頻率fs和采樣點數(shù)。所下載音樂的屬性設置如圖3.2：圖3.2 語音信號參數(shù)設置3.3語音信號的頻譜分析首先使用MATLAB分別畫出原始語音信號的時域波形，再對原始語音信號進行快速傅里葉變換，畫出原始信號的頻率幅度譜，然后再原始語音信號的基礎上加入單頻噪聲后，噪聲頻率

13、為fn=3000Hz,再所得到的音樂做時間幅度譜，同樣對它做傅里葉變換，得到此信號的頻率頻譜特性?？傻玫綀D形3.3和圖3.4。其中傅里葉變換表示為：X=abs(fft(x); Y=abs(fft(y)圖3.3原始信號與加噪后的程序運行圖由上圖可知，單頻噪聲頻率fn=3000Hz，此數(shù)據(jù)是在對原始語音信號的頻率幅度譜得出，此語音信號的主頻帶在2000Hz內(nèi)，我們?nèi)?000Hz是為了達到更好的干擾效果。因此帶上耳塞后可以聽到原始語音和加噪后的帶有尖銳單頻嘯叫聲的語音?？梢詮穆犛X上感知噪聲的存在，經(jīng)過試驗當fn取2000Hz以內(nèi)的頻率時，聽到的干擾信號無明顯的尖銳單頻嘯叫聲。圖3.4原始信號與加噪后

14、的波形圖與頻譜圖 由圖我們可以看出，在加入噪聲后，原始語音信號和加噪后有不同的時間幅度譜和頻率幅度譜，證明噪聲對原始信號造成了影響。3.4 濾波器設計將數(shù)字濾波器的設計指標設為通帶截止頻率fp=2600Hz,阻帶頻率fc=2950Hz,通帶波紋Rp=3dB，阻帶波紋As=16dB,要求確定H(z)。設計步驟如下：（1） 選取T=1。（2） 利用設參數(shù)，Rp和設計一個模擬濾波器。（3） 利用部分分式將展開。 （4） 現(xiàn)在將模擬極點pk變換為數(shù)字極點epkT，得到數(shù)字濾波器。并作化簡得出作為z-1有理函數(shù)的H(z)。根據(jù)上述步驟，按要求我們用脈沖響應不變法設計好了巴特沃斯濾波器。根據(jù)原理課設計濾波

15、器，可得到如下圖3.5： 3.5 利用脈沖響應不變法設計的數(shù)字巴特沃斯低通濾波器程序運行圖由圖可知，設置fn=3000Hz，即確定了通帶截止頻率為，fp=2600Hz，阻帶截止頻率為，fc=2950Hz，過渡帶寬為450Hz，時間T=1s，為過渡帶寬和時間T的值對于濾波器的設計十分重要。 cs,ds=afd_butt(OmegaP,OmegaS,Rp,As); * Butterworth Filter Order = 15 由此句可得出：濾波器為16階。 >> wp/pians = 0.6500 >> ws/pians = 0.7375 圖3.6 利用脈沖響應不變法設計

16、的數(shù)字巴特沃斯濾波器 由濾波器幅度響應圖可以看出，設計的通帶截止頻率wp/pi=0.6500,阻帶截止頻率ws/pi=0.7375，縱坐標差不多剛好在As=16db處，所以設計的濾波器達到要求。3.5 信號的濾波 前面已經(jīng)用脈沖響應不變法設計好了我們要的巴特沃斯低通濾波器，接著就對語音信號進行濾波處理，看自己設計的巴特沃斯濾波器有沒有對我們的加噪后的語音信號進行處理。所以就用filter函數(shù)進行濾波，即y_fil=filter(b,a,y);。我們將濾波前后的時域波形進行比較，并對其進行快速傅里葉變換，即Y_fil=abs(fft(y_fil);，目的是對比前后的頻域頻譜，具體分析設計的濾波器

17、是否達到設計要求，再將濾波后的信號與原始信號進行對比，可看出濾波效果。如圖3.7是原始語音信號、加干擾后的信號和濾波后信號的時間幅度譜與頻率譜圖：圖3.7原始信號、加噪后信號和濾波后信號的波形圖和頻譜圖觀察分析圖3.7，由原始信號與加噪后的頻率幅度譜分析可看出在fn=3000Hz的頻率上，加了單頻噪聲，由于加入噪聲，對原始信號造成了干擾，在運行sound(y,fs,bits);后，可以明顯聽出有尖銳的單頻嘯叫聲，對加噪后的信號和濾波后的信號進行頻率幅度譜分析，可得出，設計的巴特沃斯低通濾波器濾除了大部分的噪聲，運行sound(y_fil,fs,bits);后，原來的尖銳聲微乎其微，由此證明了設

18、計的合理性。最后對原始和濾波后的頻率幅度譜進行分析，濾波后的信號基本與原始信號無異，因此濾波成功。3.6 結(jié)果分析首先采集語音信號，再按照步驟用脈沖響應不變法設計巴特沃斯濾波器，得到圖3.4。并且由圖3.4可知，設計的通帶截止頻率wp/pi=0.6500，阻帶截止頻率ws/pi=0.7375，縱坐標在As=16db之下，所以設計的濾波器達到要求。我們觀察到圖3.7濾波前后語音信號的波形對比圖，發(fā)現(xiàn)時域波形中的變化不明顯，可能是因為我們采集的語音信號噪聲不是很大，但是還是有濾去噪聲的；但是可以看到在頻域波形中，很明顯地反應出設計的濾波器濾去了我們采集的語音信號中的噪聲。所以，運用脈沖響應不變法設

19、計的巴特沃斯濾波器達到了設計要求。帶上耳塞，運行程序，第一段為原音，第二段為加噪后的音樂，第三段為濾波后的音樂，通過反復聽這三段音樂，能夠聽出第一段音樂和第二段音樂確有不同，能感覺到噪聲，第二段與第三段相對比，幾乎聽不到噪聲，第一段與第三段相比，幾乎無差距。圖3.8，與相同濾波器（脈沖響應不變法的巴特沃斯低通濾波器）不同音樂（人工語音），通過觀察自己的圖形與同學的設計結(jié)果圖，在交流中發(fā)現(xiàn)，我們?nèi)绻O置相同的fn=3000Hz,我們?yōu)V波器階數(shù)相同，都為15階，但是我們的頻率幅度譜的幅度值有差異，但是濾波效果是相同的。圖3.8語音信號的脈沖響應不變法巴特沃斯濾波去噪圖4 出現(xiàn)的問題及解決方法 在這

20、次的課程設計中我們設計的濾波器對語音信號的濾噪處理或多或少會因為我們的疏忽或者實踐能力不夠還是會遇到問題的，當遇到問題時我們不但運用自己學到的知識還得到同學和老師的幫助將其解決，不斷完善和修改。設計中出現(xiàn)的問題及解決方法如下： （1）設計的濾波器達不到要求，因為在采集語音信號時，錄音機的音頻格式應該是PCM 8000Hz，8位，單聲道。 （2）在運用Matlab設計濾波器時，當編輯完前面兩條程序時無法放出聲音，后來發(fā)現(xiàn)我們應當把采集的語音信號wav文件放到Matlab的work文件夾中，如果要放在其他文件夾在調(diào)用此音樂時需標明路徑，如： x,fs,bits=wavread('d:mus

21、ic.wav'); 則標明名字為music的音樂信號時放在d盤。 （3）還要在濾波器性能曲線的wp、ws、Ap處畫一根豎線，這樣更方便看出結(jié)果，其中通帶截止頻率處線的值由wp/pi，阻帶截止頻率處的值為ws/pi。 （4）所有的時間波形橫坐標都要化為時間，濾波前后頻譜的橫坐標應是頻率，這樣在觀察通帶截止頻率和阻帶截止頻率時更加精確，誤差較小。 （5）由于通帶截止頻率fp=fn-400和阻帶截止頻率fc=fn-50的值會影響濾波器的過渡帶寬，而過渡帶寬如果太小會造成濾波器無法達到指標，所以適當?shù)脑谝阎猣n的條件下，基本保證過度帶寬在0.1*pi附近。 （6）時間T也會影響濾波器的設計，開

22、始時我設定T=1/fs，fs=8000Hz，此時采樣頻率過大，造成了音樂信號放大，經(jīng)改正將T設為1s，此時就達到了預知的效果。（7）wp=fp/fs*2*pi;ws=fc/fs*2*pi; OmegaP=wp/T; OmegaS=ws/T; 在用脈沖響應不變法設計低通濾波器時，這些公式的錯誤運用也會對濾波器設計造成嚴重的影響。（8）對于設計的濾波器我們?nèi)绾螜z驗其是否可用，我們要對濾波器的頻率幅度譜進行分析，在ws/pi=0.7375與圖中曲線的交接處以下的值是否大于As=16db，如果大于，則此濾波器設計完成。 （9）在觀察中我們發(fā)現(xiàn)，時間幅度譜對于一些細節(jié)的變化并不明顯，因此我們可以從頻率幅

23、度譜著手分析。5 結(jié)束語 在這次的課程設計中，我們用一周的時間來做設計，其余時間是自己寫與修改報告，在開始前我們需要復習課本上濾波器設計方面的知識，后來在設計中我們又發(fā)現(xiàn)了很多的問題，比如，在課本學習過程中我們以為不太重要的環(huán)節(jié)在實現(xiàn)濾波器效果時，其實是很重要的，任何一個小小的錯誤，比如一個參數(shù)設置的合理性都會造成無法設計出符合標準的濾波器。并且在這次DSP課程設計，不但鞏固了以前學的知識，還學習了新內(nèi)容，在實踐中收獲知識。原來濾去語音中的噪聲我們運用自己學到的知識能做到，并更加鞏固了MATLAB軟件的使用。 MATLAB軟件的熟練使用是很重要的，整個課程設計過程中，我們都要用到它。當然由于我

24、們的實踐能力有限，在過程中，我們需要老師的幫助，在老師的幫助下問題就迎刃而解啦。 這次做的濾波器要濾去語音信號中的噪聲，覺得很有意思，做了1次的課程設計了，我發(fā)現(xiàn)其實在這樣的課程設計中一些細微之處往往是決定成敗的關鍵。其實用MATLAB軟件做實驗是要細心的，因為很多的語法和常量變量的定義我們都要仔細，一個不小心看錯了或者輸入不認真是容易出錯誤。在發(fā)現(xiàn)錯誤時，問了很多同學也問了老師，當發(fā)現(xiàn)錯誤原因時，很多次都是因為自己的不細心造成。但經(jīng)過思考并仔細檢查，最后得以編譯成功。在做實驗是一定要小心謹慎，也許一個小小的語法錯誤和常量變量的定義的錯誤就造成整個程序出現(xiàn)問題，得不到所需的波形，導致實驗結(jié)果不

25、正確。實驗中，同學間的相互交流也是很重要的，比如，一個同學是跟我用相同音樂的，但是她課題為“語音信號濾波去噪使用脈沖響應不變法設計的切比雪夫I型濾波器”我們只是濾波器不同，在結(jié)束后我們彼此交流了兩個不同濾波器的性能，明顯能夠看出切比雪夫I型濾波器濾波效果更好。在這樣的環(huán)境下，我們能夠收獲更多。這次實驗能夠順利完成，有自己的努力，也有同學和老師的幫助，因此，在這里我要衷心地感謝他們的幫助，也感謝學校給我們的機會，讓我們能夠?qū)⒆约簩W到的知識運用到實際中！參考文獻1 張圣勤MATLAB7.0實用教程M北京：機械工程出版社，2006 2 (美)維納K英格爾，（美）約翰G普羅克斯（著）劉樹棠（譯）.數(shù)字

26、信號處理（MATLAB版）M第2版西安：西安交通大學出版社，2008 3 程佩青數(shù)字信號處理教程M北京：清華大學出版社，2002 4 劉敏，魏玲Matlab通信仿真與應用M北京：國防工業(yè)出版社，20015 Sanjit K.Mitra，孫洪，余翔宇數(shù)字信號處理實驗指導書M電子工業(yè)出版社，20056 陳懷琛數(shù)字信號處理教程MATLAB釋義與實現(xiàn)M電子工業(yè)出版社，20087 丁玉美, 高西全數(shù)字信號處理M西安：電子科技大學出版社，2010附錄1：語音信號濾波去噪設計源程序清單% 程序名稱：untitled.m% 程序功能：采用基于脈沖響應不變法設計法，設計巴特沃斯濾波器對含噪語音信號進行濾波去噪處

27、理。% 程序作者：王雄杰/% 最后修改日期：2015-3-27x,fs,bits=wavread('d:music.wav');% 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放N=length(x); % 計算信號x的長度t=0:1/fs:(N-1)/fs;% 計算時間范圍，樣本數(shù)除以采樣頻率deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率范圍X=abs(

28、fft(x); % 對原始信號和加噪信號進行fft變換,取幅度譜X=X(1:N/2);% 截取前半部分subplot(2,2,1) ;plot(t,x) %畫出原始語音信號的時間幅度譜圖xlabel('時間（單位：s）');ylabel('幅度');title('原始語音信號')subplot(2,2,2) ;plot(f,X) %畫出原始信號的頻率幅度譜圖xlabel('頻率（單位：HZ）');ylabel('幅度譜');title('原始語音信號幅度譜圖')fn=3000;% 單頻噪聲頻率x=x

29、' y=x+0.02*sin(fn*2*pi*t); %在原始語音信號基礎上加單頻噪聲后的每個樣本的值Y=abs(fft(y);% 對原始信號和加噪信號進行fft變換,取幅度譜Y=Y(1:N/2); %截取前半部分subplot(2,2,3) ;plot(t,y) %畫出加入單頻干擾后的的時間幅度譜圖axis(0 10 -1 1)% 橫坐標從0到10 ，縱坐標從-1到1xlabel('時間（單位：s）');ylabel('幅度');title('加入單頻干擾后的語音信號')subplot(2,2,4) ;plot(f,Y)%畫出加入單頻干

30、擾后的頻率幅度譜圖xlabel('頻率（單位：HZ）');ylabel('幅度譜');title('加入單頻干擾后的語音信號幅度譜圖')sound(y,fs,bits); % 可以明顯聽出有尖銳的單頻嘯叫聲% 程序名稱：filter.mfp=fn-400;fc=fn-50; %定義通帶和阻帶截止頻率T=1; %定義采樣間隔Rp=3;As=16; % 定義通帶波紋和阻帶衰減wp=fp/fs*2*pi;ws=fc/fs*2*pi; %計算對應的數(shù)字頻率OmegaP=wp/T; OmegaS=ws/T; %截止頻率線性變換 cs,ds=afd_butt

31、(OmegaP,OmegaS,Rp,As); %計算濾波器階數(shù)和截止頻率和濾波器系統(tǒng)函數(shù)分子分母系數(shù)b,a=imp_invr(cs,ds,T);% 脈沖響應不變法得到數(shù)字濾波器系統(tǒng)函數(shù)分子分母系數(shù)db,mag,pha,grd,w=freqz_m(b,a); %驗證濾波器是否達到指定性能delta=1,zeros(1,99);ha=filter(b,a,delta);% 計算脈沖響應Subplot(221);plot(w/pi,db) xlabel('w/pi');ylabel('db');title('濾波器幅度響應圖')Subplot(222)

32、;plot(w/pi,mag)xlabel('w/pi');ylabel('幅度mag');title('濾波器幅度響應圖')Subplot(223);plot(w/pi,pha)xlabel('w/pi');ylabel('相位pha');title('濾波器相位響應圖')Subplot(224);plot(ha)xlabel('n');ylabel('h(n)');title('濾波器脈沖響應圖')% 程序名稱：voisenoisefilter.mx,fs,bits=wavread('d:music.wav');% 輸入?yún)?shù)為文件的全路徑和文件名，輸出的第一個參數(shù)是每個樣本的值，fs是生成該波形文件時的采樣率，bits是波形文件每樣本的編碼位數(shù)。sound(x,fs,bits); % 按指定的采樣率和每樣本編碼位數(shù)回放N=length(x); % 計算信號x的長度t=0:1/fs:(N-1)/fs;% 計算時間范圍，樣本數(shù)除以采樣頻率deltaf=fs/N; % 計算頻譜的譜線間隔f=0:deltaf:fs/2-deltaf; % 計算頻譜頻率范圍X=abs(fft(x); % 對原始信號和加噪信號進行fft變換,