數位訊號處理課設，我們使用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');%男生聲音%輸出頻率fsfs1fs2%音樂語音信號分聲道處理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');%男生聲音%輸出頻率fsfs1fs2%音樂語音信號分聲道處理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為間隔做減抽樣endN2=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為間隔做減抽樣endN3=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截取長度200000y2=y1(1:200000);%信號1截取長度200000N=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;elseysim=Y(ypos+simix);rxy=zeros(1,Kmax+1);rxx=zeros(1,Kmax+1);Kmin=0;for k=Kmin:kdecim:Kmaxxsim=padX(Wsim+xpos+k+simix);rxx(k+1)=norm(xsim);rxy(k+1)=(ysim*xsim');endRxy=(rxx~=0).*rxy./(rxx+(rxx==0));km=min(find(Rxy==max(Rxy))-1);endxabs=xpos+km;Y(ypos+ovix)=((1-xfwin).*Y(ypos+ovix))+(xfwin.*padX(Wsim+xabs+ovix));Y(ypos+newix)=padX(Wsim+xabs+newix);endend
8.4運行結果
仿真結果分析：該程序可以很好的實現對聲音信號的變聲處理，變聲後語速基本沒有改變，可以聽清講述者想傳達的內容，只是音色有了變化。