數字信號處理期末綜合實驗報告.doc

《數字信號處理期末綜合實驗報告.doc》由會員分享，可在線閱讀，更多相關《數字信號處理期末綜合實驗報告.doc（21頁珍藏版）》請在裝配圖網上搜索。

數字信號處理綜合實驗報告 實驗題目：基于Matlab的語音信號去噪及仿真 專業(yè)名稱： 學 號： 姓 名： 日 期： 報告內容： 一、 實驗原理 1、 去噪的原理 1.1 采樣定理 在進行模擬/數字信號的轉換過程中，當采樣頻率fs.max大于信號中，最高頻率fmax的2倍時，即：fs.max>=2fmax,則采樣之后的數字信號完整地保留了原始信號中的信息，一般實際應用中保證采樣頻率為信號最高頻率的5～10倍；采樣定理又稱奈奎斯特定理。 1924年奈奎斯特(Nyquist)就推導出在理想低通信道的最高大碼元傳輸速率的公式: 理想低通信道的最高大碼元傳輸速率=2W*log2 N (其中W是理想低通信道的帶寬,N是電平強度)為什么把采樣頻率設為8kHz?在數字通信中，根據采樣定理, 最小采樣頻率為語音信號最高頻率的2倍 頻帶為F的連續(xù)信號 f(t)可用一系列離散的采樣值f(t1),f(t1Δt)，f(t12Δt)，...來表示,只要這些采樣點的時間間隔Δt≤1/2F，便可根據各采樣值完全恢復原來的信號f(t)。 這是時域采樣定理的一種表述方式。 時域采樣定理的另一種表述方式是：當時間信號函數f(t)的最高頻率分量為fM時,f(t)的值可由一系列采樣間隔小于或等于1/2fM的采樣值來確定,即采樣點的重復頻率f≥2fM。圖為模擬信號和采樣樣本的示意圖。 時域采樣定理是采樣誤差理論、隨機變量采樣理論和多變量采樣理論的基礎。對于時間上受限制的連續(xù)信號f(t)（即當│t│>T時,f(t)=0,這里T=T2-T1是信號的持續(xù)時間），若其頻譜為F（ω）,則可在頻域上用一系列離散的采樣值 （1-1） 采樣值來表示,只要這些采樣點的頻率間隔 （1-2） 。 1.2 采樣頻率 采樣頻率，也稱為采樣速度或者采樣率，定義了每秒從連續(xù)信號中提取并組成離散信號的采樣個數，它用赫茲（Hz）來表示。采樣頻率的倒數是采樣周期或者叫作采樣時間，它是采樣之間的時間間隔。通俗的講采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質、音調，衡量聲卡、聲音文件的質量標準。 采樣頻率只能用于周期性采樣的采樣器，對于非周期性采樣的采樣器沒有規(guī)則限制。 采樣頻率的常用的表示符號是 fs。 通俗的講采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質、音調，衡量聲卡、聲音文件的質量標準。采樣頻率越高，即采樣的間隔時間越短，則在單位時間內計算機得到的聲音樣本數據就越多，對聲音波形的表示也越精確。采樣頻率與聲音頻率之間有一定的關系，根據采樣定理，只有采樣頻率高于聲音信號最高頻率的兩倍時，才能把數字信號表示的聲音還原成為原來的聲音。這就是說采樣頻率是衡量聲卡采集、記錄和還原聲音文件的質量標準。 采樣位數和采樣率對于音頻接口來說是最為重要的兩個指標，也是選擇音頻接口的兩個重要標準。無論采樣頻率如何，理論上來說采樣的位數決定了音頻數據最大的力度范圍。每增加一個采樣位數相當于力度范圍增加了6dB。采樣位數越多則捕捉到的信號越精確。對于采樣率來說你可以想象它類似于一個照相機，44.1kHz意味著音頻流進入計算機時計算機每秒會對其拍照達441000次。顯然采樣率越高，計算機攝取的圖片越多，對于原始音頻的還原也越加精確 2、 數字濾波器設計的基本原理 1.1 FIR數字濾波器的設計 FIR：有限脈沖響應濾波器。有限說明其脈沖響應是有限的。與IIR相比，它具有線性相位、容易設計的優(yōu)點。這也就說明，IIR濾波器具有相位不線性，不容易設計的缺點。而另一方面，IIR卻擁有FIR所不具有的缺點，那就是設計同樣參數的濾波器，FIR比IIR需要更多的參數。這也就說明，要增加DSP的計算量。DSP需要更多的計算時間，對DSP的實時性有影響。FIR濾波器的設計比較簡單，就是要設計一個數字濾波器去逼近一個理想的低通濾波器。通常這個理想的低通濾波器在頻域上是一個矩形窗。根據傅里葉變換我們可以知道，此函數在時域上是一個采樣函數。通常此函數的表達式為： sa（n）＝sin（n）/n（1-3） 但是這個采樣序列是無限的，計算機是無法對它進行計算的。故我們需要對此采樣函數進行截斷處理。也就是加一個窗函數。就是傳說中的加窗。也就是把這個時域采樣序列去乘一個窗函數，就把這個無限的時域采樣序列截成了有限個序列值。但是加窗后對此采樣序列的頻域也產生了影響：此時的頻域便不在是一個理想的矩形窗，而是成了一個有過渡帶，阻帶有波動的低通濾波器。通常根據所加的窗函數的不同，對采樣信號加窗后，在頻域所得的低通濾波器的阻帶衰減也不同。通常我們就是根據此阻帶衰減去選擇一個合適的窗函數。如矩形窗、漢寧窗、漢明窗、BLACKMAN窗、凱撒窗等。 1.2 IIR數字濾波器的設計 對于數字高通、帶通濾波器的設計，通用方法為雙線性變換法?？梢越柚谀M濾波器的頻率轉換設計一個所需類型的過渡模擬濾波器，再經過雙線性變換將其轉換策劃那個所需的數字濾波器。具體設計步驟如下： （1）確定所需類型數字濾波器的技術指標。 （2）將所需類型數字濾波器的邊界頻率轉換成相應的模擬濾波器的邊界頻率，轉換公式為Ω=2/T tan(0.5ω) （1-4） (3)將相應類型的模擬濾波器技術指標轉換成模擬低通濾波器技術指標。 （4）設計模擬低通濾波器。 （5）通過頻率變換將模擬低通轉換成相應類型的過渡模擬濾波器。 （6）采用雙線性變換法將相應類型的過渡模擬濾波器轉換成所需類型的數字濾波器。 我們知道，脈沖響應不變法的主要缺點是會產生頻譜混疊現象，使數字濾波器的頻響偏離模擬濾波器的頻響特性。為了克服之一缺點，可以采用雙線性變換法。 二、 實驗步驟 其大概流程框圖可如下表示：（圖2-1） 語音信號采集 效果顯示、對比 語音信號錄入 語音信號變換 信號加噪 語音信號濾波 圖2-1 1. 濾波器的設計。 利用窗函數法設計FIR濾波器的步驟。如下： （1）根據對阻帶衰減及過渡帶的指標要求，選擇串窗數類型（矩形窗、三角窗、漢寧窗、哈明窗、凱塞窗等），并估計窗口長度N。先按照阻帶衰減選擇窗函數類型。原則是在保證阻帶衰減滿足要求的情況下，盡量選擇主瓣的窗函數。 （2）構造希望逼近的頻率響應函數。 （3）計算h(n).。 （4）加窗得到設計結果。 接下來，我們根據語音信號的特點給出有關濾波器的技術指標： 低通濾波器的性能指標： fp=1000Hz，fc=1200Hz，As=50db ,Ap=1dB 高通濾波器的性能指標： fp=3500Hz，fc=4000Hz，As=50dB，Ap=1dB 在Matlab中，可以利用函數fir1設計FIR濾波器，利用Matlab中的函數freqz畫出各步步器的頻率響應。 MATLAB信號處理工具箱函數cheblap,cheblord和cheeby1是切比雪夫I型濾波器設計函數。我們用到的是cheeby1函數，其調用格式如下： [B,A]=cheby1(N,Rp,wpo,’ftypr’) [B,A]=cheby1(N,Rp,wpo,’ftypr’,’s’) 利用模擬濾波器設計IIR數字低通濾波器的步驟：如下： （1）確定數字低通濾波器的技術指標：通帶邊界頻率、通帶最大衰減，阻帶截止頻率、阻帶最小衰減。 （2）將數字低通濾波器的技術指標轉換成相應的模擬低通濾波器的技術指標。 （3）按照模擬低通濾波器的技術指標設計及過渡模擬低通濾波器。 （4）用雙線性變換法，模擬濾波器系統(tǒng)函數轉換成數字低通濾波器系統(tǒng)函數。 MATLAB信號處理工具箱函數cheblap,cheblord和cheeby1是切比雪夫I型濾波器設計函數。我們用到的是cheeby1函數，其調用格式如下： [B,A]=cheby1(N,Rp,wpo,’ftypr’) [B,A]=cheby1(N,Rp,wpo,’ftypr’,’s’) 函數butter,cheby1和ellip設計IIR濾波器時都是默認的雙線性變換法，所以 在設計濾波器時只需要代入相應的實現函數即可。 2. 語音信號的錄制。 單擊自己的電腦開始程序，選擇所有程序，接著選擇附件，再選擇錄音。自己錄入語音信號，然后保存在MATLAB文件夾里面，命名為“chenghaijie1.wav”。 3. 在MATLAB平臺上讀入語音信號、繪制頻譜圖并回放原始語音信號。 利用MATLAB中的wavread命令來讀入（采集）語音信號，將它賦值給某一向量。 [y,fs,bits]=wavread( [N1 N2]);用于讀取語音，采樣值放在向量y中，fs表示采樣頻率(Hz)，bits表示采樣位數。[N1 N2]表示讀取從N1點到N2點的值（若只有一個N的點則表示讀取前N點的采樣值。 4. 利用matlab函數randn編程加入一段隨機噪音信號，再利用設計的FIR和IIR濾波器去噪，分別繪制去噪后的頻譜圖、回放語音信號與原始信號的頻譜圖、原始語音信號比較，并且比較兩種濾波器的優(yōu)缺點和得出語音信號的頻段。 三、 實驗結果 下面我們將給出設計FIR數字濾波器的主要程序和圖像。 FIR低通濾波器程序見附錄1； FIR低通濾波器圖像：（圖2—2） 圖2—2 FIR低通濾波器 FIR高通濾波程序見附錄2； FIR高通濾波圖像：（圖2—3） 圖2—3 FIR高通濾波器 下面我們將給出IIR數字濾波器的主要程序。 IIR低通濾波器程序見附錄3； IIR低通濾波器圖像：（圖2—4） 圖2—4 IIR低通濾波器 IIR濾波器高通程序見附錄4； IIR濾波器高通圖像：（圖2—5） 圖2—5 IIR高通濾波器 原始語音信號的回放及頻譜分析程序：見附錄5； 原始語音信號及其頻譜分析圖像、運行程序可以聽到原始信號的回放：（圖2—6） 圖2—6 原始語音信號及其頻譜分析圖像 利用randn函數把一段隨機噪音信號加入原始語音信號的信號處理過程程序：見附錄6； 加噪后語音信號的時域波形、頻譜圖：（圖2—7） （圖2—7）加噪后語音信號的時域波形、頻譜圖 利用FIR低通濾波器法去噪程序：見附錄7； FIR濾波器去噪前與去噪后語音信號的時域波形、頻譜圖：（圖2—8） （圖2—8）FIR濾波器去噪前與去噪后語音信號的時域波形、頻譜圖 利用IIR低通濾波器法去噪程序：見附錄8； IIR濾波器去噪前與去噪后語音信號的時域波形、頻譜圖：（圖2—9） （圖2—9）IIR濾波器去噪前與去噪后語音信號的時域波形、頻譜圖 四、 實驗分析 1. 由（圖2—8）FIR濾波器去噪前與去噪后語音信號的時域波形、頻譜圖和（圖2—9）IIR濾波器去噪前與去噪后語音信號的時域波形、頻譜圖可看出原始語音信號和加噪語音信號時域波形和頻譜圖的區(qū)別。加噪后的語音信號的時域波形比原始語音信號要模糊得多，頻譜圖則是在頻率5000Hz以后出現了明顯的變化。 再通過濾波前的信號波形和頻譜圖的對比，可以明顯看出濾波后的波形開始變得清晰了，有點接近原始信號的波形圖了。濾波后信號的頻譜圖也在5000Hz以后開始逐漸接近原始語音信號的頻譜圖。 再從對語音信號的回放，人耳可以明顯辨別出加噪后的語音信號比較渾濁，還有很明顯嘎吱嘎吱的雜音在里面。濾波后，語音信號較加噪后的信號有了明顯的改善，基本可以聽清楚了，而且雜音也沒有那么強烈，但是聲音依然沒有原始語音信號那么清晰脆耳。 2. 結合去噪后的頻譜圖對比可知FIR濾波器和IIR濾波器的優(yōu)缺點 IIR數字濾波器采用遞歸型結構，即結構上帶有反饋環(huán)路。IIR濾波器運算結構通常由延時、乘以系數和相加等基本運算組成，可以組合成直接型、正準型、級聯(lián)型、并聯(lián)型四種結構形式，都具有反饋回路。由于運算中的舍入處理，使誤差不斷累積，有時會產生微弱的寄生振蕩。 （1）IIR數字濾波器的相位特性不好控制，對相位要求較高時，需加相位校準網絡。FIR濾波器則要求較低。 （2）IIR濾波器運算誤差大，有可能出現極限環(huán)振蕩，FIR相比之下運算誤差較小，不會出現極限環(huán)振蕩。 （3）IIR幅頻特性精度很高，不是線性相位的，可以應用于對相位信息不敏感的音頻信號上； （4）與FIR濾波器的設計不同，IIR濾波器設計時的階數不是由設計者指定，而是根據設計者輸入的各個濾波器參數（截止頻率、通帶濾紋、阻帶衰減等），由軟件設計出滿足這些參數的最低濾波器階數。在MATLAB下設計不同類型IIR濾波器均有與之對應的函數用于階數的選擇。 （5）IIR單位響應為無限脈沖序列FIR單位響應為有限的 （6）FIR幅頻特性精度較之于iir低，但是線性相位，就是不同頻率分量的信號經過FIR濾波器后他們的時間差不變。這是很好的性質。 （7）IIR濾波器有噪聲反饋，而且噪聲較大，FIR濾波器噪聲較小。 FIR幅頻特性精度較之于IIR低，但是線性相位，就是不同頻率分量的信號經過FIR濾波器后他們的時間差不變，這是很好的性質。 3. 通過對語音信號做頻譜分析和通過比較加噪前后，語音的頻譜和語音回放，能明顯的感覺到加入噪聲后回放的聲音與原始的語音信號有很大的不同,前者隨較尖銳的干擾嘯叫聲。從含噪語音信號的頻譜圖中可以看出含噪聲的語音信號頻譜,在整個頻域范圍內分是布均勻。其實，這正是干擾所造成的。通過濾波前后的對比，低通濾波后效果最好，高通濾波后的效果最差。由此可見，語音信號主要分布在低頻段，而噪聲主要分布在高頻段。 五、 實驗總結 基于Matlab的語音信號去噪及仿真課題的特色在于它將語音信號看作一個向量，于是就把語音數字化了。那么，就可以完全利用數字信號處理的知識來解決語音及加噪處理問題。我們可以像給一般信號做頻譜分析一樣，來對語音信號做頻譜分析，也可以較容易的用數字濾波器來對語音進行濾波處理。通過比較加噪前后，語音的頻譜和語音回放，能明顯的感覺到加入噪聲后回放的聲音與原始的語音信號有很大的不同,前者隨較尖銳的干擾嘯叫聲。從含噪語音信號的頻譜圖中可以看出含噪聲的語音信號頻譜,在整個頻域范圍內分是布均勻。其實，這正是干擾所造成的。通過兩種濾波器濾波前后的對比，可知FIR濾波器和IIR濾波器的優(yōu)缺點，低通濾波后效果最好，高通濾波后的效果最差。由此可見，語音信號主要分布在低頻段，而噪聲主要分布在高頻段。 六、 參考文獻 [1] Sanjit K.Mitra，數字信號處理實驗指導書，電子工業(yè)出版社，2005 年 1月 [2] 胡航，語音信號處理，哈爾濱工業(yè)大學出版社，2000 年 5 月 [3］姚天任.數字語音處理[M].武漢：華中科技大學出版社，2005 年 3 月 七、 附錄 附錄 1 Ft=8000;%設置抽樣頻率為8000Hz Fp=1000;%設置通帶邊界為1000Hz Fs=1200;%設置阻帶邊界為1200Hz wp=2*Fp/Ft;%計算歸一化通帶邊界角頻率 ws=2*Fs/Ft;%計算歸一化阻帶邊界角頻率 rp=1;%設置通帶紋波為1dB rs=50;%設置阻帶衰減為50dB p=1-10.^(-rp/20);%通帶阻帶波紋 s=10.^(-rs/20); fpts=[wp ws]; mag=[1 0]; dev=[p s]; [n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%kaiserord求階數截止頻率 b21=fir1(n21,wn21,Kaiser(n21+1,beta));%由fir1設計濾波器 [h,w]=freqz(b21,1);%得到頻率響應 plot(w/pi,abs(h));%繪制FIR低通濾波器 title(FIR低通濾波器);%標題 附錄2 Ft=8001;%設置抽樣頻率為8001Hz Fp=4000;%設置通帶邊界為4000Hz Fs=3500;%設置阻帶邊界為3500Hz wp=2*Fp/Ft;%計算歸一化通帶邊界角頻率 ws=2*Fs/Ft;%計算歸一化阻帶邊界角頻率 rp=1;%設置通帶紋波為1dB rs=50;%設置阻帶衰減為50dB p=1-10.^(-rp/20);%通帶阻帶波紋 s=10.^(-rs/20); fpts=[ws wp]; mag=[0 1]; dev=[p s]; [n23,wn23,beta,ftype]=kaiserord(fpts,mag,dev);%kaiserord求階數截止頻率 b23=fir1(n23,wn23,high,Kaiser(n23+1,beta));%由fir1設計濾波器 [h,w]=freqz(b23,1);%得到頻率響應 plot(w*12000*0.5/pi,abs(h));%繪制FIR高通濾波器 title(FIR高通濾波器);%標題 axis([3000 6000 0 1.2]);%確定橫縱坐標的范圍 附錄3 Ft=8000;%設置抽樣頻率為8000Hz Fp=1000;%設置通帶邊界為1000Hz Fs=1200;%設置阻帶邊界為1200Hz wp=2*pi*Fp/Ft;%計算歸一化通帶邊界角頻率 ws=2*pi*Fs/Ft;%計算歸一化阻帶邊界角頻率 fp=2*Ft*tan(wp/2);%計算通帶邊界頻率 fs=2*Fs*tan(wp/2);%計算阻帶邊界頻率 [n11,wn11]=buttord(wp,ws,1,50,s);%求低通濾波器的階數和截止頻率 [b11,a11]=butter(n11,wn11,s);%利用butter求S域的頻率響應的參數 [num11,den11]=bilinear(b11,a11,0.5);%雙線性變換實現S域到Z域的變換 [h,w]=freqz(num11,den11);%根據參數求出頻率響應 plot(w*8000*0.5/pi,abs(h));%繪制IIR低通濾波器 title(IIR低通濾波器);%標題 legend(用butter設計);%添加圖例的標注 grid;%使繪制的圖有網格 附錄4 Ft=8001;%設置抽樣頻率為8001Hz Fp=4000;%設置通帶邊界為4000Hz Fs=3500;%設置阻帶邊界為3500Hz wp1=tan(pi*Fp/Ft);%通帶到低通濾波器參數轉換 ws1=tan(pi*Fs/Ft);%阻帶到低通濾波器參數轉換 wp=1;%計算歸一化通帶邊界角頻率 ws=wp1*wp/ws1;%計算歸一化阻帶邊界角頻率 [n13,wn13]=cheb1ord(wp,ws,1,50,s);%求模擬的低通濾波器階數和截止頻 [b13,a13]=cheby1(n13,1,wn13,s);%利用cheby1求S域的頻率響應的參數 [num,den]=lp2hp(b13,a13,wn13); %將S域低通參數轉為高通的 [num13,den13]=bilinear(num,den,0.5); %利用雙線性變換實現S域到Z域轉 [h,w]=freqz(num13,den13);%得到頻率響應 plot(w*21000*0.5/pi,abs(h));%繪制IIR高通濾波器 title(IIR高通濾波器);%標題 legend(用cheby1設計);%添加圖例的標注 附錄5 [x,fs,bits]=wavread(chenghaijie1.wav);%利用MATLAB中的wavread命令來讀入（采集）語音信號，將它賦值給某一向量 sound(x,fs,bits);%播放語音信號 X=fft(x,100000);%對信號做50000點FFT變換 magX=abs(X);%求幅度 angX=angle(X);%求相位 plot(x);title(原始信號波形);%繪制原始信號的波形 plot(X); title(原始語音信號采樣后的頻譜圖)%繪制原始語音信號采樣后的頻譜圖 附錄6 [y,fs,bits]=wavread(chenghaijie1.wav);%利用MATLAB中的wavread命令來讀入（采集）語音信號，將它賦值給某一向量 sound(y,fs)%播放語音信號 n=length(y)%計算向量y的長度 y_p=fft(y,n);%對向量y進行fft變換 f=fs*(0:n/2-1)/n;%對應點的頻率 figure(1)%建立第一幅圖 subplot(2,1,1);%分配繪制的位置為第一塊位置 plot(y);%繪制原始語音信號采樣后時域波形 title(原始語音信號采樣后時域波形);%標題 xlabel(時間軸)%標注X軸 ylabel(幅值 A)%標注Y軸 subplot(2,1,2);%分配繪制的位置為第二塊位置 plot(f,abs(y_p(1:n/2)));%繪制原始語音信號采樣后的頻譜圖 title(原始語音信號采樣后的頻譜圖); xlabel(頻率Hz);%標注X軸 ylabel(頻率幅值);%標注Y軸 L=length(y)%計算向量y的長度 noise=0.1*randn(L,2);%設置加入的隨機噪聲函數 y_z=y+noise;%將原始信號與隨機噪聲函數相疊加 sound(y_z,fs)%試聽加入隨機噪聲函數的信號 n=length(y);%計算向量y的長度 y_zp=fft(y_z,L);%對加入正弦噪聲向量y進行fft變換 f=fs*(0:L/2-1)/L;%計算對應點的頻率 figure(2)%建立第二幅圖 subplot(2,1,1);%分配繪制的位置為第一塊位置 plot(y_z);%繪制加噪語音信號時域波形 title(加噪語音信號時域波形);%標題 xlabel(時間軸)%標注X軸 ylabel(幅值 A)%標注Y軸 subplot(2,1,2);%分配繪制的位置為第二塊位置 plot(f,abs(y_zp(1:L/2)));%繪制加噪語音信號頻譜圖 title(加噪語音信號頻譜圖);%標題 xlabel(頻率Hz);%標注X軸 ylabel(頻率幅值);%標注Y軸 附錄7 [y,fs,bits]=wavread(chenghaijie1.wav);%利用MATLAB中的wavread命令來讀入（采集）語音信號，將它賦值給某一向量 sound(y,fs)%播放語音信號 n=length(y)%計算向量y的長度 y_p=fft(y,n);%對向量y進行fft變換 f=fs*(0:n/2-1)/n;%對應點的頻率 figure(1)%建立第一幅圖 subplot(2,1,1);%分配繪制的位置為第一塊位置 plot(y);%繪制原始語音信號采樣后時域波形 title(原始語音信號采樣后時域波形);%標題 xlabel(時間軸)%標注X軸 ylabel(幅值 A)%標注Y軸 subplot(2,1,2);%分配繪制的位置為第二塊位置 plot(f,abs(y_p(1:n/2)));%繪制原始語音信號采樣后的頻譜圖 title(原始語音信號采樣后的頻譜圖); xlabel(頻率Hz);%標注X軸 ylabel(頻率幅值);%標注Y軸 L=length(y)%計算向量y的長度 noise=0.1*randn(L,2);%設置加入的隨機噪聲函數 y_z=y+noise;%將原始信號與隨機噪聲函數相疊加 sound(y_z,fs)%試聽加入隨機噪聲函數的信號 n=length(y);%計算向量y的長度 y_zp=fft(y_z,L);%對加入正弦噪聲向量y進行fft變換 f=fs*(0:L/2-1)/L;%計算對應點的頻率 figure(2)%建立第二幅圖 subplot(2,1,1);%分配繪制的位置為第一塊位置 plot(y_z);%繪制加噪語音信號時域波形 title(加噪語音信號時域波形);%標題 xlabel(時間軸)%標注X軸 ylabel(幅值 A)%標注Y軸 subplot(2,1,2);%分配繪制的位置為第二塊位置 plot(f,abs(y_zp(1:L/2)));%繪制加噪語音信號頻譜圖 title(加噪語音信號頻譜圖);%標題 xlabel(頻率Hz);%標注X軸 ylabel(頻率幅值);%標注Y軸 Ft=5000;%設置抽樣頻率為5000Hz Fp=1000;%設置通帶邊界為1000Hz Fs=1200;%設置阻帶邊界為1200Hz wp=2*Fp/Ft;%計算歸一化通帶邊界角頻率 ws=2*Fs/Ft;%計算歸一化阻帶邊界角頻率 rp=1;%設置通帶紋波為1dB rs=50;%設置阻帶衰減為50dB p=1-10.^(-rp/20);%通帶阻帶波紋 s=10.^(-rs/20); fpts=[wp ws]; mag=[1 0]; dev=[p s]; [n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%kaiserord求階數截止頻率 b21=fir1(n21,wn21,Kaiser(n21+1,beta));%由fir1設計濾波器 [h,w]=freqz(b21,1);%得到頻率響應 plot(w/pi,abs(h));%繪制FIR低通濾波器 title(FIR低通濾波器);%標題 x=fftfilt(b21,y_z);%重疊相加法實現線性卷積的計算 X=fft(x,n);%計算x的fft figure(3);%建立第三幅圖 subplot(2,2,1);plot(f,abs(y_zp(1:n/2)));%繪制濾波前信號的頻譜 title(濾波前信號的頻譜);%標題 subplot(2,2,2);plot(f,abs(X(1:n/2)));%繪制濾波后信號的頻譜 title(濾波后信號的頻譜);%標題 subplot(2,2,3);plot(y_z);%繪制濾波前信號的時域波形 title(濾波前信號的時域波形)%標題 subplot(2,2,4);plot(x);%繪制濾波后信號的時域波形 title(濾波后信號的時域波形)%標題 sound(x,fs,bits)%重新聽取濾噪后的信號 附錄8 Ft=8000;%設置抽樣頻率為8000Hz Fp=1000;%設置通帶邊界為1000Hz Fs=1200;%設置阻帶邊界為1200Hz wp=2*pi*Fp/Ft;%計算歸一化通帶邊界角頻率 ws=2*pi*Fs/Ft;%計算歸一化阻帶邊界角頻率 fp=2*Ft*tan(wp/2);%計算通帶邊界頻率 fs=2*Fs*tan(wp/2);%計算阻帶邊界頻率 [n11,wn11]=buttord(wp,ws,1,50,s);%求低通濾波器的階數和截止頻率 [b11,a11]=butter(n11,wn11,s);%利用butter求S域的頻率響應的參數 [num11,den11]=bilinear(b11,a11,0.5);%雙線性變換實現S域到Z域的變換 [h,w]=freqz(num11,den11);%根據參數求出頻率響應 plot(w*8000*0.5/pi,abs(h));%繪制IIR低通濾波器 title(IIR低通濾波器);%標題 legend(用butter設計);%添加圖例的標注 grid;%使繪制的圖有網格 [y,fs,nbits]=wavread (chenghaijie1.wav);%利用MATLAB中的wavread命令來讀入（采集）語音信號，將它賦值給某一向量 n = length (y) ;%求出語音信號的長度 noise=0.01*randn(n,2);%隨機函數產生噪聲 s=y+noise;%語音信號加入噪聲 S=fft(s); %傅里葉變換 z11=filter(num11,den11,s);%濾噪聲 sound(z11);%聽取濾噪聲后的信號 m11=fft(z11);%求濾波后的信號 subplot(2,2,1);%分配繪制的位置為第一塊位置 plot(abs(S),g);%繪制濾波前信號的頻譜 title(濾波前信號的頻譜);%標題 grid;%網格 subplot(2,2,2);%分配繪制的位置為第二塊位置 plot(abs(m11),r);%繪制濾波后信號的頻譜 title(濾波后信號的頻譜);%標題 grid;%網格 subplot(2,2,3);%分配繪制的位置為第三塊位置 plot(s);%繪制濾波前信號的波形 title(濾波前信號的波形);%標題 grid;%網格 subplot(2,2,4);%分配繪制的位置為第四塊位置 plot(z11);%繪制濾波后的信號波形 title(濾波后的信號波形);%標題