MATLAB語音信號處理

數字信號處理課設，咱們使用MATLAB對語音信號進行了一系列處理，並將其全部功能集中於下圖界面中：

這個界面涉及功能衆多，其中包括語音信號的觀察分析、音色變換、AM調製解調、減抽樣、加噪去噪、相頻分析和幅頻濾波等，最重要的是對MATLAB中函數的掌握，經過不一樣函數的組合實現你想要實現的功能。

本篇不會給出整個界面的程序，下面會分塊給出每一個功能的程序，整個界面只需GUI設計界面文件、定義結構體並把對應鍵程序打進去便可。

一、語音信號的採集

1.1題目要求

使用windows下的錄音機錄製一段語音信號、音樂信號或者採用其餘軟件截取一段音樂信號（要求：時間不超過5s，文件格式爲WAV。）

① 請每位同窗都參與錄音，內容自定。

② 使用wavread語句讀取語音/音樂信號獲取抽樣率；（注意：讀取的信號是雙聲道信號，即爲雙列向量，須要分列處理）；

③ 輸出時域語音/音樂信號的波形。

④ 實現對錄音信號的聲音大小的調節。

⑤ 實現對兩種語音/音樂信號的混音音效。

⑥ 實現音樂信號的迴音音效。

1.2設計內容及方案

① 讀取音頻信號：我是經過wavread函數讀取.wav文件的方式來得到，固然首先要本身建立一個.wav音頻，我是經過電腦錄音生成.mp3而後格式工廠轉成.wav的，需保存到同一文件夾下。

② 分聲道處理：通常音樂和語音信號都是雙聲道信號，時域和頻譜圖會有兩個顏色，因此要取單列來分析，經過x1=x(:,1)語句來實現。

③ 畫時域波形圖：用plot函數來畫圖，注意橫座標爲時間t。

④ 音量大小調節：經過將音頻直接乘一個係數來實現調音量。

⑤ 混音和回聲：混音即將兩個音頻相加，要相加就得保證矩陣同樣，因此要經過截取並補零矩陣來實現；回聲是把三個信號疊加，這三個信號在不一樣位置補零音量也逐漸變小，就能夠實現回聲。

⑥ 播放聲音：本題我使用wavplay來播放聲音，會有警告，後面的題我用sound比較好。

1.3程序源碼及註釋

clear [x,fs] = wavread('beautiful.wav');%音樂信號 [y,fs1]= wavread('1.wav');%女生聲音 [z,fs2]= wavread('2.wav');%男生聲音 %輸出頻率 fs fs1 fs2 %音樂語音信號分聲道處理 x1=x(:,1); y1=y(:,1); z1=z(:,1); %畫音樂信號時域圖 n1=length(x1);%length取數列長度即元素個數 figure(1) t1=(0:(n1-1))/fs; plot(t1,x1); axis([0,5,-1,1]); xlabel('時間t'); ylabel('幅度'); title('音樂信號時域波形'); %畫語音信號時域圖 n2=length(y1); figure(2) subplot(2,1,1); t2=(0:(n2-1))/fs1; plot(t2,y1);%女生 axis([0,4,-0.5,0.5]); xlabel('時間t'); ylabel('幅度'); title('女生語音信號時域波形'); n3=length(z1); subplot(2,1,2); t3=(0:(n3-1))/fs2; plot(t3,z1);%男生 axis([0,4,-0.5,0.5]); xlabel('時間t'); ylabel('幅度'); title('男生語音信號時域波形'); %對語音信號聲音大小調節 wavplay(y,fs1); %播放原語音 y11=10*y; wavplay(y11,fs1); %加大音量播放 y22=0.5*y; wavplay(y22,fs1); %減少音量播放 %兩種語音信號的混音 [m,n]=size(y1);%size取矩陣的行列數 [m0,n0]=size(z1); a=zeros(abs(m-m0),n);%兩矩陣行數差爲零矩陣行數 if length(y1)<length(z1) y2=[y1;a]; y3=y2+z1;%兩個矩陣行數同樣才能相加，因此要補零 else y2=[z1;a]; y3=y2+y1;%y1和z1中長的那個不變，短的那個補零 end; wavplay(y3,fs1) ;%播放混音語音 %畫混音波形 figure(3) subplot(2,1,1); t4=(0:(max(n2,n3)-1))/fs1; plot(t4,y3); axis([0,4.5,-0.5,0.5]); xlabel('時間'); ylabel('幅度'); title('兩語音信號疊加後時域波形'); %音樂信號的迴音 x11=x1(1:200000);%截取部分 x11=x11';%由於輸出爲一列因此要轉置成一行
int0=zeros(1,20000);%1行2000列的零矩陣 temp1=[x11,int0,int0]; temp2=[int0,0.6*x11,int0]; temp3=[int0,int0,0.3*x11];%經過補零實現延時，同時聲音一個比一個小 hui=temp1+temp2+temp3;%三重聲音相加實現回聲 N=length(hui); wavplay(hui,fs1);%播放回音音樂 %畫回聲波形 subplot(2,1,2); t1=(0:(N-1))/fs; plot(t1,hui); axis([0,4.5,-1,1]); xlabel('時間'); ylabel('幅度'); title('回聲時域波形');

1.4運行結果

仿真結果分析：我聽到了原聲和音量放大減少的聲音，也聽到了混音和回聲的效果，變化明顯；本題我畫了音樂和兩個語音信號的時域波形以及混音回聲的時域波形，音樂信號幅度比語音信號高且連貫性高，混音以後幅度疊加，回聲以後幅度也增大，波形有很明顯的變化。

二、語音/音樂信號的頻譜和音譜的觀察

2.1題目要求

① 輸出語音/音樂信號的波形和頻譜，觀察現象；

② 使用sound語句播放語音/音樂信號，注意不一樣抽樣率下音調變化，解釋現象。

2.2設計內容及方案

本題讀取音頻信號、畫時域波形和播放原理和上題同樣，涉及的新內容有：

① 畫頻譜圖：我將橫座標設爲頻率f，縱座標須要用fft函數求傅里葉變換而後利用abs函數求幅值畫幅度譜，再用plot畫出頻譜圖。

② 調節頻率：我將頻率fs乘一個係數放大縮小並播放，感覺頻率對語音的影響。

2.3程序源碼及註釋

clear [x,fs] = wavread('beautiful.wav');%音樂信號 [y,fs1]= wavread('1.wav');%女生聲音 [z,fs2]= wavread('2.wav');%男生聲音 %輸出頻率 fs fs1 fs2 %音樂語音信號分聲道處理 x1=x(:,1); y1=y(:,1); z1=z(:,1); %畫音樂信號時域圖 n1=length(x1);%length取數列長度即元素個數 figure(1) subplot(2,1,1); t1=(0:(n1-1))/fs; plot(t1,x1); axis([0,5,-1,1]); xlabel('時間t'); ylabel('幅度'); title('音樂信號時域波形'); %畫音樂信號頻域圖 X1=fft(x1,n1); subplot(2,1,2); f1=0:fs/n1:fs*(n1-1)/n1; plot(f1,abs(X1));%也可使用fftshift函數使fft後的結果以fs/2爲中心左右互換 axis([0,44100,0,6000]); xlabel('頻率f'); ylabel('幅度'); title('音樂信號頻譜'); %畫女生語音信號時域圖 n2=length(y1); figure(2) subplot(2,2,1); t2=(0:(n2-1))/fs1; plot(t2,y1); axis([0,4,-0.5,0.5]); xlabel('時間t'); ylabel('幅度'); title('女生語音信號時域波形'); %畫女生語音信號頻域圖 Y=fft(y1,n2); subplot(2,2,2); f2=0:fs1/n2:fs1*(n2-1)/n2; plot(f2,abs(Y)); axis([0,48000,0,700]); xlabel('頻率f'); ylabel('幅度'); title('女生語音信號頻譜'); %畫男生語音信號時域圖 n3=length(z1); subplot(2,2,3); t3=(0:(n3-1))/fs2; plot(t3,z1); axis([0,4,-0.5,0.5]); xlabel('時間t'); ylabel('幅度'); title('男生語音信號時域波形'); %畫男生語音信號頻域圖 Z=fft(z1,n3); subplot(2,2,4); f3=0:fs2/n3:fs2*(n3-1)/n3; plot(f3,abs(Z)); axis([0,48000,0,700]); xlabel('頻率f'); ylabel('幅度'); title('男生語音信號頻譜'); %不一樣抽樣率語音信號 sound(y1,fs);%播放原音樂 pause(5); sound(y1,2*fs);%播放高頻率音樂 pause(2); sound(y1,0.5*fs);%播放低頻率音樂

2.4運行結果

仿真結果分析：本題中畫了頻譜圖，能夠更直觀的看到信號的特徵，包括截至頻率和音頻範圍，個人音樂和語音信號主要集中在低頻段，並且音樂信號的幅度要比語音信號的幅度高。經過放大縮小頻率的聲音變化，即頻率大時語調變高且語速加快，讓我感覺到了頻率對信號的做用。

三、語音/音樂信號的抽取（減抽樣）

3.1題目要求

① 觀察語音/音樂信號頻率的上限，選擇適當的抽取間隔對信號進行減抽樣（給出兩種抽取間隔，表明混疊與非混疊）；

② 輸出減抽樣語音/音樂信號的波形和頻譜，觀察現象，給出理論解釋；

③ 播放減抽樣語音/音樂信號，注意抽樣率改變，比較不一樣抽取間隔下的聲音，解釋現象。

3.2設計內容及方案

減抽樣：本題研究的惟一內容就是對信號以不一樣間隔抽樣，我抽取了200000長度的信號，便於以必定間隔抽樣。首先我以小間隔2抽樣，信號聲音基本和原聲差很少，沒有發生混疊；然後我又用大間隔20抽樣，信號聲音有了很明顯的變化，即發生了混疊。

3.3程序源碼及註釋

clear [x,fs]=wavread('beautiful.wav');%音樂信號 x1=x(:,1);%取單列 x2=x1(1:200000);%截取長度 N1=length(x2)%求信號長度 %畫原信號時域波形和頻譜 t1=(0:(N1-1))/fs; figure(1); subplot(2,1,1); plot(t1,x2); title('原信號時域波形'); xlabel('時間t'); ylabel('幅度'); f1=0:fs/N1:fs*(N1-1)/N1; Fx1=fft(x2,N1); %求傅里葉變換 subplot(2,1,2); plot(f1,abs(Fx1)); title('原信號的頻譜');%畫原信號頻譜 xlabel('頻率f'); ylabel('幅度'); %非混疊間隔減抽樣 d1=2;j=0; for i=1:d1:200000 j=j+1; x22(j)=x2(i);%對信號以d1爲間隔作減抽樣 end N2=length(x22);%對信號以d1爲間隔作減抽樣的長度 %畫d1間隔抽樣圖 t2=(0:(N2-1))/fs; figure(2); subplot(2,1,1); plot(t2,x22); title('以d1爲間隔減抽樣時域波形'); xlabel('時間t'); ylabel('幅度'); subplot(2,1,2); f2=0:(fs/d1)/N2:(fs/d1)*(N2-1)/N2; Fx2=fft(x22,N2);%求傅里葉變換畫頻譜 plot(f2,abs(Fx2)); title('以d1爲間隔減抽樣後的頻譜(非混疊)'); xlabel('頻率f'); ylabel('幅度'); %混疊間隔減抽樣 d2=20;j=0; for i=1:d2:200000 j=j+1; x3(j)=x2(i);%對信號以d2爲間隔作減抽樣 end N3=length(x3) ;%對信號以d2爲間隔作減抽樣長度 %畫d2間隔抽樣圖 t3=(0:(N3-1))/fs; figure(3) subplot(2,1,1); plot(t3,x3); title('以d2爲間隔減抽樣後的時域波形'); xlabel('時間t'); ylabel('幅度'); subplot(2,1,2); f3=0:(fs/d2)/N3:(fs/d2)*(N3-1)/N3; Fx3=fft(x3,N3);%求傅里葉變換畫頻譜 plot(f3,abs(Fx3)); title('以d2爲間隔減抽樣後的頻譜(混疊)'); xlabel('頻率f'); ylabel('幅度'); %播放減抽樣信號 sound(x2,fs);%原信號播放 pause(5); sound(x22,fs/d1);%以d1爲間隔減抽樣信號播放(非混疊) pause(8); sound(x3,fs*(1/d2));%以d2爲間隔減抽樣信號播放(混疊)

3.4運行結果

仿真結果分析：從信號的時域波形和頻譜圖能夠看出抽樣後時間範圍和頻率範圍都有所減少，小間隔抽樣頻譜波形變化比較小，沒有發生混疊，大間隔抽樣頻譜波形有了很大的變化，發生了混疊。抽取間隔越小，聲音越清晰，時間間隔越大，聲音越不清晰，混疊現象越明顯。未混疊時，聲音尖銳，混疊時，聲音輕，只有淡淡的音調，基本沒有起伏，不清晰。

緣由：採樣頻率fs必須大於等於2fc，採樣頻率小於2fc時，發生混疊，若採樣頻率越小，則混疊越明顯。個人圖中fc不到5000Hz,間隔2採樣頻率fs爲22000Hz大於兩倍的fc因此不混疊；間隔20採樣頻率fs爲2200Hz小於兩倍的fc因此發生混疊。

4 語音/音樂信號的AM調製

4.1題目要求

① 觀察語音/音樂信號的頻率上限，選擇適當調製頻率對信號進行調製（給出高、低兩種調製頻率）；

② 輸出調製信號的波形和頻譜，觀察現象，給出理論解釋；

③ 播放調製語音/音樂信號，注意不一樣調製頻率下的聲音，解釋現象。

4.2設計內容及方案

在這道題中我統一使用了頻率歸一化，便於從原信號中讀取截止頻率和設置載波頻率。

① 取低頻載波對信號進行AM調製：這裏取了0.1pi爲低頻調製載波頻率，與原信號相乘實現AM調製，這裏用點乘轉置矩陣實現。

② 取高頻載波對信號進行AM調製：這裏取了0.7pi爲低頻調製載波頻率，與原信號相乘實現AM調製，這裏用點乘轉置矩陣實現。

③ 播放調製後信號：分別播放低頻和高頻調製時的音樂，用sound函數播放。

4.3程序源碼及註釋

clear [x,fs]=wavread('beautiful.wav'); %讀取音樂信號 x1=x(:,1);%取單列 N=length(x1);%求信號長度 n=0:N-1;%全部元素 t=(0:(N-1))/fs;%時間 f=0:fs/N:fs*(N-1)/N;%頻率 w=2*f/fs;%歸一化 %畫原信號時域波形和頻譜 figure(1); subplot(2,1,1); plot(t,x1); title('原信號時域波形'); xlabel('時間t'); ylabel('幅度'); Fx1=fft(x1,N); %求傅里葉變換 subplot(2,1,2); plot(w,abs(Fx1)); title('原信號的頻譜');%畫原信號頻譜 xlabel('w'); ylabel('幅度'); %低頻調製 x2=cos(n*0.1*pi);%低頻時餘弦信號 x22=x1.*x2';%低頻調製信號
X2=fft(x2); X22=fft(x22); %高頻調製 x3=cos(n*0.7*pi);%高頻時餘弦信號 x33=x1.*x3';%高頻調製信號
X3=fft(x3); X33=fft(x33); %畫出低頻時載波和調製信號波形 figure(2) subplot(2,2,1); plot(t,x2); axis([0,0.001,-1,1]); title('低頻時餘弦信號時域波形'); xlabel('時間t'); ylabel('幅度'); subplot(2,2,2); plot(w,abs(X2)); title('低頻時餘弦信號頻譜'); xlabel('w'); ylabel('幅度'); subplot(2,2,3); plot(t,x22); title('低頻調製信號時域波形'); xlabel('時間t'); ylabel('幅度'); subplot(2,2,4); plot(w,abs(X22)); title('低頻調製信號頻譜'); xlabel('w'); ylabel('幅度'); %畫出高頻時載波和調製信號波形 figure(3) subplot(2,2,1); plot(t,x3); axis([0,0.001,-1,1]); title('高頻時餘弦信號時域波形'); xlabel('時間t'); ylabel('幅度'); subplot(2,2,2); plot(w,abs(X3)); title('高頻時餘弦信號頻譜'); xlabel('w'); ylabel('幅度'); subplot(2,2,3); plot(t,x33); title('高頻調製信號時域波形'); xlabel('時間t'); ylabel('幅度'); subplot(2,2,4); plot(w,abs(X33)); title('高頻調製信號頻譜'); xlabel('w'); ylabel('幅度'); %播放調製信號 sound(x1,fs);%播放原音樂 pause(5); sound(x22,fs);%低頻調製播放 pause(5); sound(x33,fs);%高頻調製播放

4.4運行結果

仿真結果分析：從圖中能夠看出原音樂信號的截止頻率爲0.2pi，這裏設置了0.1pi和0.7pi兩種頻率對信號進行AM調製，原信號的調製至關於頻譜搬移， 左移一個右移一個，調製的目的是便於信號在信道中傳輸。當調製頻率較高時，聲音響度低，幾乎只能聽見茲茲的聲音， 信號幾乎徹底失真，當調製頻率較低時，聲音很尖銳，響度較大，能聽出調子，但也有茲茲的聲音。

五、AM調製語音/音樂信號的同步解調

5.1題目要求

① 設計巴特沃斯濾波器完成同步解調，觀察濾波器頻率響應曲線；

② 窗函數法設計FIR濾波器完成同步解調，觀察濾波器頻率響應曲線（要求：分別使用矩形窗和布萊克曼窗，進行比較）；

③ 輸出解調信號的波形和頻譜，觀察現象，給出理論解釋；

④ 播放解調語音/音樂信號，比較不一樣濾波器下的聲音，解釋現象。

5.2設計內容及方案

① 對調製後的信號進行解調：將調製後的信號與調製時相同的載波相乘實現解調，這裏用點乘轉置矩陣實現。

② 用巴特沃斯濾波器對解調信號進行濾波：首先求巴特沃斯濾波器的頻率響應，其中用到了buttord求知足性能指標的濾波器階數N和3dB截止頻率wc、用butter計算模擬濾波器的傳輸函數Ha(s)、用freqz求頻響。而後用filter實現濾波。

③ 用矩形窗FIR濾波器對解調信號進行濾波：首先求矩形窗FIR濾波器的頻率響應，其中先求理想低通單位脈衝響應hd，而後加矩形窗截斷求模擬脈衝響應，再利用freqz求頻率響應。而後利用卷積conv實現濾波。

④ 用布萊克曼窗FIR濾波器對解調信號進行濾波：首先求布萊克曼窗FIR濾波器的頻率響應，其中先求理想低通單位脈衝響應hd，而後加布萊克曼窗截斷求模擬脈衝響應，再利用freqz求頻率響應。而後利用卷積conv實現濾波。

⑤ 播放調製、解調和濾波後的聲音：經過聲音變化感覺調製、解調和濾波。

5.3程序源碼及註釋

clear [y,fs]=audioread('beautiful.wav');%讀取音樂信號 y1=y(:,1);%取單列 N1=length(y1);%求信號長度 n=0:N1-1;%全部元素 t=n/fs;%時間 f=0:fs/N1:fs*(N1-1)/N1;%頻率 w=2*f/fs;%歸一化 y2=cos(n*pi*0.12);%0.12pi時餘弦信號 Y=y1.*y2';%音樂信號AM調製
Y2=Y.*y2';%解調相乘器
Y22=fft(Y2);%傅里葉變換 %畫解調信號時域波形和頻譜 figure(1); subplot(2,1,1); plot(t,Y2); xlabel('時間t'); ylabel('幅度'); title('解調信號時域波形'); subplot(2,1,2); plot(w,abs(Y22)); xlabel('w'); ylabel('幅度'); title('解調信號頻譜'); %求巴特沃斯濾波器頻率響應 wp=0.12;ws=0.3;rp=1;rs=50; %設計巴特沃斯IIR濾波器參數 [N,wc]=buttord(wp,ws,rp,rs); %求知足性能指標的濾波器階數N和3dB截止頻率wc [b,a]=butter(N,wc); %計算模擬濾波器的傳輸函數Ha(s) [Hd,w]=freqz(b,a);%求頻響 figure(2) subplot(3,1,1); plot(w/pi,abs(Hd)); axis([0,1,0,1.5]); xlabel('w/π'); ylabel('幅度'); title('巴特沃斯頻率響應曲線'); %求巴特沃斯濾波器濾波信號 Y3=filter(b,a,Y2);%濾波音樂信號 N3=length(Y3) t1=(0:(N3-1))/fs;%時間 f1=0:fs/N3:fs*(N3-1)/N3;%頻率 w1=2*f1/fs;%歸一化 Y33=fft(Y3);%傅里葉變換 subplot(3,1,2); plot(t1,Y3); xlabel('時間t'); ylabel('幅度'); title('巴特沃斯濾波信號時域波形'); subplot(3,1,3); plot(w1,abs(Y33)); xlabel('w'); ylabel('幅度'); title('巴特沃斯濾波信號頻譜'); %理想低通單位脈衝響應hd計算(窗函數法要用) N=11;n1=[0:1:(N-1)];wc1=0.2*pi; alpha=(N-1)/2; m=n1-alpha+eps;%加一個小數避免0做除數 hd=sin(wc1*m)./(pi*m); %求矩形窗FIR濾波器頻率響應 w_han1=(boxcar(N))';%矩形窗
h1=hd.*w_han1;%加窗截斷得脈衝響應 %如下爲頻率響應計算 [H,w]=freqz(h1,1,1000,'whole'); mag=abs(H);%mag爲幅度 db=20*log10((mag+eps)/max(mag)); pha=angle(H);%pha爲相位 grd=grpdelay(h1,1,w);%grd爲羣延遲 %如下爲畫圖 figure(3); subplot(3,2,1); plot(w/pi,db); axis([0,2,-200,30]); xlabel('w/π'); ylabel('幅度'); title('矩形窗頻響圖'); %求矩形窗FIR濾波器濾波信號 Y4=conv(Y2,h1);%卷積濾波 N4=length(Y4); t2=(0:(N4-1))/fs;%時間 f2=0:fs/N4:fs*(N4-1) /N4;%頻率 w2=2*f2/fs;%歸一化 Y44=fft(Y4,N4);%傅里葉變換 subplot(3,2,3); plot(t2,Y4); xlabel('時間t'); ylabel('幅度'); title('矩形窗濾波信號時域波形'); subplot(3,2,5); plot(w2,abs(Y44)); xlabel('w'); ylabel('幅度'); title('矩形窗濾波信號頻譜'); %求布萊克曼窗FIR濾波器頻率響應 w_han2=(blackman(N))';%布萊克窗
h2=hd.*w_han2;%加窗截斷得脈衝響應 %如下爲頻率響應計算 [H,w]=freqz(h2,1,1000,'whole'); mag=abs(H);%mag爲幅度 db=20*log10(mag+eps)/max(mag); pha=angle(H);%pha爲相位 grd=grpdelay(h2,1,w);%grd爲羣延遲 %如下爲畫圖 subplot(3,2,2); plot(w/pi,db); axis([0,2,-200,30]); xlabel('w/π'); ylabel('幅度'); title('布萊克曼窗頻響圖'); %求布萊克曼窗FIR濾波器濾波信號 Y5=conv(Y2,h2);%卷積濾波 N5=length(Y5); t3=(0:(N4-1))/fs;%時間 f3=0:fs/N4:fs*(N4-1)/N4;%頻率 w3=2*f3/fs;%歸一化 Y55=fft(Y5,N5);%傅里葉變換 subplot(3,2,4); plot(t3,Y5); xlabel('時間t'); ylabel('幅度'); title('布萊克曼濾波信號時域波形'); subplot(3,2,6); plot(w3,abs(Y55)); xlabel('w'); ylabel('幅度'); title('布萊克曼濾波信號頻譜圖'); %播放解調後聲音 sound(y1,fs);%原聲音 pause(5); sound(Y,fs);%AM調製後的聲音 pause(5); sound(Y2,fs);%解調後的聲音 pause(5); sound(Y3,fs);%巴特沃斯濾波器濾波後的聲音 pause(5); sound(Y4,fs);%矩形窗濾波器濾波後的聲音 pause(5); sound(Y5,fs);%布萊克曼窗濾波後的聲音

5.4運行結果

 仿真結果分析：首先獲得瞭解調後的時域波形和頻譜圖，能夠看出解調後的信號並無徹底恢復原信號，會夾雜一點調製過程當中的載波，經過濾波後信號頻譜有了很大改善。播放聲音發現：巴特沃斯濾波後聲音清晰，基本和原來的音樂差很少，可是音樂稍微低沉。巴特沃斯濾波器的特色是通頻帶的頻率響應曲線平滑。矩形窗濾波聲音較爲沉悶，也伴有雜音。緣由是矩形窗主瓣窄，旁瓣大，頻率識別精度最高，幅值識別精度最低 。布萊克曼窗濾波聲音與原信號較爲接近，聲音有些尖銳，也有輕微雜音。緣由是布萊克曼窗主瓣寬，旁瓣小，頻率識別精度最低，但幅值識別精度最高。

六、語音/音樂信號的濾波去噪

6.1題目要求

① 原始信號疊加幅度爲0.05，頻率爲3kHz,5kHz,8kHz的三餘弦混合噪聲，觀察噪聲頻譜以及加噪後語音/音樂信號的音頻和頻譜，並播放音樂，感覺噪聲對語音/音樂信號的影響；

② 給原始語音/音樂信號疊加幅度爲0.5的隨機白噪聲（可用rand語句產生），觀察噪聲頻譜以及加噪後語音/音樂信號的音頻和頻譜，並播放音樂，感覺噪聲對語音/音樂信號的影響；

③ 根據步驟①、②觀察到的頻譜，選擇合適指標設計濾波器進行濾波去噪，關注去噪後信號音譜和頻譜，並播放音樂，解釋現象。

6.2設計內容及方案

① 給信號加三餘弦混合噪聲：首先求得三餘弦混合噪聲，幅值爲0.05，分別給頻率3kHz,5kHz,8kHz，疊加獲得噪聲，而後將噪聲和原信號疊加。

② 給信號加隨機白噪聲：首先求得隨機白噪聲，幅值爲0.5，用randn語句產生噪聲，而後將噪聲和原信號疊加。

③ 濾掉噪聲：我使用了巴特沃斯濾波器來濾噪，其中用到buttord求知足性能指標的濾波器階數N和3dB截止頻率wc、用butter求s域的頻率響應的參數、用bilinear函數即利用雙線性變換實現頻率響應s域到z域的變換，而後用filter求濾波後信號。

6.3程序源碼及註釋

clear [y,fs]=audioread('beautiful.wav');%讀取音樂信號 y1=y(:,1);%取單列 N=length(y1);%求信號長度 Y=fft(y1,N);%傅里葉變換 t=(0:(N-1))/fs;%時間 f=0:fs/N:fs*(N-1)/N;%頻率 w=2*f/fs;%歸一化 %原信號時域波形和頻譜圖 figure(1); subplot(2,1,1); plot(t,y1); xlabel('時間t'); ylabel('幅度'); title('原信號時域波形'); subplot(2,1,2); plot(w,abs(Y)); xlabel('w'); ylabel('幅度'); title('原信號頻譜圖'); %加三餘弦混合噪聲 A0=0.05;%噪聲幅度 d1=A0*cos(2*pi*3000*t);%3kHz的餘弦噪音 d2=A0*cos(2*pi*5000*t);%5kHz的餘弦噪音 d3=A0*cos(2*pi*8000*t);%8kHz的餘弦噪音 x1=d1+d2+d3;%噪聲疊加 X1=fft(x1);%噪聲傅里葉變換 y2=y1+(x1)';%加三餘弦混合噪聲後信號
Y1=fft(y2);%加噪信號傅里葉變換 figure(2); subplot(3,2,1); plot(w,abs(X1)); xlabel('w'); ylabel('幅度'); title('三餘弦噪聲頻譜'); subplot(3,2,3); plot(t,y2); xlabel('時間t'); ylabel('幅度'); title('加三餘弦噪聲信號時域波形'); subplot(3,2,5); plot(w,abs(Y1)); xlabel('w'); ylabel('幅度'); title('加三餘弦噪聲信號頻譜圖'); %加隨機白噪聲 x2=0.5*randn(N,1);%噪聲 X2=fft(x2);%噪聲傅里葉變換 y3=y1+x2;%加隨機噪聲後信號 Y2=fft(y3);%加噪信號傅里葉變換 subplot(3,2,2); plot(w,abs(X2)); xlabel('w'); ylabel('幅度'); title('隨機白噪聲頻譜圖'); subplot(3,2,4); plot(t,y3); xlabel('時間t'); ylabel('幅度'); title('加隨機噪聲信號時域波形'); subplot(3,2,6); plot(w,abs(Y2)); xlabel('w'); ylabel('幅度'); title('加隨機噪聲信號頻譜圖'); %巴特沃斯濾波器去噪 wp=0.1;ws=0.16;Rp=1;Rs=15; [n1,Wn1]=buttord(wp,ws,Rp,Rs,'s');%求低通濾波器的階數和截止頻率 [b1,a1]=butter(n1,Wn1,'s');%求s域的頻率響應的參數 [b2,a2]=bilinear(b1,a1,1);%利用雙線性變換實現頻率響應s域到z域的變換 z1=filter(b2,a2,y2);%求濾三餘弦混合噪聲後的信號 z2=filter(b2,a2,y3);%求濾隨機噪聲後的信號 %畫濾噪後的圖 m1=fft(z1); m2=fft(z2); figure(3); subplot(2,2,1); plot(t,z1); xlabel('時間t'); ylabel('幅度'); title('濾三餘弦噪聲後信號時域波形'); subplot(2,2,3); plot(w,abs(m1)); xlabel('w'); ylabel('幅度'); title('濾三餘弦噪聲後信號頻譜'); subplot(2,2,2); plot(t,z2); xlabel('時間t'); ylabel('幅度'); title('濾隨機噪聲後信號時域波形'); subplot(2,2,4); plot(w,abs(m2)); xlabel('w'); ylabel('幅度'); title('濾隨機噪聲後信號頻譜'); %播放 sound(y1,fs);%原音樂 pause(5); sound(y2,fs);%加三餘弦噪聲 pause(5); sound(y3,fs);%加隨機噪聲 pause(5); sound(z1,fs);%濾三餘弦噪聲 pause(5); sound(z2,fs);%濾隨機噪聲

6.4運行結果

仿真結果分析：從圖中能夠看到三餘弦混合噪聲的頻譜圖以及信號加三餘弦混合噪聲後的時域波形和頻譜，經過播放感覺到了噪聲對信號的影響；濾波以後對噪聲的改善很大，噪聲變小，原聲音更加清晰，只是巴特沃斯濾波會把一部分原信號頻率濾掉，聲音會有點低沉。

七、語音/音樂信號的幅頻濾波及相頻分析

7.1題目要求

① 設計低通濾波器（可自行選擇不一樣的截止頻率），濾除原始語音/音樂信號的高頻信息，觀察濾波先後的幅度頻譜，並比較濾波先後的音樂效果，感覺高頻信息對語音/音樂信號的影響；

② 設計高通濾波器（可自行選擇不一樣的截止頻率），濾除原始語音/音樂信號的低頻信息，觀察濾波先後的幅度頻譜，並比較濾波先後的音樂效果，感覺低頻信息對語音/音樂信號的影響；

③ 選取兩段不一樣的語音/音樂信號，分別將其幅度譜與相位譜交叉組合構成新的語音/音樂信號，播放比較組合後的音樂與原始音樂，感覺相頻信息對語音/音樂信號的影響。

7.2設計內容及方案

① 低通濾波器設計：我這裏用了巴特沃斯低通濾波器，其中用buttord求低通濾波器的階數和截止頻率，用butter求s域的頻率響應的參數，用bilinear即雙線性變換法實現頻率響應s域到z域的變換，用filter實現濾波。

② 高通濾波器設計：我這裏用了巴特沃斯低通濾波器轉高通，其中用buttord求低通濾波器的階數和截止頻率，用buttap建立巴特沃斯低通濾波器原型，用zp2tf將模擬低通變高通，用bilinear即雙線性變換法實現頻率響應s域到z域的變換，用filter實現濾波。

③ 幅度譜和相位譜交叉：用函數abs和angle分別取信號的幅度和相位，而後將其交叉建立新的信號。

7.3程序源碼及註釋

clear [y,fs]=audioread('忽然好想你.wav');%讀取音樂信號 y1=y(:,1);%取單列 N=length(y1);%求信號長度 Y=fft(y1,N);%傅里葉變換 t=(0:(N-1))/fs;%時間 f=0:fs/N:fs*(N-1)/N;%頻率 w=2*f/fs;%歸一化 %IIR低通濾波器 wp=0.1;ws=0.2;Rp=1;Rs=15; [n1,Wn1]=buttord(wp,ws,Rp,Rs,'s');%求低通濾波器的階數和截止頻率 [b1,a1]=butter(n1,Wn1,'s');%求s域的頻率響應的參數 [b2,a2]=bilinear(b1,a1,1);%利用雙線性變換實現頻率響應s域到z域的變換 z1=filter(b2,a2,y1);%濾波後信號 m1=fft(z1); figure(1); subplot(2,1,1); plot(w,abs(Y)); xlabel('w'); ylabel('幅度'); title('原信號頻譜'); subplot(2,1,2); plot(w,abs(m1)); xlabel('w'); ylabel('幅度'); title('低通濾波信號頻譜'); %高通濾波器 wp=0.1;ws=0.2;rp=1;rs=15; [n,Wn]=buttord(wp,ws,rp,rs,'s');%設計模擬濾波器 [z,p,k]=buttap(n);%建立巴特沃斯低通濾波器原型 [b,a]=zp2tf(z,p,k);%由零極點轉換爲傳遞函數 [b11,a11]=lp2hp(b,a,Wn);%模擬低通變高通 [b22,a22]=bilinear(b11,a11,0.5);%利用雙線性變換實現頻率響應S域到Z域的變換 z2=filter(b22,a22,y1);%濾波後信號 m2=fft(z2); figure(2); subplot(2,1,1); plot(w,abs(Y)); xlabel('w'); ylabel('幅度'); title('原信號頻譜'); subplot(2,1,2); plot(w,abs(m2)); xlabel('w'); ylabel('幅度'); title('高通濾波信號頻譜'); %幅度譜相位譜交叉 [x,fs0]=audioread('beautiful.wav'); x0=x(:,1); x1=x0(1:200000);%信號2截取長度200000 y2=y1(1:200000);%信號1截取長度200000 N=length(x1); X=fft(x1,N); Y=fft(y2,N); t0=(0:(N-1))/fs0; f0=0:fs0/N:fs0*(N-1)/N; w0=2*f0/fs0;%數字角頻率 %畫兩信號的幅度譜和相位譜 figure(3); subplot(2,2,1); plot(w0,abs(Y)); xlabel('w'); ylabel('幅度'); title('信號1的幅度譜'); subplot(2,2,2); plot(w0,abs(X)); xlabel('w'); ylabel('幅度'); title('信號2的幅度譜'); subplot(2,2,3); plot(w0,angle(Y)); xlabel('w'); ylabel('相位'); title('信號1的相位譜'); subplot(2,2,4); plot(w0,angle(X)); xlabel('w'); ylabel('相位'); title('信號2的相位譜'); F1=abs(Y).*exp(1i*angle(X));%音樂信號1的幅度譜和音樂信號2的相位譜交叉 F2=abs(X).*exp(1i*angle(Y));%音樂信號2的幅度譜和音樂信號1的相位譜交叉y3=real(ifft(F1,N));%對交叉獲得的頻域信號作傅里葉逆變換 y4=real(ifft(F2,N));%對交叉獲得的頻域信號作傅里葉逆變換 %繪製交叉獲得的信號的波形和頻譜 figure(4) subplot(2,2,1); plot(t0,y3); xlabel('時間'); ylabel('幅度'); title('1幅度譜和2相位譜交叉信號時域波形'); subplot(2,2,2); plot(w0,abs(F1)); xlabel('w'); ylabel('幅度'); title('1幅度譜和2相位譜交叉信號頻譜'); subplot(2,2,3); plot(t0,y4); xlabel('時間'); ylabel('幅度') title('2幅度譜和1相位譜交叉信號時域波形'); subplot(2,2,4); plot(w0,abs(F2)); xlabel('w'); ylabel('幅度'); title('2幅度譜和1相位譜交叉信號頻譜'); %播放 sound(y1,fs);%原音樂 pause(5); sound(z1,fs);%低通濾波聲音 pause(5); sound(z2,fs);%高通濾波聲音 pause(5); sound(y3,fs);%幅度譜1與相位譜2交叉 pause(5); sound(y4,fs);%相位譜1與幅度譜2交叉

7.4運行結果

仿真結果分析：經過觀察原信號的歸一化頻譜，肯定巴特沃斯高通濾波器的參數wp,ws的值並實現濾波，從低通濾波和高通濾波的頻譜圖中能夠看出：低通濾波器濾掉了信號的高頻部分，聲音變得低沉；高通濾波器濾掉了信號的低頻部分，聲音變得尖銳。幅度譜與相位譜交叉時，經過聽交叉後的語音讓我感覺到了相頻特性對一個信號的影響，音樂幅度譜沒變相位譜變還會有原聲，只是總體節奏改變。

八、語音信號變聲處理系統

8.1題目要求

① 電視臺常常針對某些事件的知情者進行採訪，爲了保護知情者，常常改變說話人的聲音，請利用所學的知識，將其實現；

② 要求處理後的語音信號基本不影響正常收聽與理解。

8.2設計內容及方案

變聲主要經過調用voice函數來實現，改變voice函數的參數，小於1時偏女聲，大於1時偏男聲。

8.3程序源碼及註釋

 

clear [y1,fs]=audioread('beautiful.wav');%讀取音樂信號 y1=y1(:,1);%取單列 N1=length(y1);%求信號長度 n=0:N1-1;%全部元素 t=n/fs;%時間 f=0:fs/N1:fs*(N1-1)/N1;%頻率 Y1=fft(y1,N1);%傅里葉變換 %畫原信號圖 figure(1) subplot(2,2,1); plot(t,y1); xlabel('時間t'); ylabel('幅度'); title('原信號時域波形圖'); subplot(2,2,2); plot(f,abs(Y1)); xlabel('頻率f'); ylabel('幅度'); title('原信號頻譜圖'); y2=voice(y1,0.5);  %變聲 N2=length(y2); k=0:1:N2-1; t=k/fs; f=k*fs/N2; Y2=fft(y2,N2); %畫變音後信號圖 subplot(2,2,3); plot(t,y2); xlabel('時間t'); ylabel('幅度'); title('變聲信號時域波形圖'); subplot(2,2,4); plot(f,abs(Y2)); xlabel('頻率f'); ylabel('幅度'); title('變聲信號頻譜圖'); sound(y2,fs);

 

voice函數：

function Y=voice(y1,f)      %更改採樣率使基頻改變 f>1下降;f<1升高? f=round(f*1000); d=resample(y1,f,1000);      %時長整合使語音文件恢復原來時長 W=400; Wov=W/2; Kmax=W*2; Wsim=Wov; xdecim=8; kdecim=2; X=d';
F=f/1000; Ss=W-Wov; xpts=size(X,2); ypts=round(xpts/F); Y=zeros(1,ypts); xfwin=(1:Wov)/(Wov+1); ovix=(1-Wov):0; newix=1:(W-Wov); simix=(1:xdecim:Wsim)-Wsim; padX=[zeros(1,Wsim),X,zeros(1,Kmax+W-Wov)]; Y(1:Wsim)=X(1:Wsim); lastxpos=0; km=0; for ypos=Wsim:Ss:(ypts-W) xpos=round(F*ypos); kmpred=km+(xpos-lastxpos); lastxpos=xpos; if(kmpred<=Kmax) km=kmpred; else ysim=Y(ypos+simix); rxy=zeros(1,Kmax+1); rxx=zeros(1,Kmax+1); Kmin=0; for k=Kmin:kdecim:Kmax xsim=padX(Wsim+xpos+k+simix); rxx(k+1)=norm(xsim); rxy(k+1)=(ysim*xsim');
end Rxy=(rxx~=0).*rxy./(rxx+(rxx==0)); km=min(find(Rxy==max(Rxy))-1); end xabs=xpos+km; Y(ypos+ovix)=((1-xfwin).*Y(ypos+ovix))+(xfwin.*padX(Wsim+xabs+ovix)); Y(ypos+newix)=padX(Wsim+xabs+newix); end end

8.4運行結果

仿真結果分析：該程序能夠很好的實現對聲音信號的變聲處理，變聲後語速基本沒有改變，能夠聽清講述者想傳達的內容，只是音色有了變化。

小結：一次收穫不少的課設，但願對你有用。