matlab處理音頻信號

一、 問題的提出：

數字語音是信號的一種，我們處理數字語音信號，也就是對一種信號的處理，那信號是什么呢？

信號是傳遞信息的函數。離散時間信號%26mdash;%26mdash;序列%26mdash;%26mdash;可以用圖形來表示。

按信號特點的不同，信號可表示成一個或幾個獨立變量的函數。例如，圖像信號就是空間位置（二元變量）的亮度函數。一維變量可以是時間，也可以是其他參量，習慣上將其看成時間。信號有以下幾種：

（1）連續時間信號：在連續時間范圍內定義的信號，但信號的幅值可以是連續數值，也可以是離散數值。當幅值為連續這一特點情況下又常稱為模擬信號。實際上連續時間信號與模擬信號常常通用，用以說明同一信號。

（2）離時間信號：時間為離散變量的信號，即獨立變量時間被量化了。而幅度仍是連續變化的。

（3）數字信號：時間離散而幅度量化的信號。

語音信號是基于時間軸上的一維數字信號，在這里主要是對語音信號進行頻域上的分析。在信號分析中，頻域往往包含了更多的信息。對于頻域來說，大概有8種波形可以讓我們分析：矩形方波，鋸齒波，梯形波，臨界阻尼指數脈沖波形，三角波，余旋波，余旋平方波，高斯波。對于各種波形，我們都可以用一種方法來分析，就是傅立葉變換：將時域的波形轉化到頻域來分析。

于是，本課題就從頻域的角度對信號進行分析，并通過分析頻譜來設計出合適的濾波器。當然，這些過程的實現都是在MATLAB軟件上進行的，MATLAB軟件在數字信號處理上發揮了相當大的優勢。

二、 設計方案：

利用MATLAB中的wavread命令來讀入（采集）語音信號，將它賦值給某一向量。再將該向量看作一個普通的信號，對其進行FFT變換實現頻譜分析，再依據實際情況對它進行濾波。對于波形圖與頻譜圖（包括濾波前后的對比圖）都可以用 MATLAB畫出。我們還可以通過sound命令來對語音信號進行回放，以便在聽覺上來感受聲音的變化。

選擇設計此方案，是對數字信號處理的一次實踐。在數字信號處理的課程學習過程中，我們過多的是理論學習，幾乎沒有進行實踐方面的運用。這個課題正好是對數字語音處理的一次有利實踐，而且語音處理也可以說是信號處理在實際應用中很大眾化的一方面。

這個方案用到的軟件也是在數字信號處理中非常通用的一個軟件%26mdash;%26mdash;MATLAB軟件。所以這個課題的設計過程也是一次數字信號處理在MATLAB中應用的學習過程。課題用到了較多的MATLAB語句，而由于課題研究范圍所限，真正與數字信號有關的命令函數卻并不多。

三、 主體部分：

（一）、語音的錄入與打開：

[y,fs,bits]=wavread('Blip',[N1 N2]);用于讀取語音，采樣值放在向量y中，fs表示采樣頻率(Hz)，bits表示采樣位數。[N1 N2]表示讀取從N1點到N2點的值（若只有一個N的點則表示讀取前N點的采樣值）。

sound(x,fs,bits); 用于對聲音的回放。向量y則就代表了一個信號（也即一個復雜的%26ldquo;函數表達式%26rdquo;）也就是說可以像處理一個信號表達式一樣處理這個聲音信號。

在MATLAB的信號處理工具箱中函數FFT和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進行處理。

經函數fft求得的序列y一般是復序列，通常要求其幅值和相位。MATLAB提供求復數的幅值和相位函數：abs，angle，這些函數一般和FFT同時使用。

函數abs(x)用于計算復向量x的幅值，函數angle(x)用于計算復向量的相角，介于 和 之間，以弧度表示。

函數unwrap(p)用于展開弧度相位角p ,當相位角絕對變化超過 時，函數把它擴展至 。

用MATLAB工具箱函數fft進行頻譜分析時需注意：

（1） 函數fft返回值y的數據結構對稱性

若已知序列x=[4,3,2,6,7,8,9,0],求X(k)=DFT[x(n)]。

利用函數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.2929i

0 - 5.0000i

4.7782 - 7.7071i

5.0000

4.7782 + 7.7071i

0 + 5.0000i

-10.7782 - 6.2929i

由程序運行所得結果可見，X(k)和x(n)的維數相同，共有8個元素。X(k)的第一行元素對應頻率值為0，第五行元素對應頻率值為Nyquist頻率，即標準頻率為1.因此第一行至第五行對應的標準頻率為0～1。而第五行至第八行對應的是負頻率，其X(k)值是以Nyquist頻率為軸對稱。（注：通常表示為Nyquist頻率外擴展，標以正值。）

一般而言，對于N點的x(n)序列的FFT是N點的復數序列，其點n=N/2+1對應Nyquist頻率，作頻譜分析時僅取序列X(k)的前一半，即前N/2點即可。X(k)的后一半序列和前一半序列時對稱的。

（2） 頻率計算

若N點序列x(n)(n=0,1,…,N-1)是在采樣頻率 下獲得的。它的FFT也是N點序列，即X(k)(k=0,1,2,…,N-1)，則第k點所對應實際頻率值為f=k*f /N.

2、設計內容：

（1）下面的一段程序是語音信號在MATLAB中的最簡單表現，它實現了語音的讀入打開，以及繪出了語音信號的波形頻譜圖。

[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('原始信號波形');

subplot(222);plot(X); title('原始信號頻譜');

subplot(223);plot(magX);title('原始信號幅值');

subplot(224);plot(angX);title('原始信號相位');

程序運行可以聽到聲音，得到的圖形為：

（2）定點分析：已知一個語音信號，數據采樣頻率為100Hz，試分別繪制N＝128點DFT的幅頻圖和N＝1024點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);

subplot(221);plot(f,magy);

xlabel('頻率（Hz）');ylabel('幅值');

title('N=128(a)');grid

subplot(222);plot(f(1:N/2),magy(1:N/2));

xlabel('頻率（Hz）');ylabel('幅值');

title('N=128(b)');grid

fs=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）');ylabel('幅值');

title('N=1024(c)');grid

subplot(224);plot(f(1:N/2),magy(1:N/2));

xlabel('頻率（Hz）');ylabel('幅值');

title('N=1024(d)');grid

運行結果如圖：

上圖(a)、(b)為N=128點幅頻譜圖，(c)、(d)為N=1024點幅頻譜圖。由于采樣頻率f =100Hz，故Nyquist頻率為 50Hz。(a)、(c)是0～100Hz頻譜圖，(b)、(d)是0～50Hz頻譜圖。由(a)或(c)可見，整個頻譜圖是以Nyquist頻率為軸對稱的。因此利用fft對信號作頻譜分析，只要考察0～Nyquist頻率（采樣頻率一半）范圍的幅頻特性。比較(a)和(c)或(b)和(d)可見，幅值大小與fft選用點數N有關，但只要點數N足夠不影響研究結果。從上圖幅頻譜可見，信號中包括15Hz和40Hz的正弦分量。

（3）若信號長度T=25.6s，即抽樣后x(n)點數為T/Ts=256，所得頻率分辨率為 Hz，以此觀察數據長度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));grid

xlabel('Hz'),ylabel('|H(ejw)|')

%數據長度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));grid

xlabel('Hz'),ylabel('|H(ejw)|')

結果如圖：

（三）、濾波器設計：

1、相關原理：

設計數字濾波器的任務就是尋求一個因果穩定的線性時不變系統，并使系統函數H(z)具有指定的頻率特性。

數字濾波器從實現的網絡結構或者從單位沖激響應分類，可以分成無限長單位沖激響應（IIR）數字濾波器和有限長單位沖激響應（FIR）數字濾波器。

數字濾波器頻率響應的三個參數：

（1） 幅度平方響應：

（2） 相位響應

其中，相位響應

（3） 群時延響應

IIR數字濾波器：

IIR數字濾波器的系統函數為 的有理分數，即

IIR數字濾波器的逼近問題就是求解濾波器的系數 和 ，使得在規定的物理意義上逼近所要求的特性的問題。如果是在s平面上逼近，就得到模擬濾波器，如果是在z平面上逼近，則得到數字濾波器。

FIR數字濾波器：

設FIR的單位脈沖響應h(n)為實數，長度為N，則其z變換和頻率響應分別為

按頻域采樣定理FIR數字濾波器的傳輸函數H(z)和單位脈沖響應h(n)可由它的N個頻域采樣值H(k)唯一確定。

MATLAB中提供了幾個函數，分別用于實現IIR濾波器和FIR濾波器。

（1）卷積函數conv

卷積函數conv的調用格式為 c=conv(a,b)

該格式可以計算兩向量a和b的卷積，可以直接用于對有限長信號采用FIR濾波器的濾波。

（2）函數filter

函數filter的調用格式為 y=filter(b,a,x)

該格式采用數字濾波器對數據進行濾波，既可以用于IIR濾波器，也可以用于FIR濾波器。其中向量b和a分別表示系統函數的分子、分母多項式的系數，若a＝1，此時表示FIR濾波器，否則就是IIR濾波器。該函數是利用給出的向量b和a，對x中的數據進行濾波，結果放入向量y。

（3）函數fftfilt

該格式是利用基于FFT的重疊相加法對數據進行濾波，這種頻域濾波技術只對FIR濾波器有效。該函數是通過向量b描述的濾波器對x數據進行濾波。

關于用butter函數求系統函數分子與分母系數的幾種形式。

[b,a]=butter(N,wc,'high'):設計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的單位為 。

如果在這個函數輸入變元的最后，加一個變元%26ldquo;s%26rdquo;，表示設計的是模擬濾波器。這里不作討論。

為了設計任意的選項巴特沃斯濾波器，必須知道階數N和3dB邊緣頻率矢量wc。這可以直接利用信號處理工具箱中的buttord函數來計算。如果已知濾波器指標 ， ， 和 ，則調用格式為

[N,wc]=buttord(wp,ws,Rp,As)

對于不同類型的濾波器，參數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;ws1%26lt;ws2%26lt;wp2。

2、設計內容：

（1）濾波器示例：

在這里為了說明如何用MATLAB來實現濾波，特舉出一個簡單的函數信號濾波實例（對信號x(n)=sin( n/4)+5cos( n/2)進行濾波，信號長度為500點），從中了解濾波的實現過程。程序如下：

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('濾波前信號的頻譜');

y=filter(b,a,x);

Y=fft(y,4096);

subplot(223);plot(y);title('濾波后信號的波形');

subplot(224);plot(Y);title('濾波后信號的頻譜');

在這里，是采用了butter命令，設計出一個巴特沃斯低通濾波器，從頻譜圖中可以很明顯的看出來。下面，也就是本課題的主要內容，也都是運用到了butter函數，以便容易的得到系統函數的分子與分母系數，最終以此來實現信號的濾波。

（2）N階高通濾波器的設計（在這里，以5階為例，其中wc為其3dB邊緣頻率，以 為單位），程序設計如下：

x=wavread('ding.wav');

sound(x);

N=5;wc=0.3;

[b,a]=butter(N,wc,'high');

X=fft(x);

subplot(322);plot(X);title('濾波前信號的頻譜');

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('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('濾波前信號的頻譜');

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('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('濾波前信號的波形');

subplot(322);plot(X);title('濾波前信號的頻譜');

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濾波后信號的波形');

得到結果如圖：

（5）2N階帶阻濾波器的設計（在這里，以10階為例，其中wc為其3dB邊緣頻率，以 為單位，wc=[w1,w2],w1 wc w2），程序設計如下：

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('濾波前信號的波形');

subplot(322);plot(X);title('濾波前信號的頻譜');

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('FIR濾波后信號的頻譜');

得到結果如圖：

（6）小結：以上幾種濾波，我們都可以從信號濾波前后的波形圖以及頻譜圖上看出變化。當然，也可以用sound()函數來播放濾波后的語音，從聽覺上直接感受語音信號的變化，但由于人耳聽力的限制，有些情況下我們是很難聽出異同的。

同樣，通過函數的調用，也可以將信號的頻譜進行%26ldquo;分離觀察%26rdquo;，如顯出信號的幅值或相位。下面，通過改變系統函數的分子與分母系數比，來觀察信號濾波前后的幅值與相位。并且使結果更加明顯，使人耳得以很容易的辨聽。

x=wavread('ding.wav');

sound(x);

b=100;a=5;

y=filter(b,a,x);

X=fft(x,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));title('濾波前信號相位');

subplot(212);plot(angle(Y));title('濾波后信號相位');

得圖：

（四）、界面設計：

直接用M文件編寫GUI程序很繁瑣，而使用GUIDE設計工具可以大大提高工作效率。GUIDE相當于一個控制面板，從中可以調用各種設計工具以輔助完成界面設計任務，例如控件的創建和布局、控件屬性的編輯和菜單設計等。

使用GUIDE設計GUI程序的一般步驟如下：

1. 將所需控件從控件面板拖拽到GUIDE的設計區域；

2. 利用工具條中的工具（或相應的菜單和現場菜單），快速完成界面布局；

3. 設置控件的屬性。尤其是tag屬性，它是控件在程序內部的唯一標識；

4. 如果需要，打開菜單編輯器為界面添加菜單或現場菜單；

5. 保存設計。GUIDE默認把GUI程序保存為兩個同名文件：一個是.fig文件，用來保存窗體布局和所有控件的界面信息；一個是.m文件，該文件的初始內容是GUIDE自動產生的程序框架，其中包括了各個控件回調函數的定義。該M文件與一般的M文件沒有本質區別，但是鑒于它的特殊性，MATALAB把這類文件統稱為GUI-M文件。保存完后GUI-M文件自動在編輯調試器中打開以供編輯。

6. 為每個回調函數添加代碼以實現GUI程序的具體功能。這一步與一般函數文件的編輯調試過程相同。

設計過程及內容：

在MATLAB版面上，通過鍵入GUIDE彈出一個菜單欄進入gui制作界面（或者在File到new來進入gui）,從而開始應用界面的制作。

該界面主要實現了以下幾個功能：

①打開wav格式的音頻文件，并將該音頻信號的值讀取并賦予某一向量；

②播放音頻文件，可以選擇性的顯示該音頻信號的波形、頻譜、幅值以及相位；

③對音頻信號進行IIR與FIR的5階固定濾波處理，可以選擇性的顯示濾波前后信號的波形、頻譜、幅值以及相位，以及播放濾波后的聲音。

界面如圖所示：

通過該界面，可以方便用戶進行語音信號的處理。

界面主程序見附件。

（五）、校驗：

1、本設計圓滿的完成了對語音信號的讀取與打開，與課題的要求十分相符；

2、本設計也較好的完成了對語音信號的頻譜分析，通過fft變換，得出了語音信號的頻譜圖；

3、在濾波這一塊，課題主要是從巴特沃斯濾波器入手來設計濾波器，也從一方面基本實現了濾波；

4、初略的完成了界面的設計，但也存在相當的不足，只是很勉強的達到了打開語音文件、顯示已定濾波前后的波形等圖。

四、 結論：