：基于MATLAB的IIR濾波器的語音信號去噪

1、摘 要濾波器設(shè)計在數(shù)字信號處理中占有極其重要的地位，本次課程設(shè)計主要是錄制一段語音信號對其進行加噪處理，然后利用IIR低通濾波器對加有隨機噪聲的語音信號進行濾波處理及時頻譜分析，畫出濾波之后的頻譜圖與時域波形，并對信號濾波處理前后進行分析比較，分析信號的變化。通過對對所設(shè)計濾波器的仿真和頻率特性分析，由仿真結(jié)果可以看出，所設(shè)計的濾波器能夠?qū)崿F(xiàn)對語音信號的語音有效去噪，并對濾波前后的語音信號進行對比。 關(guān)鍵詞： 去噪；濾波器；MATLAB目 錄一 語音信號去噪的設(shè)計任務(wù)3二 語音信號去噪的基本原理42.1 數(shù)字濾波器的基本設(shè)計方法42.2 雙線性變換法52.3數(shù)字濾波器設(shè)計基本思想62.4 數(shù)字

2、濾波器的設(shè)計步驟72.5采樣定理8三 基于MATLAB的仿真結(jié)果及結(jié)果分析103.1 IIR高通濾波器的仿真103.2 原始語音信號的錄制103.3 語音信號的時頻域分析113.4 加隨機噪聲后的時頻域分析123.5 濾波前后的時頻域比較15總結(jié)17參考文獻18致謝19附錄20 一 基本原理1.1 數(shù)字濾波器的基本設(shè)計方法IIR數(shù)字濾波器的設(shè)計一般有兩種方法：一個是借助模擬濾波器的設(shè)計方法進行。其設(shè)計步驟是，先設(shè)計模擬濾波器，再按照某種方法轉(zhuǎn)換成數(shù)字濾波器。這種方法比較容易一些，因為模擬濾波器的設(shè)計方法已經(jīng)非常成熟，不僅有完整的設(shè)計公式，還有完善的圖表供查閱；另外一種直接在頻率或者時域內(nèi)進行，

3、由于需要解聯(lián)立方程，設(shè)計時需要計算機做輔助設(shè)計。其設(shè)計步驟是：先設(shè)計過渡模擬濾波器得到系統(tǒng)函數(shù)，然后將按某種方法轉(zhuǎn)換成數(shù)字濾波器的系統(tǒng)函數(shù)1。為了保證轉(zhuǎn)換后的穩(wěn)定且滿足技術(shù)指標要求，對轉(zhuǎn)換關(guān)系提出兩點要求:（1）因果穩(wěn)定的模擬濾波器轉(zhuǎn)換成數(shù)字濾波器，仍是因果穩(wěn)定的。（2）數(shù)字濾波器的頻率相應(yīng)模仿模擬濾波器的頻響特性，s平面的虛軸映射為z平面的單位圓，相應(yīng)的頻率之間呈線性關(guān)系。利用模擬濾波器成熟的理論設(shè)計IIR數(shù)字濾波器的過程是：(1)確定數(shù)字低通濾波器的技術(shù)指標：通帶邊界頻率、通帶最大衰減、阻帶截止頻率、阻帶最小衰減。(2)將數(shù)字低通濾波器的技術(shù)指標轉(zhuǎn)換成相應(yīng)的模擬低通濾波器的技術(shù)指標。(3)

4、按照模擬低通濾波器的技術(shù)指標設(shè)計過渡模擬低通濾波器。(4)用所選的轉(zhuǎn)換方法，將模擬濾波器轉(zhuǎn)換成數(shù)字低通濾波器系統(tǒng)函數(shù)。IIR數(shù)字濾波器的設(shè)計流程圖2-1如下： 模擬濾波器 技術(shù)指標數(shù)字濾波器技術(shù)指標 變換 模擬濾波器 設(shè)計方法 =（ ） 模擬濾波器 數(shù)字濾波器 變換 S=f（Z）圖2-1 IIR數(shù)字濾波器的設(shè)計步驟流程圖1成熟的模擬濾波器設(shè)計方法主要有脈沖響應(yīng)不變法和雙線性變換法。2.2 雙線性變換法脈沖響應(yīng)不變法的主要缺點是產(chǎn)生頻率響應(yīng)的混疊失真。這是因為從S平面到平面是多值的映射關(guān)系所造成的。為了克服這一缺點，可以采用非線性頻率壓縮方法，將整個頻率軸上的頻率范圍壓縮到-/T/T之間，再用z

5、=esT轉(zhuǎn)換到Z平面上。也就是說，第一步先將整個S平面壓縮映射到S1平面的-/T/T一條橫帶里；第二步再通過標準變換關(guān)系z=es1T將此橫帶變換到整個Z平面上去。這樣就使S平面與Z平面建立了一一對應(yīng)的單值關(guān)系，消除了多值變換性，也就消除了頻譜混疊現(xiàn)象，映射關(guān)系如圖2-2所示2。圖2-2雙線性變換的映射關(guān)系3為了將S平面的整個虛軸j壓縮到S1平面j1軸上的-/T到/T段上，可以通過以下的正切變換實現(xiàn) （2-1）式中,T仍是采樣間隔。當1由-/T經(jīng)過0變化到/T時，由-經(jīng)過0變化到+，也即映射了整個j軸。將式（2-1）寫成 （2-2） 將此關(guān)系解析延拓到整個S平面和S1平面，令j=s，j1=s1，

6、則得 （2-3）再將S1平面通過以下標準變換關(guān)系映射到Z平面 （2-4）從而得到S平面和Z平面的單值映射關(guān)系為： （2-5） （2-6）式（2-4）與式（2-5）是S平面與Z平面之間的單值映射關(guān)系4，這種變換都是兩個線性函數(shù)之比，因此稱為雙線性變換式（2-5）與式（2-6）的雙線性變換符合映射變換應(yīng)滿足的兩點要求。首先,把z=ej，可得 （2-7）即S平面的虛軸映射到Z平面的單位圓。其次，將s=+j代入式（2-7），得 （2-8）因此 （2-9）由此看出，當0時，|z|0時，|z|1。也就是說，S平面的左半平面映射到Z平面的單位圓內(nèi)，S平面的右半平面映射到Z平面的單位圓外，S平面的虛軸映射到Z

7、平面的單位圓上。因此，穩(wěn)定的模擬濾波器經(jīng)雙線性變換后所得的數(shù)字濾波器也一定是穩(wěn)定的。 1.3數(shù)字濾波器設(shè)計基本思想 一個數(shù)字濾波器可用它的系統(tǒng)函數(shù)H(z)來描述，或者用一個N階差分方程來描述，即6因此，設(shè)計一個數(shù)字濾波器，實質(zhì)上是尋找一組系數(shù)ak,br ，使其性能滿足預(yù)定的技術(shù)要求，它是一個數(shù)學逼近問題,顯然它與模擬濾波器的設(shè)計方法是完全一致的，只不過模擬濾波器的設(shè)計是在Z平面上用數(shù)學逼近方法尋找近似于所需特性的H(s)，而數(shù)字濾波器的設(shè)計則在Z平面上尋找合適的H(z)。確定了ak,br，剩下的問題是設(shè)計一個具體的網(wǎng)絡(luò)結(jié)構(gòu)去實現(xiàn)它?？梢姅?shù)字濾波器設(shè)計的基本步驟如下：（1）確定指標在設(shè)計一個數(shù)字

8、濾波器之前，必須首先根據(jù)工程實際的需要確定數(shù)字濾波器的技術(shù)指標。在很多實際應(yīng)用中，數(shù)字濾波器常常用來實現(xiàn)選頻操作。因此，指標一般在頻域中給出，諸如通帶截止頻率wp、阻帶截止頻率ws、阻帶內(nèi)允許的最大衰減ap、阻帶內(nèi)允許的最小衰減as 等。此外還必須確定采樣周期T或采樣頻率Fs。 （2）逼近確定了技術(shù)指標后，就可以建立一個目標數(shù)字濾波器模型。通常采用理想的數(shù)字濾波器模型。之后，利用數(shù)字濾波器的設(shè)計方法，設(shè)計出一個實際濾波器模型來逼近給定的目標。 （3）性能分析和計算機仿真上兩步的結(jié)果是得到以系統(tǒng)函數(shù)H(z)或單位沖激響應(yīng)h(n)描述的數(shù)字濾波器。根據(jù)這個描述就可以分析其頻率特性和相位特性，以驗證

9、設(shè)計結(jié)果是否滿足指標要求；或者利用計算機仿真實現(xiàn)設(shè)計的濾波器，再分析濾波結(jié)果來判斷。數(shù)字濾波器根據(jù)其單位沖激響應(yīng)函數(shù)的時域特性，可分為兩種，即無限長沖激響應(yīng)（IIR）濾波器和有限長沖激響應(yīng)（FIR）濾波器。IIR濾波器的特征是具有無限持續(xù)時間沖激響應(yīng)。這種濾波器一般需要用遞歸模型來實現(xiàn)，因而有時也稱之為遞歸濾波器。FIR濾波器的沖激響應(yīng)只能延續(xù)一定時間，在工程實際中可以采用遞歸的方式實現(xiàn)，也可以采用非遞歸的方式實現(xiàn)。數(shù)字濾波器的設(shè)計方法有多種，如脈沖響應(yīng)不變法、雙線性變換法、窗函數(shù)設(shè)計法、插值逼近法和Chebyshev逼近法等等 5。1.4 數(shù)字濾波器的設(shè)計步驟（1）確定所需類型數(shù)字濾波器的技

10、術(shù)指標：通帶邊界頻率Fp、通帶最大衰減As，阻帶截止頻率Fc、阻帶最小衰減Ap。（2）將所需類型數(shù)字濾波器的邊界頻率轉(zhuǎn)換成相應(yīng)的模擬濾波器的邊界頻率，轉(zhuǎn)換公式為=2/T tan(0.5)(3)將相應(yīng)類型的模擬濾波器技術(shù)指標轉(zhuǎn)換成模擬低通濾波器技術(shù)指標。（4）設(shè)計模擬濾波器。（5）通過頻率變換將模擬低通轉(zhuǎn)換成相應(yīng)類型的過渡模擬濾波器。（6）采用雙線性變換法將相應(yīng)類型的過渡模擬濾波器轉(zhuǎn)換成所需類型的數(shù)字濾波器。 MATLAB信號處理工具箱函數(shù)buttp buttor butter是巴特沃斯濾波器設(shè)計函數(shù)，其有5種調(diào)用格式，本課程設(shè)計中用到的是N,wc=butter(N,wc,Rp,As,s),該格

11、式用于計算巴特沃斯模擬濾波器的階數(shù)N和3dB截止頻率wc。MATLAB信號處理工具箱函數(shù)cheblap,cheblord和cheeby1是切比雪夫I型濾波器設(shè)計函數(shù)。我們用到的是cheeby1函數(shù)，其調(diào)用格式如下：B,A=cheby1(N,Rp,wpo,ftypr)B,A=cheby1(N,Rp,wpo,ftypr,s)函數(shù)butter,cheby1和ellip設(shè)計IIR濾波器時都是默認的雙線性變換法，所以在設(shè)計濾波器時只需要代入相應(yīng)的實現(xiàn)函數(shù)即可.1.5采樣定理(1)時域采樣定理 頻帶為F的連續(xù)信號f(t)可用一系列離散的采樣值f(t1),f(t1t)，f(t12t)，.來表示,只要這些采樣

12、點的時間間隔t1/2F，便可根據(jù)各采樣值完全恢復(fù)原來的信號f(t)8。 這是時域采樣定理的一種表述方式。 時域采樣定理的另一種表述方式是：當時間信號函數(shù)f(t)的最高頻率分量為fM時,f(t)的值可由一系列采樣間隔小于或等于1/2fM的采樣值來確定,即采樣點的重復(fù)頻率f2fM。圖為模擬信號和采樣樣本的示意圖。 時域采樣定理是采樣誤差理論、隨機變量采樣理論和多變量采樣理論的基礎(chǔ)。(2)頻域采樣定理 對于時間上受限制的連續(xù)信號f(t)（即當tT時,f(t)=0,這里T=T2-T1是信號的持續(xù)時間），若其頻譜為F（）,則可在頻域上用一系列離散的采樣值 來表示,只要這些采樣點的頻率間隔 / tm 。(

13、3) 采樣頻率采樣頻率，也稱為采樣速度或者采樣率，單位為赫茲（Hz）。定義了每秒從連續(xù)信號中提取并組成離散信號的采樣個數(shù)，常用的表示符號是 fs采樣時間，采樣頻率的倒數(shù)是采樣周期，是采樣之間的時間間隔。通俗的講采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質(zhì)、音調(diào)，衡量聲卡、聲音文件的質(zhì)量標準。通俗的講采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質(zhì)、音調(diào)，衡量聲卡、聲音文件的質(zhì)量標準。采樣頻率越高，即采樣的間隔時間越短，則在單位時間內(nèi)計算機得到的聲音樣本數(shù)據(jù)就越多，對聲音波形的表示也越精確。根據(jù)采樣定理，只有采樣頻率高于聲音信號最高頻率的兩倍時，才能把數(shù)字信號

14、表示的聲音還原成為原來的聲音。這就是說采樣頻率是衡量聲卡采集、記錄和還原聲音文件的質(zhì)量標準。采樣位數(shù)和采樣率對于音頻接口來說是最為重要的兩個指標，也是選擇音頻接口的兩個重要標準。無論采樣頻率如何，理論上來說采樣的位數(shù)決定了音頻數(shù)據(jù)最大的力度范圍。每增加一個采樣位數(shù)相當于力度范圍增加了6dB。采樣位數(shù)越多則捕捉到的信號越精確。 二 實現(xiàn)框圖本次課程設(shè)計主要任務(wù)是對語音信號的簡單處理。運用數(shù)字信號學基本原理實現(xiàn)語音信號的處理.其大概流程框圖可如下表示：語音信號采集效果顯示、對比語音信號錄入語音信號變換信號加噪語音信號濾波三 基于MATLAB的仿真結(jié)果及結(jié)果分析3.1 IIR高通濾波器的仿真程序見附

15、錄（1）圖3-1 高通濾波器3.2 原始語音信號的錄制1.基本步驟（1）設(shè)置采樣頻率Fs=11025；（2）利用函數(shù)“wavrecord（15*Fs,Fs,int16）”命令進行錄音；（3）利用函數(shù)wavplay()重新播放，看看是否滿足所要的效果；（4）利用“wavwrite”進行保存，應(yīng)記住保存路徑，以便后面進行調(diào)用時能快速找到。需要指出的是利用MATLAB錄的音已經(jīng)是通過采樣得到的數(shù)字信號2.程序 Fs = 11025;%采樣頻率為11025 y = wavrecord(15*Fs, Fs, int16);%進行15秒-16bit的錄音 wavplay(y,Fs);%重新播放錄音 wav

16、write(y,Fs,xxx.wav);保存該錄音，且名字為xxx,格式為wav wavwrite(y,xxx.wav);注：Wavread 函數(shù)幾種調(diào)用格式。（1） y=wavread（file）功能說明：讀取file所規(guī)定的wav文件，返回采樣值放在向量y中。（2）y,fs,nbits=wavread(file) 功能說明：采樣值放在向量y中，fs表示采樣頻率（hz），nbits表示采樣位數(shù)。（3）y=wavread（file，N）功能說明：讀取錢N點的采樣值放在向量y中。（4）y=wavread（file，N1,N2）功能說明：讀取從N1到N2點的采樣值放在向量y中。 3.3 語音信號的

17、時頻域分析畫出語音信號的時域波形，再對語音信號進行頻譜分析。MATLAB提供了快速傅里葉變換算法FFT計算DFT的函數(shù)fft,其調(diào)用格式如下：Xk=fft(xn,N)9參數(shù)xn為被變換的時域序列向量，N是DFT變換區(qū)間長度，當N大于xn的長度時，fft函數(shù)自動在xn后面補零。，當N小于xn的長度時，fft函數(shù)計算xn的前N個元素，忽略其后面的元素。在本次設(shè)計中，我們利用fft對語音信號進行快速傅里葉變換，就可以得到信號的頻譜特性。程序見附錄（2）仿真圖：圖3-2 原始信號的時域頻譜圖3.4 加隨機噪聲后的時頻域分析對比本次試驗中利用MATLAB中的隨機函數(shù)(rand或randn)產(chǎn)生隨機噪聲加

18、入到語音信號中，模仿語音信號被污染，并對其頻譜分析。 matlab函數(shù)randn：產(chǎn)生正態(tài)分布的隨機數(shù)或矩陣的函數(shù)產(chǎn)生均值為0，方差 2 = 1，標準差 = 1的正態(tài)分布的隨機數(shù)或矩陣的函數(shù)。用法7： （1）Y = randn(n)：返回一個n*n的隨機項的矩陣。如果n不是個數(shù)量，將返回錯誤信息。 （2）Y = randn(m,n) 或 Y = randn(m n) ：返回一個m*n的隨機項矩陣。 （3）Y = randn(m,n,p,.) 或 Y = randn(m n p.)：產(chǎn)生隨機數(shù)組。 （4） Y = randn(size(A)：返回一個和A有同樣維數(shù)大小的隨機數(shù)組。 （5）rand

19、n：返回一個每次都變化的數(shù)量。下面一段程序?qū)崿F(xiàn)了利用randn函數(shù)把一段隨機噪音信號加入原始語音信號的信號處理過程： 程序見附錄（3）仿真圖如下：圖3-3 原始信號時域頻域圖圖3-4 原始信號相位幅度圖圖3-5 加隨機噪聲后語音信號的時域頻譜圖仿真結(jié)果分析：通過對兩張圖片的對比，很明顯可以看加噪后的語音信號時域波形比原始語音信號渾濁了許多，在時間軸上可以明顯看出00.5S的幅值增大了；通過對原始語音信號的頻譜圖與加噪后的語音信號頻譜圖的對比，也可以看出在頻率5000Hz以后的頻率幅值發(fā)生了明顯的增加5。再通過對原始語音信號的回放效果與加噪后的語音信號回放的效果的對比，人耳可以明顯辨別出兩種語音

20、信號不一樣了，加噪后的語音信號在聽覺上比原始語音信號要渾濁很多，而且還有吱吱嘎嘎的混雜音。3.5 濾波前后的時頻域比較程序見附錄（4）仿真圖 圖3-6 濾波前后比較結(jié)果分析：由圖3-6中濾波前后波形比較可看出，經(jīng)過濾波后的波形比原波形的振幅有所減小，去除了很多由于噪聲所產(chǎn)生的干擾；從濾波前后的頻譜比較可以看出經(jīng)過濾波后除了原本的聲音外，中間由于噪聲產(chǎn)生的頻譜波形已經(jīng)濾除；由圖3-4濾波前后相位比較圖可看出由于經(jīng)過濾波，相位變得稀疏；經(jīng)過MATLAB仿真，聽濾波前后的聲音，可以聽出有明顯的濾波效果。因此利用雙線性法設(shè)計的巴特沃斯濾波器已經(jīng)達到了設(shè)計的要求總結(jié)語音信號處理是語音學與數(shù)字信號處理技術(shù)

21、相結(jié)合的交叉學科，課題在這里不討論語音學，而是將語音當做一種特殊的信號，即一種“復(fù)雜向量”來看待。也就是說，課題更多的還是體現(xiàn)了數(shù)字信號處理技術(shù)。本設(shè)計采用了高效快捷的開發(fā)工具MATLAB，實現(xiàn)了語音信號的采集，對語音信號加噪聲及設(shè)計濾波器濾除噪聲的一系列工作。從頻率響應(yīng)圖中可以看出:巴特沃斯濾波器具有單調(diào)下降的幅頻特性，通帶內(nèi)是平滑的。論文初步完成了設(shè)計任務(wù)，由于本人能力有限，還存在許多不足的地方，比如濾波器的設(shè)計種類還比較單一，沒有做更多的濾波效果比較等。在以后的工作和學習中會更加努力來完善設(shè)計任務(wù)。從課題的中心來看，課題“基于MATLAB的IIR濾波器語音信號去噪”是希望將數(shù)字信號處理技

22、術(shù)應(yīng)用于某一實際領(lǐng)域，這里就是指對語音及加噪處理。作為存儲于計算機中的語音信號，其本身就是離散化了的向量，我們只需將這些離散的量提取出來，就可以對其進行處理了。這一過程的實現(xiàn)，用到了處理數(shù)字信號的強有力工具MATLAB。通過MATLAB里幾個命令函數(shù)的調(diào)用，很輕易的在實際語音與數(shù)字信號的理論之間搭了一座橋。參考文獻1高西全，丁玉美數(shù)字信號處理第3版北京：西安電子科技大學出版社，2008：55-612劉泉，闕大順數(shù)字信號處理原理與實現(xiàn)北京：電子工業(yè)出版社，2005：333張磊，畢靖，郭蓮英MATLAB實用教程北京：人民郵電出版社，2008：894張威MATLAB基礎(chǔ)與編程入門西安：西安電子科技大

23、學出版社，2006：42-455張葛祥，李娜MATLAB仿真技術(shù)與應(yīng)用北京：清華大學出版社，2003：10-136胡廣書 數(shù)字信號處理、理論、算法與實現(xiàn)M 北京：清華大學出版社，1997：33-367陳希林，肖明清一種LabWindows/CVI與MATLAB混合編程的實現(xiàn)方法J.微計算機信息，2005：578劉波MATLAB信號處理北京：電子工業(yè)出版社，2006：239施陽等MATLAB語言工具箱西安：西北工業(yè)大學出版社，1999：48-51致謝經(jīng)過幾個星期的努力和完善，系統(tǒng)終于可以正常的運行。求學歷程是艱苦的，但又是快樂的。在這里首先感謝我的課程設(shè)計指導(dǎo)教師陳海燕老師，在這段時間一直給我的

24、支持與鼓勵。認真負責的監(jiān)督我們畢業(yè)設(shè)計的進度，耐心的指導(dǎo)我們使我們能夠按時的完成任務(wù)。在這段時間學到了很多，雖然由于自身的不足沒有能夠為系統(tǒng)提出更好的解決方案。但這對我來說絕對是一個非常寶貴的歷練。從中我切身體會到了理論和現(xiàn)實的差距，只有真正動手去做才能發(fā)現(xiàn)問題。在此，我還要感謝在一起愉快的度過畢業(yè)設(shè)計生活的同學，正是由于你們的幫助和支持，我才能克服一個一個的困難和疑惑，直至本文的順利完成。特別感謝我的同學，她們對本課題做了不少工作，給予我不少的幫助。在論文即將完成之際，我的心情無法平靜，從開始進入課題到論文的順利完成，有多少可敬的師長、同學、朋友給了我無言的幫助，在這里請接受我誠摯的謝意!最

25、后我還要感謝培，謝謝你們！附錄（1）. IIR高通濾波器的仿真clf;Ft=8000;Fp=4000;Fs=3500;wp1=tan(pi*Fp/Ft); %高通到低通濾波器參數(shù)轉(zhuǎn)換ws1=tan(pi*Fs/Ft);wp=1;ws=wp1*wp/ws1;n13,wn13=cheb1ord(wp,ws,1,50,s); %求模擬的低通濾波器階數(shù)和截止頻率b13,a13=cheby1(n13,1,wn13,s); %求S域的頻率響應(yīng)的參數(shù)num,den=lp2hp(b13,a13,wn13); %將S域低通參數(shù)轉(zhuǎn)為高通的num13,den13=bilinear(num,den,0.5); %利用

26、雙線性變換實現(xiàn)S域到Z域轉(zhuǎn)換h,w=freqz(num13,den13);plot(w*21000*0.5/pi,abs(h);title( IIR高通濾波器); legend(用cheby1設(shè)計);end （2）.原始信號時域頻譜y,fs,nbits=wavread (xxx.wav);sound(y,fs,nbits); %回放語音信號N= length (y) ; %求出語音信號的長度Y=fft(y,N); %對語音信號進行傅里葉變換subplot(2,1,1);plot(y);title(原始信號時域圖 );subplot(2,1,2);plot(abs(Y);title( 原始信號頻

27、域圖)（3）.加隨機噪聲后的時頻域分析fs=22050; %語音信號采樣頻率為22050x,fs,bits=wavread( D:桌面新建文件夾functionxj.wav); sound(x,fs,bits); %播放語音信號X=fft(x,4096);magX=abs(X);angX=angle(X);y1=fft(x,1024); %對信號做1024點FFT變換f=fs*(0:511)/1024;figure(1)subplot(211);plot(magX);title(原始信號幅值);grid on;subplot(212);plot(angX);title(原始信號相位);grid

28、 on;figure(2)subplot(211);plot(x); %繪制原始語音信號的時域波形圖title(原始語音信號時域波形圖);xlabel(time n);ylabel(fuzhi n); grid on; Subplot(212); %繪制原始語音信號的頻率響應(yīng)圖 plot(f,abs(y1(1:512);title(原始語音信號頻譜圖)xlabel(Hz);ylabel(fudu); grid on; %添加隨機噪聲及添加噪聲后的時域圖和頻譜圖noise_mu=0;noise_var=0.05; x0=randn(size(x).*sqrt(noise_var)+noise_mu; x1=x+x0; ts=1/fs; % 繪制在原始信號上加隨機噪聲的信號圖 ta=(length(x)-1)/fs; t=0:ts:ta;figure(3);subplot(211); plot(t,x1); t