如何利用matlab處理音頻信號

1、Matlab處理音頻信號一、問題的提出：數(shù)字語音是信號的一種，我們處理數(shù)字語音信號，也就是對一 種信號的處理，那信號是什么呢？ 信號是傳遞信息的函數(shù)。一、問題的提出：數(shù)字語音是信號的一種，我們處理數(shù)字語音信號，也就是對一種信號的處理，那信 號是什么呢？信號是傳遞信息的函數(shù)。離散時間信號 26mdash;%26mdas市 列26mdash;%26mdash；以用圖形來表示。按信號特點的不同，信號可表示成一個或幾個獨立變量的函數(shù)。例如，圖像信號就是空間位置(二元變量)的亮度函數(shù)。一維變量可以是時間，也可以是其他參量，習慣上將其看成時間。信號有以下幾種：(1)連續(xù)時間信號：在連續(xù)時間范圍內(nèi)定義的信號

2、，但信號的幅值可以是連續(xù)數(shù)值，也可以是離散數(shù)值。當幅值為連續(xù)這一特點情況下又常稱為模擬信號。實際上 連續(xù)時間信號與模擬信號常常通用，用以說明同一信號。(2)離時間信號：時間為離散變量的信號，即獨立變量時間被量化了。而幅度仍 是連續(xù)變化的。(3)數(shù)字信號：時間離散而幅度量化的信號。語音信號是基于時間軸上的一維數(shù)字信號，在這里主要是對語音信號進行頻域上的分析。在信號分析中，頻域往往包含了更多的信息。對于頻域來說，大概有8種波形可以讓我們分析：矩形方波，鋸齒波，梯形波，臨界阻尼指數(shù)脈沖波形，三角 波，余旋波，余旋平方波，高斯波。對于各種波形，我們都可以用一種方法來分 析，就是傅立葉變換：將時域的波形

3、轉(zhuǎn)化到頻域來分析。于是，本課題就從頻域的角度對信號進行分析，并通過分析頻譜來設計出合適的濾 波器。當然,這些過程的實現(xiàn)都是在 MATLAB軟件上進行的，MATLAB軟件在數(shù) 字信號處理上發(fā)揮了相當大的優(yōu)勢。二、設計方案：利用MATLAB中的wavrea而令來讀入（采集）語音信號，將它賦值給某一向量。再將該向量看作一個普通的信號，對其進行 FFT變換實現(xiàn)頻譜分析，再依據(jù) 實際情況對它進行濾波。對于波形圖與頻譜圖（包括濾波前后的對比圖）都可以用 MATLAB畫出。我們還可以通過sound命令來對語音信號進行回放，以便在聽覺 上來感受聲音的變化。選擇設計此方案，是對數(shù)字信號處理的一次實踐。在數(shù)字信號

4、處理的課程學習過程 中，我們過多的是理論學習，幾乎沒有進行實踐方面的運用。這個課題正好是對數(shù) 字語音處理的一次有利實踐，而且語音處理也可以說是信號處理在實際應用中很大 眾化的一方面。這個方案用到的軟件也是在數(shù)字信號處理中非常通用的一個軟件26mdash;%26mdash;MATLA啾件。所以這個課題的設計過程也是一次數(shù)字信號處理在MATLAB中應用的學習過程。課題用到了較多的 MATLAB語句，而由于 課題研究范圍所限，真正與數(shù)字信號有關的命令函數(shù)卻并不多。三、主體部分：（一）、語音的錄入與打開：y,fs,bits=wavread（Blip,N1 N2）用于讀取語音，采樣值放在向量 y中，fs

5、表示采樣 頻率（Hz）, bits表示采樣位數(shù)。N1 N2表示讀取從N1點到N2點的值（若只有一 個N的點則表示讀取前N點的采樣值）。sound（x,fs,bits）用于對聲音的回放。向量y則就代表了一個信號（也即一個復雜 的26ldquo函數(shù)表達式26rdquo）也就是說可以像處理一個信號表達式一樣處理 這個聲音信號。FFT的MATLAB實現(xiàn)在MATLAB的信號處理工具箱中函數(shù)FFT和IFFT用于快速傅立葉變換和逆變 換。下面介紹這些函數(shù)。函數(shù)FFT用于序列快速傅立葉變換 函數(shù)的一種調(diào)用格式為y=fft(x)其中，x是序列，y是序列的FFT, x可以為一向量或矩陣，若x為一向量，y是x 的F

6、FT。且和x相同長度。若x為一矩陣，則y是對矩陣的每一列向量進行FFT。如果x長度是2的幕次方，函數(shù)fft執(zhí)行高速基2FFT算法；否則fft執(zhí)行一種混 合基的離散傅立葉變換算法，計算速度較慢。函數(shù)FFT的另一種調(diào)用格式為y=fft(x,N)式中，x, y意義同前，N為正整數(shù)。函數(shù)執(zhí)行N點的FFT。若x為向量且長度小于N,則函數(shù)將x補零至長度No若 向量x的長度大于N,則函數(shù)截短x使之長度為No若x為矩陣，按相同方法對x 進行處理。經(jīng)函數(shù)fft求得的序列y 一般是復序列，通常要求其幅值和相位。MATLAB提供求復數(shù)的幅值和相位函數(shù)：abs angle這些函數(shù)一般和FFT同時使用。函數(shù)abs(x用

7、于計算復向量x的幅值，函數(shù)angle(x用于計算復向量的相角，介于 和之間，以弧度表示。函數(shù)unwrap(p歸于展開弧度相位角p，當相位角絕對變化超過時，函數(shù)把它擴展 至。用MATLAB工具箱函數(shù)fft進行頻譜分析時需注意：(1)函數(shù)fft返回值y的數(shù)據(jù)結構對稱性若已知序列 x=4,3,2,6,7,8,9,0脹 X(k尸DFTx(n)。利用函數(shù)fft計算，用MATLAB編程如下：N=8;n=0:N-1;xn=4 3 2 6 7 8 9 0;XK=fft(xn)結果為：XK =39.0000-10.7782 + 6.2929i0 - 5.0000i4.7782 - 7.7071i5.00004.

8、7782 + 7.7071i0 + 5.0000i-10.7782 - 6.2929i由程序運行所得結果可見，X(k)和x(n)的維數(shù)相同，共有8個元素。X(k)的第一行 元素對應頻率值為0,第五行元素對應頻率值為 Nyquist頻率，即標準頻率為1.因 此第一行至第五行對應的標準頻率為01。而第五行至第八行對應的是負頻率，其X(k)值是以Nyquist頻率為軸對稱。(注：通常表示為Nyquist頻率外擴展，標以正值。)一般而言，對于N點的x(n)序列的FFT是N點的復數(shù)序列，其點n=N/2+1對應 Nyquist頻率，作頻譜分析時僅取序列 X(k)的前一半，即前N/2點即可。X(k)的后 一

9、半序列和前一半序列時對稱的。(2)頻率計算若N點序列x(n)(n=0,1，,N-1)是在采樣頻率 下獲得的。它的FFT也是N點序列, 即X(k)(k=0,1,2,N-1),則第k點所對應實際頻率值為f=k*f/N.(3)作FFT分析時，幅值大小與FFT選擇點數(shù)有關，但不影響分析結果。2、設計內(nèi)容：(1)下面的一段程序是語音信號在MATLAB中的最簡單表現(xiàn)，它實現(xiàn)了語音的讀入打開，以及繪出了語音信號的波形頻譜圖。x,fs,bits=wavread(ding.wav,1024 5120);sound(x,fs,bits);X=fft(x,4096);magX=abs(X);angX=angle(X

10、);subplot(221);plot(x);title媒始信號波形);subplot(222);plot(X); title(原始信號頻譜);subplot(223);plot(magX);title原始信號幅值);subplot(224);plot(angX);title原始信號相位);程序運行可以聽到聲音，得到的圖形為：100Hz,試分別繪制 N = 128(2)定點分析：已知一個語音信號，數(shù)據(jù)采樣頻率為 點DFT的幅頻圖和N = 1024點DFT幅頻圖。編程如下:x=wavread(ding.wav);sound(x);fs=100;N=128;y=fft(x,N);magy=abs(

11、y);f=(0:length(y)-1)*fs/length(y);subplot(221);plot(f,magy);xlabel(頻率(Hz) );ylabel(fi值)；title(N=128(a);gridsubplot(222);plot(f(1:N/2),magy(1:N/2);xlabel(頻率(Hz) );ylabel(fi值)；title(N=128(b);gridfs=100;N=1024;y=fft(x,N);magy=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(223);plot(f,magy);xlabel（頻率（Hz）

12、 ）;ylabel（|fi值）;title(N=1024(c);gridsubplot(224);plot(f(1:N/2),magy(1:N/2);xlabel(頻率(Hz) );ylabel(lfi值)；title(N=1024(d);grid運行結果如圖：上圖(a) (b)為N=128點幅頻譜圖，(c)、(d)為N=1024點幅頻譜圖。由于采樣頻率f =100Hz,故 Nyquist頻率為 50Hz。(a)、(c)是 0100Hz頻譜圖，(b)、(d)是 050Hz 頻譜圖。由(a域(c)可見，整個頻譜圖是以Nyquist頻率為軸對稱的。因此利用fft對 信號作頻譜分析，只要考察0Nyq

13、uist頻率(采樣頻率一半)范圍的幅頻特性。比較(a和或(b)和(d)可見，幅值大小與fft選用點數(shù)N有關，但只要點數(shù)N足夠不影響研究結果。從上圖幅頻譜可見，信號中包括15Hz和40Hz的正弦分量。(3)若信號長度T=25.6s,即抽樣后x(n)點數(shù)為T/Ts=256 ,所得頻率分辨率為Hz,以此觀察數(shù)據(jù)長度N的變化對DTFT分辨率的影響：編程如下：x,fs,bits=wavread(ding.wav);N=256;f=0:fs/N:fs/2-1/N;X=fft(x);X=abs(X);subplot(211)plot(f(45:60),X(45:60);gridxlabel(Hz),ylab

14、el(|H(ejw)|)%數(shù)據(jù)長度N擴大4倍后觀察信號頻譜N=N*4;f=0:fs/N:fs/2-1/N;X=fft(x);X=abs(X);subplot(212)plot(f(45*4:4*60),X(4*45:4*60);gridxlabel(Hz),ylabel(|H(ejw)|)結果如圖：(三)、濾波器設計：1、相關原理：設計數(shù)字濾波器的任務就是尋求一個因果穩(wěn)定的線性時不變系統(tǒng)，并使系統(tǒng)函數(shù) H(z)具有指定的頻率特性。數(shù)字濾波器從實現(xiàn)的網(wǎng)絡結構或者從單位沖激響應分類，可以分成無限長單位沖激 響應(IIR)數(shù)字濾波器和有限長單位沖激響應(FIR)數(shù)字濾波器。數(shù)字濾波器頻率響應的三個參

15、數(shù)：(1)幅度平方響應:(2)相位響應其中，相位響應(3)群時延響應IIR數(shù)字濾波器：IIR數(shù)字濾波器的系統(tǒng)函數(shù)為的有理分數(shù)，即IIR數(shù)字濾波器的逼近問題就是求解濾波器的系數(shù) 和，使得在規(guī)定的物理意義上 逼近所要求的特性的問題。如果是在 s平面上逼近，就得到模擬濾波器，如果是在 z平面上逼近，則得到數(shù)字濾波器。FIR數(shù)字濾波器：設FIR的單位脈沖響應h(n)為實數(shù)，長度為N,則其z變換和頻率響應分別為按頻域采樣定理FIR數(shù)字濾波器的傳輸函數(shù) H(z)和單位脈沖響應h(n)可由它的N 個頻域采樣值H(k)唯一確定。MATLAB中提供了幾個函數(shù)，分別用于實現(xiàn)IIR濾波器和FIR濾波器。(1)卷積函

16、數(shù)conv卷積函數(shù)conv的調(diào)用格式為c=conv(a,b)該格式可以計算兩向量a和b的卷積，可以直接用于對有限長信號采用FIR濾波器的濾波。(2)函數(shù) filter函數(shù)filter的調(diào)用格式為y=filter(b,a,x)該格式采用數(shù)字濾波器對數(shù)據(jù)進行濾波，既可以用于 IIR濾波器，也可以用于FIR 濾波器。其中向量b和a分別表示系統(tǒng)函數(shù)的分子、分母多項式的系數(shù)，若a=1,此時表示FIR濾波器，否則就是IIR濾波器。該函數(shù)是利用給出的向量b和a,對x中的數(shù)據(jù)進行濾波，結果放入向量 V。3 3)函數(shù) fftfilt函數(shù)fftfilt的調(diào)用格式為y=fftfilt(b,x)該格式是利用基于FFT

17、的重疊相加法對數(shù)據(jù)進行濾波，這種頻域濾波技術只對 FIR濾波器有效。該函數(shù)是通過向量 b描述的濾波器對x數(shù)據(jù)進行濾波。關于用butter函數(shù)求系統(tǒng)函數(shù)分子與分母系數(shù)的幾種形式。b,a=butter(N,wGhigh)設計N階高通濾波器，wc為它的3dB邊緣頻率，以 為單 位，故。b,a=butter(N,wc):當wc為具有兩個元素的矢量 wc=w1,w2時，它設計2N階帶通濾 波器，3dB通帶為，w的單位為。b,a=butter(N,wc,stop)若 wc=w1,w2,則它設計 2N 階帶阻濾波器，3dB通帶為， w的單位為。如果在這個函數(shù)輸入變元的最后，加一個變元26ldquo;s%26

18、rdquq;表示設計的是模擬濾波器。這里不作討論。為了設計任意的選項巴特沃斯濾波器，必須知道階數(shù) N和3dB邊緣頻率矢量wc。 這可以直接利用信號處理工具箱中的 buttord函數(shù)來計算。如果已知濾波器指 標，和，則調(diào)用格式為N,wc=buttord(wp,ws,Rp,As)對于不同類型的濾波器，參數(shù) wp和ws有一些限制：對于低通濾波器，wp%26lt;ws；對于高通濾波器，wp%26gt;ws；對于帶通濾波器，wp和ws分別為具有兩個元素的矢量，wp=wp1,wp2和ws=ws1,ws2,并且ws1%26lt;wp1%26lt;wp2%26lt;ws2 對于帶阻濾波器wp1%26lt;ws

19、1%26lt;ws2%26lt;wp22、設計內(nèi)容:(1)濾波器示例：在這里為了說明如何用 MATLAB來實現(xiàn)濾波，特舉出一個簡單的函數(shù)信號濾波實 例(對信號x(n尸sin( n/4)+5cos( n/2)進行濾波，信號長度為500點)，從中了解濾 波的實現(xiàn)過程。程序如下：Wn=0.2*pi;N=5;b,a=butter(N,Wn/pi);n=0:499;x=sin(pi*n/4)+5*cos(pi*n/2);X=fft(x,4096);subplot(221);plot(x);title波前信號的波形);subplot(222);plot(X);title(據(jù)波前信號的頻譜);y=filte

20、r(b,a,x);Y=fft(y,4096);subplot(223);plot(y);title就波后信號的波形);subplot(224);plot(Y);title波后信號的頻譜);在這里，是采用了 butter命令，設計出一個巴特沃斯低通濾波器，從頻譜圖中可以 很明顯的看出來。下面，也就是本課題的主要內(nèi)容，也都是運用到了butter函數(shù),以便容易的得到系統(tǒng)函數(shù)的分子與分母系數(shù)，最終以此來實現(xiàn)信號的濾波。(2) N階高通濾波器的設計(在這里，以 5階為例，其中wc為其3dB邊緣頻 率，以為單位)，程序設計如下：x=wavread(ding.wav);sound(x);N=5;wc=0.3

21、;b,a=butter(N,wc,high);X=fft(x);subplot(321);plot(x);title波前信號的波形);subplot(322);plot(X);title(據(jù)波前信號的頻譜);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR濾波后信號的波形);subplot(324);plot(Y);title(IIR濾波后信號的頻譜);z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title(FIR濾波后信號的波形);subplot(326);plot(Z);title(

22、FIR濾波后信號的頻譜);得到結果如圖:(3) N階低通濾波器的設計(在這里，同樣以5階為例，其中wc為其3dB邊緣頻率，以為單位)，程序設計如下：x=wavread(ding.wav);sound(x);N=5;wc=0.3;b,a=butter(N,wc);X=fft(x);subplot(321);plot(x);title波前信號的波形);subplot(322);plot(X);title(據(jù)波前信號的頻譜);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR濾波后信號的波形);subplot(324);plot(Y);ti

23、tle(IIR濾波后信號的頻譜);z=fffilt(b,x);Z=fft(z);subplot(325);plot(z);title(FIR濾波后信號的波形);subplot(325);plot(z);title(FIR濾波后信號的波形);得到結果如圖：(4) 2N階帶通濾波器的設計(在這里，以10階為例，其中wc為其3dB邊緣頻 率，以 為單位，wc=w1,w2,w1 wc w2 ,程序設計如下：x=wavread(ding.wav);sound(x);N=5;wc=0.3,0.6;b,a=butter(N,wc);X=fft(x);subplot(321);plot(x);title波前信

24、號的波形);subplot(322);plot(X);title(據(jù)波前信號的頻譜);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR濾波后信號的波形);subplot(324);plot(Y);title(IIR濾波后信號的頻譜);z=ff情lt(b,x);Z=fft(z);subplot(326);plot(Z);title(FIR濾波后信號的頻譜);得到結果如圖：(5) 2N階帶阻濾波器的設計(在這里，以10階為例，其中wc為其3dB邊緣頻率，以 為單位，wc=w1,w2,w1 wc w2 ,程序設計如下：x=wavread(

25、ding.wav);sound(x);N=5;wc=0.2,0.7;b,a=butter(N,wc,stop);X=fft(x);subplot(321);plot(x);title波前信號的波形);subplot(322);plot(X);title(據(jù)波前信號的頻譜);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR濾波后信號的波形);subplot(324);plot(Y);title(IIR濾波后信號的頻譜);z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title(FIR濾波后信

26、號的波形);得到結果如圖：(6)小結：以上幾種濾波，我們都可以從信號濾波前后的波形圖以及頻譜圖上看 出變化。當然，也可以用sounder數(shù)來播放濾波后的語音，從聽覺上直接感受語 音信號的變化，但由于人耳聽力的限制，有些情況下我們是很難聽出異同的。同樣，通過函數(shù)的調(diào)用，也可以將信號的頻譜進行26ldquo分離觀察26rdquo;,如顯出信號的幅值或相位。下面，通過改變系統(tǒng)函數(shù)的分子與分母系數(shù)比，來觀察 信號濾波前后的幅值與相位。并且使結果更加明顯，使人耳得以很容易的辨聽。x=wavread(ding.wav);sound(x);b=100;a=5;y=filter(b,a,x);X=fft(x,

27、4096);subplot(221);plot(x);title波前信號的波形);subplot(222);plot(abs(X);title諜波前信號的幅值);Y=fft(y,4096);subplot(223);plot(y);title(*波后信號的波形);subplot(224);plot(abs(Y);title濾波后信號的幅值);結果如圖：%26gt;%26gt; sound(y);可以聽到聲音明顯變得高亢了。從上面的波形與幅值(即幅頻)圖，也可看出，濾 波后的幅值變成了濾波前的20倍。%26gt;%26gt; figure,subplot(211);plot(angle(X);t

28、itle濾波前信號相位);subplot(212);plot(angle(Y);title(g波后信號相位);得圖：可以看到相位譜沒什么變化。(四)、界面設計：直接用M文件編寫GUI程序很繁瑣，而使用GUIDE設計工具可以大大提高工作 效率。GUIDE相當于一個控制面板，從中可以調(diào)用各種設計工具以輔助完成界面 設計任務，例如控件的創(chuàng)建和布局、控件屬性的編輯和菜單設計等。使用GUIDE設計GUI程序的一般步驟如下：1 .將所需控件從控件面板拖拽到 GUIDE的設計區(qū)域；2 .利用工具條中的工具(或相應的菜單和現(xiàn)場菜單)，快速完成界面布局；3 .設置控件的屬性。尤其是tag屬性，它是控件在程序內(nèi)部的唯一標識；4 .如果需要，打開菜單編輯器為界面添加菜單或現(xiàn)場菜單；5 .保存設計。GUIDE默認把GUI程序保存為兩個同名文件：一個是.fig文件，用 來保存窗體布局和所有控件的界面信息；一個是.m文件，該文件的初始內(nèi)容是 GUIDE自動產(chǎn)生的程序框架，其中包括了各個控件回調(diào)函數(shù)的定義。該 M文件與一般的M文件沒有本質(zhì)區(qū)別，但是鑒于它的特殊性，MATALAB把這類文件統(tǒng)稱為GUI-M文件。保存完后GUI-M文件自動在編輯調(diào)試器中打開以供編輯。6 .為每個回調(diào)函數(shù)添加代碼以實現(xiàn) G