數字信號處理課設，我們使用MATLAB對語音信號進行了一系列處理，並將其所有功能集中於下圖界面中：
這個界面涉及功能眾多，其中包括語音信號的觀察分析、音色變換、AM調制解調、減抽樣、加噪去噪、相頻分析和幅頻濾波等，最重要的是對MATLAB中函數的掌握，通過不同函數的組合實現你想要實現的功能。
本篇不會給出整個界面的程序，下面會分塊給出每個功能的程序，整個界面只需GUI設計界面文件、定義結構體並把對應鍵程序打進去即可。
1、語音信號的采集
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、語音/音樂信號的頻譜和音譜的觀察
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、語音/音樂信號的抽取（減抽樣）
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調制，原信號的調制相當於頻譜搬移， 左移一個右移一個，調制的目的是便於信號在信道中傳輸。當調制頻率較高時，聲音響度低，幾乎只能聽見茲茲的聲音， 信號幾乎完全失真，當調制頻率較低時，聲音很尖銳，響度較大，能聽出調子，但也有茲茲的聲音。
5、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、語音/音樂信號的濾波去噪
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、語音/音樂信號的幅頻濾波及相頻分析
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、語音信號變聲處理系統
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運行結果
仿真結果分析：該程序可以很好的實現對聲音信號的變聲處理，變聲后語速基本沒有改變，可以聽清講述者想傳達的內容，只是音色有了變化。
小結：一次收獲很多的課設，希望對你有用。