matlab音頻處理

1、matlab處理音頻信號(hào)一、           問題的提出：數(shù)字語音是信號(hào)的一種，我們處理數(shù)字語音信號(hào)，也就是對(duì)一種信號(hào)的處理，那信號(hào)是什么呢？信號(hào)是傳遞信息的函數(shù)。離散時(shí)間信號(hào)序列可以用圖形來表示。按信號(hào)特點(diǎn)的不同，信號(hào)可表示成一個(gè)或幾個(gè)獨(dú)立變量的函數(shù)。例如，圖像信號(hào)就是空間位置（二元變量）的亮度函數(shù)。一維變量可以是時(shí)間，也可以是其他參量，習(xí)慣上將其看成時(shí)間。信號(hào)有以下幾種：（1）連續(xù)時(shí)間信號(hào)：在連續(xù)時(shí)間范圍內(nèi)定義的信號(hào)，但信號(hào)的幅值可以是連續(xù)數(shù)值，也可以是離散數(shù)值。當(dāng)幅值為連續(xù)這一特點(diǎn)情況

2、下又常稱為模擬信號(hào)。實(shí)際上連續(xù)時(shí)間信號(hào)與模擬信號(hào)常常通用，用以說明同一信號(hào)。（2）離時(shí)間信號(hào)：時(shí)間為離散變量的信號(hào)，即獨(dú)立變量時(shí)間被量化了。而幅度仍是連續(xù)變化的。（3）數(shù)字信號(hào)：時(shí)間離散而幅度量化的信號(hào)。語音信號(hào)是基于時(shí)間軸上的一維數(shù)字信號(hào)，在這里主要是對(duì)語音信號(hào)進(jìn)行頻域上的分析。在信號(hào)分析中，頻域往往包含了更多的信息。對(duì)于頻域來說，大概有8種波形可以讓我們分析：矩形方波，鋸齒波，梯形波，臨界阻尼指數(shù)脈沖波形，三角波，余旋波，余旋平方波，高斯波。對(duì)于各種波形，我們都可以用一種方法來分析，就是傅立葉變換：將時(shí)域的波形轉(zhuǎn)化到頻域來分析。于是，本課題就從頻域的角度對(duì)信號(hào)進(jìn)行分析，并通過分析頻譜來設(shè)計(jì)

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

4、程中，我們過多的是理論學(xué)習(xí)，幾乎沒有進(jìn)行實(shí)踐方面的運(yùn)用。這個(gè)課題正好是對(duì)數(shù)字語音處理的一次有利實(shí)踐，而且語音處理也可以說是信號(hào)處理在實(shí)際應(yīng)用中很大眾化的一方面。這個(gè)方案用到的軟件也是在數(shù)字信號(hào)處理中非常通用的一個(gè)軟件MATLAB軟件。所以這個(gè)課題的設(shè)計(jì)過程也是一次數(shù)字信號(hào)處理在MATLAB中應(yīng)用的學(xué)習(xí)過程。課題用到了較多的MATLAB語句，而由于課題研究范圍所限，真正與數(shù)字信號(hào)有關(guān)的命令函數(shù)卻并不多。三、       主體部分：（一）、語音的錄入與打開：y,fs,bits=wavread('Blip',N1 N2);用

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

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

7、x進(jìn)行處理。經(jīng)函數(shù)fft求得的序列y一般是復(fù)序列，通常要求其幅值和相位。MATLAB提供求復(fù)數(shù)的幅值和相位函數(shù)：abs，angle，這些函數(shù)一般和FFT同時(shí)使用。函數(shù)abs(x)用于計(jì)算復(fù)向量x的幅值，函數(shù)angle(x)用于計(jì)算復(fù)向量的相角，介于 和 之間，以弧度表示。函數(shù)unwrap(p)用于展開弧度相位角p ,當(dāng)相位角絕對(duì)變化超過 時(shí)，函數(shù)把它擴(kuò)展至 。用MATLAB工具箱函數(shù)fft進(jìn)行頻譜分析時(shí)需注意：（1）    函數(shù)fft返回值y的數(shù)據(jù)結(jié)構(gòu)對(duì)稱性若已知序列x=4,3,2,6,7,8,9,0,求X(k)=DFTx(n)。利用函數(shù)fft計(jì)算，用MATLAB

8、編程如下：N=8;n=0:N-1;xn=4 3 2 6 7 8 9 0'XK=fft(xn)結(jié)果為：XK =39.0000         -10.7782 + 6.2929i        0 - 5.0000i   4.7782 - 7.7071i   5.0000          &

9、#160; 4.7782 + 7.7071i        0 + 5.0000i-10.7782 - 6.2929i由程序運(yùn)行所得結(jié)果可見，X(k)和x(n)的維數(shù)相同，共有8個(gè)元素。X(k)的第一行元素對(duì)應(yīng)頻率值為0，第五行元素對(duì)應(yīng)頻率值為Nyquist頻率，即標(biāo)準(zhǔn)頻率為1.因此第一行至第五行對(duì)應(yīng)的標(biāo)準(zhǔn)頻率為01。而第五行至第八行對(duì)應(yīng)的是負(fù)頻率，其X(k)值是以Nyquist頻率為軸對(duì)稱。（注：通常表示為Nyquist頻率外擴(kuò)展，標(biāo)以正值。）一般而言，對(duì)于N點(diǎn)的x(n)序列的FFT是N點(diǎn)的復(fù)數(shù)序列，其點(diǎn)n=N/2+1對(duì)應(yīng)

10、Nyquist頻率，作頻譜分析時(shí)僅取序列X(k)的前一半，即前N/2點(diǎn)即可。X(k)的后一半序列和前一半序列時(shí)對(duì)稱的。（2）    頻率計(jì)算    若N點(diǎn)序列x(n)(n=0,1,N-1)是在采樣頻率 下獲得的。它的FFT也是N點(diǎn)序列，即X(k)(k=0,1,2,N-1)，則第k點(diǎn)所對(duì)應(yīng)實(shí)際頻率值為f=k*f /N.（3）    作FFT分析時(shí)，幅值大小與FFT選擇點(diǎn)數(shù)有關(guān)，但不影響分析結(jié)果。2、設(shè)計(jì)內(nèi)容：（1）下面的一段程序是語音信號(hào)在MATLAB中的最簡單表現(xiàn)，它實(shí)現(xiàn)了語音的讀入打開，以及繪出了語音信

11、號(hào)的波形頻譜圖。     x,fs,bits=wavread('ding.wav',1024 5120);   sound(x,fs,bits);   X=fft(x,4096);magX=abs(X);angX=angle(X);   subplot(221);plot(x);title('原始信號(hào)波形');subplot(222);plot(X); title('原始信號(hào)頻譜');subplot(223);plot(magX);title(&#

12、39;原始信號(hào)幅值');subplot(224);plot(angX);title('原始信號(hào)相位');程序運(yùn)行可以聽到聲音，得到的圖形為：（2）定點(diǎn)分析：已知一個(gè)語音信號(hào)，數(shù)據(jù)采樣頻率為100Hz，試分別繪制N128點(diǎn)DFT的幅頻圖和N1024點(diǎn)DFT幅頻圖。     編程如下：x=wavread('ding.wav');   sound(x);fs=100;N=128;y=fft(x,N);magy=abs(y);f=(0:length(y)-1)'*fs/length(y);sub

13、plot(221);plot(f,magy);xlabel('頻率（Hz）');ylabel('幅值');title('N=128(a)');gridsubplot(222);plot(f(1:N/2),magy(1:N/2);xlabel('頻率（Hz）');ylabel('幅值');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);

14、plot(f,magy);xlabel('頻率（Hz）');ylabel('幅值');title('N=1024(c)');gridsubplot(224);plot(f(1:N/2),magy(1:N/2);xlabel('頻率（Hz）');ylabel('幅值');title('N=1024(d)');grid運(yùn)行結(jié)果如圖：上圖(a)、(b)為N=128點(diǎn)幅頻譜圖，(c)、(d)為N=1024點(diǎn)幅頻譜圖。由于采樣頻率f =100Hz，故Nyquist頻率為 50Hz。(a)、(c)是0100Hz

15、頻譜圖，(b)、(d)是050Hz頻譜圖。由(a)或(c)可見，整個(gè)頻譜圖是以Nyquist頻率為軸對(duì)稱的。因此利用fft對(duì)信號(hào)作頻譜分析，只要考察0Nyquist頻率（采樣頻率一半）范圍的幅頻特性。比較(a)和(c)或(b)和(d)可見，幅值大小與fft選用點(diǎn)數(shù)N有關(guān)，但只要點(diǎn)數(shù)N足夠不影響研究結(jié)果。從上圖幅頻譜可見，信號(hào)中包括15Hz和40Hz的正弦分量。（3）若信號(hào)長度T=25.6s，即抽樣后x(n)點(diǎn)數(shù)為T/Ts=256，所得頻率分辨率為 Hz，以此觀察數(shù)據(jù)長度N的變化對(duì)DTFT分辨率的影響：編程如下：x,fs,bits=wavread('ding.wav');N=25

16、6;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'),ylabel('|H(ejw)|')%數(shù)據(jù)長度N擴(kuò)大4倍后觀察信號(hào)頻譜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)|')結(jié)果如圖：（三）、濾波器設(shè)計(jì)：1、相關(guān)原理：設(shè)

17、計(jì)數(shù)字濾波器的任務(wù)就是尋求一個(gè)因果穩(wěn)定的線性時(shí)不變系統(tǒng)，并使系統(tǒng)函數(shù)H(z)具有指定的頻率特性。數(shù)字濾波器從實(shí)現(xiàn)的網(wǎng)絡(luò)結(jié)構(gòu)或者從單位沖激響應(yīng)分類，可以分成無限長單位沖激響應(yīng)（IIR）數(shù)字濾波器和有限長單位沖激響應(yīng)（FIR）數(shù)字濾波器。數(shù)字濾波器頻率響應(yīng)的三個(gè)參數(shù)：（1）    幅度平方響應(yīng)：（2）    相位響應(yīng)其中，相位響應(yīng)             （3）    群時(shí)延響應(yīng)I

18、IR數(shù)字濾波器：IIR數(shù)字濾波器的系統(tǒng)函數(shù)為 的有理分?jǐn)?shù)，即IIR數(shù)字濾波器的逼近問題就是求解濾波器的系數(shù) 和 ，使得在規(guī)定的物理意義上逼近所要求的特性的問題。如果是在s平面上逼近，就得到模擬濾波器，如果是在z平面上逼近，則得到數(shù)字濾波器。FIR數(shù)字濾波器：設(shè)FIR的單位脈沖響應(yīng)h(n)為實(shí)數(shù)，長度為N，則其z變換和頻率響應(yīng)分別為按頻域采樣定理FIR數(shù)字濾波器的傳輸函數(shù)H(z)和單位脈沖響應(yīng)h(n)可由它的N個(gè)頻域采樣值H(k)唯一確定。MATLAB中提供了幾個(gè)函數(shù)，分別用于實(shí)現(xiàn)IIR濾波器和FIR濾波器。（1）卷積函數(shù)conv卷積函數(shù)conv的調(diào)用格式為   &

19、#160;       c=conv(a,b)該格式可以計(jì)算兩向量a和b的卷積，可以直接用于對(duì)有限長信號(hào)采用FIR濾波器的濾波。（2）函數(shù)filter函數(shù)filter的調(diào)用格式為          y=filter(b,a,x)該格式采用數(shù)字濾波器對(duì)數(shù)據(jù)進(jìn)行濾波，既可以用于IIR濾波器，也可以用于FIR濾波器。其中向量b和a分別表示系統(tǒng)函數(shù)的分子、分母多項(xiàng)式的系數(shù)，若a1，此時(shí)表示FIR濾波器，否則就是IIR濾波器。該函數(shù)是利用給出的向量b和a，

20、對(duì)x中的數(shù)據(jù)進(jìn)行濾波，結(jié)果放入向量y。（3）函數(shù)fftfilt函數(shù)fftfilt的調(diào)用格式為         y=fftfilt(b,x)該格式是利用基于FFT的重疊相加法對(duì)數(shù)據(jù)進(jìn)行濾波，這種頻域?yàn)V波技術(shù)只對(duì)FIR濾波器有效。該函數(shù)是通過向量b描述的濾波器對(duì)x數(shù)據(jù)進(jìn)行濾波。關(guān)于用butter函數(shù)求系統(tǒng)函數(shù)分子與分母系數(shù)的幾種形式。b,a=butter(N,wc,'high'):設(shè)計(jì)N階高通濾波器，wc為它的3dB邊緣頻率，以 為單位，故 。b,a=butter(N,wc):當(dāng)wc為具有兩個(gè)元素的矢量

21、wc=w1,w2時(shí)，它設(shè)計(jì)2N階帶通濾波器，3dB通帶為 ，w的單位為 。b,a=butter(N,wc,'stop'):若wc=w1,w2，則它設(shè)計(jì)2N階帶阻濾波器，3dB通帶為 ，w的單位為 。如果在這個(gè)函數(shù)輸入變?cè)淖詈?，加一個(gè)變?cè)皊”，表示設(shè)計(jì)的是模擬濾波器。這里不作討論。為了設(shè)計(jì)任意的選項(xiàng)巴特沃斯濾波器，必須知道階數(shù)N和3dB邊緣頻率矢量wc。這可以直接利用信號(hào)處理工具箱中的buttord函數(shù)來計(jì)算。如果已知濾波器指標(biāo) ， ， 和 ，則調(diào)用格式為N,wc=buttord(wp,ws,Rp,As)對(duì)于不同類型的濾波器，參數(shù)wp和ws有一些限制：對(duì)于低通濾波器，wp&l

22、t;ws；對(duì)于高通濾波器，wp>ws；對(duì)于帶通濾波器，wp和ws分別為具有兩個(gè)元素的矢量，wp=wp1,wp2和ws=ws1,ws2，并且 ws1<wp1<wp2<ws2；對(duì)于帶阻濾波器wp1<ws1<ws2<wp2。2、設(shè)計(jì)內(nèi)容：（1）濾波器示例：在這里為了說明如何用MATLAB來實(shí)現(xiàn)濾波，特舉出一個(gè)簡單的函數(shù)信號(hào)濾波實(shí)例（對(duì)信號(hào)x(n)=sin( n/4)+5cos( n/2)進(jìn)行濾波，信號(hào)長度為500點(diǎn)），從中了解濾波的實(shí)現(xiàn)過程。程序如下：Wn=0.2*pi;N=5;b,a=butter(N,Wn/pi);n=0:499;x=sin(pi*n/4

23、)+5*cos(pi*n/2);X=fft(x,4096);subplot(221);plot(x);title('濾波前信號(hào)的波形');subplot(222);plot(X);title('濾波前信號(hào)的頻譜');y=filter(b,a,x);Y=fft(y,4096);subplot(223);plot(y);title('濾波后信號(hào)的波形');subplot(224);plot(Y);title('濾波后信號(hào)的頻譜');在這里，是采用了butter命令，設(shè)計(jì)出一個(gè)巴特沃斯低通濾波器，從頻譜圖中可以很明顯的看出來。下面，也就

24、是本課題的主要內(nèi)容，也都是運(yùn)用到了butter函數(shù)，以便容易的得到系統(tǒng)函數(shù)的分子與分母系數(shù)，最終以此來實(shí)現(xiàn)信號(hào)的濾波。（2）N階高通濾波器的設(shè)計(jì)（在這里，以5階為例，其中wc為其3dB邊緣頻率，以 為單位），程序設(shè)計(jì)如下：x=wavread('ding.wav');   sound(x);N=5;wc=0.3;b,a=butter(N,wc,'high');X=fft(x);subplot(321);plot(x);title('濾波前信號(hào)的波形');subplot(322);plot(X);title('濾波前信號(hào)的

25、頻譜');y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title('IIR濾波后信號(hào)的波形');subplot(324);plot(Y);title('IIR濾波后信號(hào)的頻譜');z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title('FIR濾波后信號(hào)的波形');subplot(326);plot(Z);title('FIR濾波后信號(hào)的頻譜');得到結(jié)果如圖：（3）N階低通濾波器的設(shè)計(jì)（在這里，同樣以5階為例，其中wc為其3dB

26、邊緣頻率，以 為單位），程序設(shè)計(jì)如下：x=wavread('ding.wav');   sound(x);N=5;wc=0.3;b,a=butter(N,wc);X=fft(x);subplot(321);plot(x);title('濾波前信號(hào)的波形');subplot(322);plot(X);title('濾波前信號(hào)的頻譜');y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title('IIR濾波后信號(hào)的波形');subplot(324);plot(Y);ti

27、tle('IIR濾波后信號(hào)的頻譜');z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title('FIR濾波后信號(hào)的波形');subplot(326);plot(Z);title('FIR濾波后信號(hào)的頻譜');得到結(jié)果如圖：（4）2N階帶通濾波器的設(shè)計(jì)（在這里，以10階為例，其中wc為其3dB邊緣頻率，以 為單位，wc=w1,w2,w1 wc w2），程序設(shè)計(jì)如下：x=wavread('ding.wav');   sound(x);N=5;wc=0.3,0.6;b,a

28、=butter(N,wc);X=fft(x);subplot(321);plot(x);title('濾波前信號(hào)的波形');subplot(322);plot(X);title('濾波前信號(hào)的頻譜');y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title('IIR濾波后信號(hào)的波形');subplot(324);plot(Y);title('IIR濾波后信號(hào)的頻譜');z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title('FI

29、R濾波后信號(hào)的波形');subplot(326);plot(Z);title('FIR濾波后信號(hào)的頻譜');得到結(jié)果如圖：（5）2N階帶阻濾波器的設(shè)計(jì)（在這里，以10階為例，其中wc為其3dB邊緣頻率，以 為單位，wc=w1,w2,w1 wc w2），程序設(shè)計(jì)如下：x=wavread('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('濾波前信號(hào)的波形');

30、subplot(322);plot(X);title('濾波前信號(hào)的頻譜');y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title('IIR濾波后信號(hào)的波形');subplot(324);plot(Y);title('IIR濾波后信號(hào)的頻譜');z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title('FIR濾波后信號(hào)的波形');subplot(326);plot(Z);title('FIR濾波后信號(hào)的頻譜');得到結(jié)果

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

32、r(b,a,x);X=fft(x,4096);subplot(221);plot(x);title('濾波前信號(hào)的波形');subplot(222);plot(abs(X);title('濾波前信號(hào)的幅值');Y=fft(y,4096);subplot(223);plot(y);title('濾波后信號(hào)的波形');subplot(224);plot(abs(Y);title('濾波后信號(hào)的幅值');結(jié)果如圖：>> sound(y);可以聽到聲音明顯變得高亢了。從上面的波形與幅值（即幅頻）圖，也可看出，濾波后的幅值變成了

33、濾波前的20倍。>> figure,subplot(211);plot(angle(X);title('濾波前信號(hào)相位');subplot(212);plot(angle(Y);title('濾波后信號(hào)相位');得圖：可以看到相位譜沒什么變化。（四）、界面設(shè)計(jì)：直接用M文件編寫GUI程序很繁瑣，而使用GUIDE設(shè)計(jì)工具可以大大提高工作效率。GUIDE相當(dāng)于一個(gè)控制面板，從中可以調(diào)用各種設(shè)計(jì)工具以輔助完成界面設(shè)計(jì)任務(wù)，例如控件的創(chuàng)建和布局、控件屬性的編輯和菜單設(shè)計(jì)等。使用GUIDE設(shè)計(jì)GUI程序的一般步驟如下：1.      將所需控件從控件面板拖拽到GUIDE的設(shè)計(jì)區(qū)域；2.      利用工具條中的工具（或相應(yīng)的菜單和現(xiàn)場(chǎng)菜單），快速完成界面布局；3.      設(shè)置控件的屬性。尤其是tag屬性，它是控件在程序內(nèi)部的唯一標(biāo)識(shí)；4.      如果需要，打開菜單編輯器為界面添加菜單或現(xiàn)場(chǎng)菜單；5.      保存設(shè)計(jì)。GUIDE默認(rèn)把GUI程序