1、數字信號處理課程設計報告 題目：語音數字信號處理與分析及Matlab實現系 別 通信工程 專業班級 學生姓名 學 號 指導教師 提交日期 摘 要 本次課程設計綜合利用數字信號處理的理論知識進行語音信號的頻譜分析，通過理論推導得出相應結論，再利用MATLAB作為編程工具進行計算機實現，從而加深對所學知識的理解，建立概念。本次課程設計要求利用MATLAB對語音信號進行分析和處理，要求學生采集語音信號后，在MATLAB軟件平臺進行頻譜分析；并對所采集的語音信號加入干擾噪聲，對加入噪聲的信號進行頻譜分析，設計合適的濾波器濾除噪聲，恢復原信號。待處理語音信號是一個在20Hz20kHz頻段的低頻信號。采用

2、了高效快捷的開發工具MATLAB，實現了語音信號的采集，對語音信號加噪聲及設計濾波器濾除噪聲的一系列工作。利用采樣原理設計了高通濾波器、低通濾波器、帶通濾波器、帶阻濾波器。同學通過查閱資料自己獲得程序進行濾波器的設計，能過得到很好的鍛煉。關鍵詞：MATLAB 濾波器 數字信號處理目錄第一章 緒論11.1 設計的目的及意義11.2 設計要求11.3 設計內容1第二章 系統方案論證32.1 設計方案分析32.2 實驗原理3第三章 信號頻譜分析63.1 原始信號及頻譜分析63.2 加入干擾噪聲后的信號及頻譜分析7第四章 數字濾波器的設計與實現114.1 高通濾波器的設計114.2 低通濾波器的設計1

3、24.3 帶通濾波器的設計154.4 帶阻濾波器的設計16第五章 課程設計總結19參考文獻20附錄I附錄II 吉林工程技術師范學院課程設計論文 第一章 緒論1.1 設計的目的及意義 綜合利用數字信號處理的理論知識進行語音信號的頻譜分析，通過理論推導得出相應結論，再利用MATLAB作為編程工具進行計算機實現，從而加深對所學知識的理解，建立概念。本設計采用了高效快捷的開發工具MATLAB，實現了語音信號的采集，對語音信號加噪聲及設計濾波器濾除噪聲的一系列工作。1.2 設計要求基本要求：本次課程設計要求利用MATLAB對語音信號進行分析和處理，要求學生采集語音信號后，在MATLAB軟件平臺進行頻譜分

4、析；并對所采集的語音信號加入干擾噪聲，對加入噪聲的信號進行頻譜分析，設計合適的濾波器濾除噪聲，恢復原信號。待處理語音信號是一個在20Hz20kHz頻段的低頻信號。1.3 設計內容選擇一個wav文件作為分析的對象，或錄制一段語音信號，對其進行頻譜分析，分別對加噪前后的語音信號進行頻譜分析，再通過不同濾波器根據信號的頻譜特點重構語音信號，選出最佳濾波方案。原理：(1)采樣定理在進行模擬/數字信號的轉換過程中，當采樣頻率fs.max大于信號中，最高頻率fmax的2倍時，即：fs.max>=2fmax,則采樣之后的數字信號完整地保留了原始信號中的信息，一般實際應用中保證采樣頻率為信號最高頻率的5

5、10倍；采樣定理又稱奈奎斯特定理： 1924年奈奎斯特(Nyquist)就推導出在理想低通信道的最高大碼元傳輸速率的公式:理想低通信道的最高大碼元傳輸速率=2W*log2 N (其中W是理想低通信道的帶寬,N是電平強度)(2)采樣頻率采樣頻率是指計算機每秒鐘采集多少個聲音樣本，是描述聲音文件的音質、音調，衡量聲卡、聲音文件的質量標準。采樣頻率越高，即采樣的間隔時間越短，則在單位時間內計算機得到的聲音樣本數據就越多，對聲音波形的表示也越精確。采樣頻率與聲音頻率之間有一定的關系，根據奎斯特理論，只有采樣頻率高于聲音信號最高頻率的兩倍時，才能把數字信號表示的聲音還原成為原來的聲音。這就是說采樣頻率是

6、衡量聲卡采集、記錄和還原聲音文件的質量標準。(3)采樣位數與采樣頻率采樣位數即采樣值或取樣值，用來衡量聲音波動變化的參數，是指聲卡在采集和播放聲音文件時所使用數字聲音信號的二進制位數。采樣頻率是指錄音設備在一秒鐘內對聲音信號的采樣次數，采樣頻率越高聲音的還原就越真實越自然。采樣位數和采樣率對于音頻接口來說是最為重要的兩個指標，也是選擇音頻接口的兩個重要標準。無論采樣頻率如何，理論上來說采樣的位數決定了音頻數據最大的力度范圍。每增加一個采樣位數相當于力度范圍增加了6dB。采樣位數越多則捕捉到的信號越精確。對于采樣率來說你可以想象它類似于一個照相機，44.1kHz意味著音頻流進入計算機時計算機每秒

7、會對其拍照達441000次。顯然采樣率越高，計算機攝取的圖片越多，對于原始音頻的還原也越加精確。第二章 系統方案論證2.1 設計方案分析信號的傅立葉表示在信號的分析與處理中起著重要的作用。因為對于線性系統來說，可以很方便地確定其對正弦或復指數和的響應，所以傅立葉分析方法能完善地解決許多信號分析和處理問題。另外，傅立葉表示使信號的某些特性變得更明顯，因此，它能更深入地說明信號的各項紅物理現象。 由于語音信號是隨著時間變化的，通常認為，語音是一個受準周期脈沖或隨機噪聲源激勵的線性系統的輸出。輸出頻譜是聲道系統頻率響應與激勵源頻譜的乘積。聲道系統的頻率響應及激勵源都是隨時間變化的，因此一般標準的傅立

8、葉表示雖然適用于周期及平穩隨機信號的表示，但不能直接用于語音信號。由于語音信號可以認為在短時間內，近似不變，因而可以采用短時分析法。本實驗要求掌握傅里葉分析原理，會利用已學的知識，編寫程序估計短時譜、倒譜，畫出語譜圖，并分析實驗結果，在此基礎上，借助頻域分析方法所求得的參數分析語音信號的基音周期或共振峰。2.2 實驗原理1、短時傅立葉變換由于語音信號是短時平穩的隨機信號，某一語音信號幀的短時傅立葉變換的定義為： （2-1）其中w(n-m)是實窗口函數序列，n表示某一語音信號幀。令n-m=k'，則得到 （2-2）于是可以得到 （2-3）假定 （2-4）則可以得到 （2-5）同樣，不同的窗

9、口函數，將得到不同的傅立葉變換式的結果。由上式可見，短時傅立葉變換有兩個變量：n和，所以它既是時序n的離散函數，又是角頻率的連續函數。與離散傅立葉變換逼近傅立葉變換一樣，如令=2k/N，則得離散的短時傅立葉吧如下： (2-6)2、語譜圖水平方向是時間軸，垂直方向是頻率軸，圖上的灰度條紋代表各個時刻的語音短時譜。語譜圖反映了語音信號的動態頻率特性，在語音分析中具有重要的實用價值。被成為可視語言。語譜圖的時間分辨率和頻率分辨率是由窗函數的特性決定的。時間分辨率高，可以看出時間波形的每個周期及共振峰隨時間的變化，但頻率分辨率低，不足以分辨由于激勵所形成的細微結構，稱為寬帶語譜圖；而窄帶語譜圖正好與之

10、相反。寬帶語譜圖可以獲得較高的時間分辨率，反映頻譜的快速時變過程；窄帶語譜圖可以獲得較高的頻率分辨率，反映頻譜的精細結構。兩者相結合，可以提供帶兩與語音特性相關的信息。語譜圖上因其不同的灰度，形成不同的紋路，稱之為“聲紋”。聲紋因人而異，因此可以在司法、安全等場合得到應用。3、復倒譜和倒譜復倒譜是x(n)的Z變換取對數后的逆Z變換，其表達式如下: （2-7）倒譜c(n)定義為x(n)取Z變換后的幅度對數的逆Z變換，即 （2-8）在時域上，語音產生模型實際上是一個激勵信號與聲道沖激響應的卷積。對于濁音，激勵信號可以由周期脈沖序列表示；對于清音，激勵信號可以由隨機噪聲序列表示。聲道系統相當于參數緩

11、慢變化的零極點線性濾波器。這樣經過同態處理后，語音信號的復倒譜，激勵信號的復倒譜，聲道系統的復倒譜之間滿足下面的關系： （2-9）由于倒譜對應于復倒譜的偶部，因此倒譜與復倒譜具有同樣的特點，很容易知道語音信號的倒譜，激勵信號的倒譜以及聲道系統的倒譜之間滿足下面關系： （2-10）濁音信號的倒譜中存在著峰值，它的出現位置等于該語音段的基音周期，而清音的倒譜中則不存在峰值。利用這個特點我們可以進行清濁音的判斷，并且可以估計濁音的基音周期。第三章 信號頻譜分析3.1 原始信號及頻譜分析 選擇一個wav文件作為分析的對象，可以利用Windows下的錄音機或其他軟件，錄制一段自己的話音，在MATLAB中

12、，y,fs,bits=wavread('Blip',N1 N2);用于讀取語音，采樣值放在向量y中，fs表示采樣頻率(Hz)，bits表示采樣位數。N1 N2表示讀取的值從N1點到N2點的值。 sound(y); 用于對聲音的回放。向量y則就代表了一個信號，也即一個復雜的“函數表達式”，也可以說像處理一個信號的表達式一樣處理這個聲音信號。 下面是語音信號在MATLAB中的語言程序，它實現了語音的讀入與打開，并繪出了語音信號時域波形，然后對語音信號進行頻譜分析。在MATLAB中，可以利用函數fft對信號進行快速傅里葉變化，得到信號的頻譜特性。 在頻譜特性中分析最大值的位置（可能有

13、幾個），它代表的頻率和時域的采樣時間有關，相鄰的兩點之間的距離為。其中，N是離散傅里葉變換用的點數，是采樣的時間，前面在讀取 wav文件時得到了采樣頻率。 既然知道了該聲波的頻譜，按頻率就可以反演它的時域值，利用以上分析的主要峰值來重構聲波。由于沒有考慮相位和其他的頻譜分量，所以波形和原來的波形相差甚大，但大體的頻率是沒錯的。如圖3-1所示。圖3-1 原始信號頻譜圖3.2 加入干擾噪聲后的信號及頻譜分析 wav語音信號加噪聲 在MATLAB軟件平臺下，給原始的語音信號疊加上噪聲，噪聲類型分為如下幾種：（1）單頻噪色（正弦干擾）；（2）高斯隨機噪聲。繪出加噪聲后的語音信號時域和頻譜圖，在視覺上與

14、原始語音信號圖形對比，也可通過Windows播放軟件從聽覺上進行對比，分析并 體會含噪語音信號頻譜和時域波形的改變。本實驗采用正弦干擾。 clc; clear; fs=22050; %語音信號采樣頻率為22050 x,fs,bits=wavread('C:Documents and SettingsAdministrator桌面語音音頻.wav'); %讀取語音信號的數據，賦給變量x y1=fft(x,4096); %對信號做4096點FFT變換 f=fs*(0:511)/4096; t=(0:length(x)-1)/22050; x1=0.05*sin(2*pi*10000

15、*t)' x2=x+x1; sound(x2,fs,bits); figure(1) subplot(2,1,1) plot(x) %做原始語音信號的時域圖形 title('原語音信號時域圖') subplot(2,1,2) plot(x2) %做原始語音信號的時域圖形 title('加高斯噪聲后語音信號時域圖') xlabel('time n'); ylabel('fudu'); y2=fft(x2,4096); figure(2) subplot(2,1,1) plot(abs(y1) title('原始語音信

16、號頻譜'); xlabel('Hz'); ylabel('fudu'); subplot(2,1,2) plot(abs(y2) title('加噪語音信號頻譜'); xlabel('Hz'); ylabel('fudu'); axis(0 4500 0 300); wavwrite(x2,fs,'C:Documents and SettingsAdministrator桌面語音加噪.wav'); 調試分析過程描述：I.語音信號的分析及處理方法a.語音的錄入與打開在MATLAB中，y,fs,

17、bits=wavread('Blip',N1 N2);用于讀取語音，采樣值放在向量y中，fs表示采樣頻率(Hz)，bits表示采樣位數。N1 N2表示讀取從N1點到N2點的值（若只有一個N的點則表示讀取前N點的采樣值）。 sound(x,fs,bits); 用于對聲音的回放。向量y則就代表了一個信號（也即一個復雜的“函數表達式”）也就是說可以像處理一個信號表達式一樣處理這個聲音信號。b.時域信號的FFT分析FFT即為快速傅氏變換，是離散傅氏變換的快速算法，它是根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。在MATLAB的信號處理工具箱中函數FFT

18、和IFFT用于快速傅立葉變換和逆變換。函數FFT用于序列快速傅立葉變換，其調用格式為y=fft(x)，其中，x是序列，y是序列的FFT，x可以為一向量或矩陣，若x為一向量，y是x的FFT且和x相同長度；若x為一矩陣，則y是對矩陣的每一列向量進行FFT。如果x長度是2的冪次方，函數fft執行高速基2FFT算法，否則fft執行一種混合基的離散傅立葉變換算法，計算速度較慢。函數FFT的另一種調用格式為y=fft(x,N)，式中，x，y意義同前，N為正整數。函數執行N點的FFT，若x為向量且長度小于N，則函數將x補零至長度N；若向量x的長度大于N，則函數截短x使之長度為N；若x 為矩陣，按相同方法對x

19、進行處理。II.數字濾波器設計原理數字濾波是數字信號分析中最重要的組成部分之一，與模擬濾波相比，它具有精度和穩定性高、系統函數容易改變、靈活性強、便于大規模集成和可實現多維濾波等優點。在信號的過濾、檢測和參數的估計等方面，經典數字濾波器是使用最廣泛的一種線性系統。數字濾波器的作用是利用離散時間系統的特性對輸入信號波形(或頻譜)進行加工處理，或者說利用數字方法按預定的要求對信號進行變換。數字濾波器的設計步驟：不論是IIR濾波器還是FIR濾波器的設計都包括三個步驟：a. 按照實際任務的要求，確定濾波器的性能指標。b. 用一個因果、穩定的離散線性時不變系統的系統函數去逼近這一性能指標。根據不同的要求

20、可以用IIR系統函數，也可以用FIR系統函數去逼近。c. 利用有限精度算法實現系統函數，包括結構選擇、字長選擇等。III.IIR濾波器與FIR濾波器的性能比較有限沖擊響應FIR:Finite Impulse response，IIR:Infinite Impulse response，無限沖擊響應從性能上來說，IIR濾波器傳輸函數的極點可位于單位圓內的任何地方，因此可用較低的階數獲得高的選擇性，所用的存貯單元少，所以經濟而效率高。但是這個高效率是以相位的非線性為代價的。選擇性越好，則相位非線性越嚴重。相反，FIR濾波器卻可以得到嚴格的線性相位，然而由于FIR濾波器傳輸函數的極點固定在原點，所以

21、只能用較高的階數達到高的選擇性；對于同樣的濾波器設計指標，FIR濾波器所要求的階數可以比IIR濾波器高510倍，結果，成本較高，信號延時也較大；如果按相同的選擇性和相同的線性要求來說，則IIR濾波器就必須加全通網絡進行相位較正，同樣要大增加濾波器的節數和復雜性。整體來看，IIR濾波器達到同樣效果階數少，延遲小，但是有穩定性問題，非線性相位；FIR濾波器沒有穩定性問題，線性相位，但階數多，延遲大。如圖3-2所示。圖3-2 加入噪聲后的信號頻譜圖第四章 數字濾波器的設計與實現4.1 高通濾波器的設計（1）模擬低通濾波器原型設計。 在此單元模塊主要依據設計參數完成模擬低通濾波器的原型設計。 模擬原型

22、濾波器指的是截止頻率為1的濾波器。此處的模擬低通濾波器原型以巴特沃斯（Butterworth）低通濾波器為模型進行逼近。MATLAB信號處理工具箱提供Butterworth模擬低通濾波器原型設計函數buttap,函數調用形式為： z,p,k=buttap(N) 式中，N為butterworth濾波器階數；z,p,k分別為濾波器的零點、極點和增益。又Butterworth低通濾波器的幅度函數只由階數N控制，假定階數為1，則語句描述為： z,p,k=buttap1; 再求解模擬低通濾波器的分子分母系數。最后利用MATLAB工具箱提供的模擬濾波器()aHs的頻率響應函數freqs和相關繪圖函數求出其

23、傳遞函數圖形。其語句描述為： H,w=freqs(B,A,n); magH2=(abs(H).2; hold on; plot(w,magH2); end xlabel('w/wc');ylabel('|H(jw)|2'); title('Butterworth模擬原型濾波器'); 至此，可得到模擬低通濾波器的原型，該濾波器的截止頻率為為1。 （2）頻率轉換，得到模擬高通濾波器。 此單元模塊可以利用MATLAB工具箱提供的模擬頻率變換函數求得。即利用函數lp2hp即可由模擬低通濾波器得到模擬高通濾波器。語句描述為： Bt,At=lp2hp(B,A

24、,wc); 其中wc為題設所要求的模擬高通濾波器的阻帶截止角頻率，即 wc=2*pi*f=2*pi*200=400*pi 至此即得到符合設計參數要求的阻帶頻率為200Hz的模擬高通濾波器。 最后利用MATLAB工具箱提供的模擬濾波器()aHs的頻率響應函數freqs和相關繪圖函數求出各響應曲線，以便對濾波器的性能有更為直觀的了解和認識。語句描述為： h,w=freqs(Bt,At); subplot(); plot(w,20*log10(abs(h); grid; xlabel(w/pi); ylabel(頻率響應/dB); title('模擬高通濾波器');如圖4-1所示。圖

25、4-1 高通濾波器頻譜圖4.2 低通濾波器的設計巴特沃斯低通濾波，對加入高斯隨機噪聲和正弦噪聲的語音信號進行濾波。用雙線性變換法設計了巴特沃斯數字低通IIR濾波器對兩加噪語音信號進行濾波，并繪制了巴特沃斯低通濾波器的幅度圖和兩加噪語音信號濾波前后的時域圖和頻譜圖。 clear all; fb = 1000; fc = 1200; fs = 22050; wp=0.1*pi; ws=0.4*pi; Rp=1; Rs=15;Fs=22050;Ts=1/Fs; wp1=2/Ts*tan(wp/2); %將模擬指標轉換成數字指標 ws1=2/Ts*tan(ws/2); N,Wn=buttord(wp1

26、,ws1,Rp,Rs,'s'); %選擇濾波器的最小階數Z,P,K=buttap(N); %創建butterworth模擬濾波器 Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn); bd,ad=bilinear(b,a,Fs); %用雙線性變換法實現模擬濾波器到數字濾波器的轉換 h,w=freqz(bd,ad); figure(1) subplot(111); plot(w*fs/(2*pi),abs(h) grid; title('濾波器的性能分析'); pause; figure(2) x,fs,bits=wavread(

27、'C:Documents and SettingsAdministrator桌面語音加噪.wav'); n=length(x); f=fs*(0:(n/2-1)/n; X=fft(x); z=filter(bd,ad,x); subplot(211); plot(x); title('原始信號的波形'); subplot(212); plot(z); title('濾波后信號的波形');pause;figure(3) sound(z,fs,bits);subplot(211); plot(f,abs(X(1:n/2); title('原始

28、信號的頻譜'); xlabel('Hz'); Z=fft(z); subplot(212); plot(f,abs(Z(1:n/2); title('濾波后的信號頻譜');xlabel('Hz'); wavwrite(z,fs,'C:Documents and SettingsAdministrator桌面語音巴濾.wav');如圖4-2所示。圖4-2 低通濾波器頻譜圖4.3 帶通濾波器的設計帶通濾波器設計步驟： 1、根據需求選擇合適的低通濾波器原型 2、把帶通濾波器帶寬作為低通濾波器的截止頻率，根據抑

29、制點的頻率距離帶通濾波器中心頻點距離的兩倍作為需要抑制的頻率，換算抑制頻率與截止頻率的比值，得出的值，然后根據m值選擇低通濾波器的原型參數值。 濾波器的時域特性  ：任何信號通過濾波器都會產生時延。Bessel filter是特殊的濾波器在于對于通帶內的所有頻率而言，引入的時延都是恒定的。這就意味著相對于輸入，輸出信號的相位變化與工作的頻率是成比例的。而其他類型的濾波器（如Butterworth, Chebyshev,inverse Chebyshev,andCauser）在輸出信號中引入的相位變化與頻率不成比例。相位隨頻率變化的速率

30、稱之為群延遲（group delay）。群延遲隨濾波器級數的增加而增加。如圖4-3所示。圖4-3帶通濾波器頻譜圖  4.4 帶阻濾波器的設計帶阻濾波器廣泛應用于無線通信系統中，用來抑制高功率發射機的雜散輸出以及非線性功放或帶通濾波器產生的寄生通帶等。將帶通濾波器的各種拓撲結構應用于帶阻濾波器是可行的，而提取各種拓撲結構濾波器的等效電路參數成為設計濾波器的關鍵。function y=bands(x,f1,f3,fsl,fsh,rp,rs,Fs) %帶阻濾波 %使用注意事項：通帶或阻帶的截止頻率與采樣率的選取范圍是不能超過采樣率的一半 %即，f1,f3,fs1,fsh,的值小于

31、 Fs/2%x:需要帶通濾波的序列 % f 1：通帶左邊界 % f 3：通帶右邊界% fs1：衰減截止左邊界% fsh：衰變截止右邊界 %rp：邊帶區衰減DB數設置 %rs：截止區衰減DB數設置 %FS：序列x的采樣頻率 % f1=300;f3=500;%通帶截止頻率上下限 % fsl=200;fsh=600;%阻帶截止頻率上下限% rp=0.1;rs=30;%通帶邊衰減DB值和阻帶邊衰減DB值 % Fs=2000;%采樣率 %wp1=2*pi*f1/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsl/Fs; wsh=2*pi*fsh/Fs; wp=wp1 wp3; ws=wsl w

32、sh; % % 設計切比雪夫濾波器； n,wn=cheb1ord(ws/pi,wp/pi,rp,rs); bz1,az1=cheby1(n,rp,wp/pi,'stop'); %查看設計濾波器的曲線 h,w=freqz(bz1,az1,256,Fs); h=20*log10(abs(h); figure;plot(w,h);title('所設計濾波器的通帶曲線');grid on; y=filter(bz1,az1,x); end 使用例子 %帶阻濾波器測試 fs=1000; t=(1:fs)/fs; y=sin(2*pi*100*t)+sin(2*pi*150*t)+sin(2*pi*200*t); figure;hua_fft(y,fs,1); z=bands(y,110,190,140,160,0.1,30,fs); figure;hua_fft(z,fs,1); 如圖4-4所示。 圖4-4 帶阻濾波器頻譜圖第五章 課程設計總結本設計圓滿的完