1、MATLAB課程設計報告課題：語音信號采集與處理目錄一、實踐目的3二、實踐原理：3三、課題要求：3四、MATLAB仿真41、頻譜分析：42、調制與解調：53、信號變化：8快放：8慢放：8倒放：8回聲：9男女變聲：94、信號加噪115、用窗函數法設計FIR濾波器12FIR低通濾波器：13FIR高通濾波器：14FIR帶通濾波：15一、實踐目的本次課程設計的課題為基于MATLAB的語音信號采集與處理，學會運用MATLAB的信號處理功能，采集語音信號，并對語音信號進行濾波及變換處理，觀察其時域和頻域特性，加深對信號處理理論的理解，并為今后熟練使用MATLAB進行系統的分析仿真和設計奠定基礎。此次實習課

2、程主要是為了進一步熟悉對matlab軟件的使用，以及學會利用matlab對聲音信號這種實際問題進行處理，將理論應用于實際，加深對它的理解。二、實踐原理：利用MATLAB對語音信號進行分析和處理，采集語音信號后，利用MATLAB軟件平臺進行頻譜分析；并對所采集的語音信號加入干擾噪聲，對加入噪聲的信號進行頻譜分析，設計合適的濾波器濾除噪聲，恢復原信號。語音信號的“ 短時譜”對于非平穩信號, 它是非周期的, 頻譜隨時間連續變化, 因此由傅里葉變換得到的頻譜無法獲知其在各個時刻的頻譜特性。如果利用加窗的方法從語音流中取出其中一個短斷, 再進行傅里葉變換, 就可以得到該語音的短時譜。三、課題要求：利用w

3、indows 自帶的錄音機或者其它錄音軟件，錄制幾段語音信號（要有幾種不同的聲音，要有男聲、女聲）。對錄制的語音信號進行頻譜分析，確定該段語音的主要頻率范圍，由此頻率范圍判斷該段語音信號的特點（低沉or 尖銳）。利用采樣定理，對該段語音信號進行采樣，觀察不同采樣頻率（過采樣、欠采樣、臨界采樣）對信號的影響。對采集到的語音信號進行調制與解調，觀測調制與解調前后信號的變化。實現語音信號的快放、慢放、倒放、回聲、男女變聲。對語音信號加噪，然后進行濾波，分析不同的濾波方式對信號的影響。利用MATLAB GUI 制作語音信號采集與分析演示系統。四、MATLAB仿真1、頻譜分析：用WINDOWS下的錄音機

4、，用單聲道錄制一段音樂或聲音，時間在5S內。然后MATLAB軟件平臺下，利用函數WAVREAD對語音信號進行采樣，記住采樣頻率和采樣點數。對語音信號進行快速傅立葉變換，在一個窗口同時畫出信號的時域波形圖和頻譜圖，分析語音信號的頻譜特點程序：fs =22050;Nbits =16;x,fs,Nbits =wavread('D:matlab22hexian.wav') ; %讀聲音文件n=length(x);t=0:1/fs:(length(x)-1)/fs; %求出語音信號的長度y1=fft(x,n) ; %傅里葉變換y2=fftshift(y1); %對頻譜圖進行平移f=0:f

5、s/n:fs*(n-1)/n; %得出頻點 subplot(2,1,1);plot(t/2,x) %做原始語音信號的時域圖形title('原始信號時域波形圖');subplot(2,1,2);plot(f,abs(y2);title('原始信號頻譜圖')仿真波形：門鈴：和弦：男女聲： 2、調制與解調：首先畫出語音信號的時域波形，然后對語音信號進行頻譜分析。在Matlab 中可以利用函數fft 對信號行快速傅里葉變換,得到信號的頻譜特性，從而加深對頻譜特性的理解。程序：clear; dt=1/44100; fs=44100; f1,fs,nbi

6、ts=wavread('D:1huan.wav'); figure(1); subplot(1,1,1); N=length(f1); t=0:1/fs:(N-1)/fs;plot(t,f1); title('信息信號的時域波形');fy1=fft(f1);w1=0:fs/(N-1):fs;figure(2);subplot(1,1,1); plot(w1,abs(fy1);title('信息信號的頻譜'); f2=cos(22000*pi*t);figure(3);subplot(1,1,1); fy2 = fft(f2);N2=length(

7、f2);w2=fs/N*0:N-1;plot(w2,abs(abs(fy2); title('載波信號的頻譜');f1=f1(:,1);f3=f1'.*f2;figure(4);subplot(1,1,1); fy3 = fft(f3);plot(w1,abs(abs(fy3); title('已調信號的頻譜');sound(f3,fs,nbits);f4=f3.*f2;figure(5);subplot(1,1,1);fy4=fft(f4);plot(w1,abs(abs(fy4);title('解調信號的頻譜');sound(f4,f

8、s,nbits);fp1=0; fs1=5000; As1=100;wp1=2*pi*fp1/fs; ws1=2*pi*fs1/fs; BF1=ws1-wp1;wc1=(wp1+ws1)/2;M1=ceil(As1-7.95)/(2.286*BF1)+1;N1=M1+1;beta1=0.1102*(As1-8.7); Window=(kaiser(N1,beta1); b1=fir1(M1,wc1/pi,Window);figure(6);subplot(1,1,1);freqz(b1,1,512); title('FIR低通濾波器的頻率響應'); f4_low = filte

9、r(b1,1, f4);plot(t,f4_low);title('濾波后的解調信號時域波形');sound(f4_low,fs,nbits); f5=fft(f4_low);figure(7);subplot(1,1,1);plot(w1,abs(f5);title('濾波后的解調信號頻譜'); 仿真波形： 3、信號變化：快放：x,fs,nbits=wavread('D:1huan.wav')w=2；M=w*fs;wavplay(x,M);慢放：x,fs,nbits=wavread('D:1huan.wav')w=0.8M=w*

10、fs;wavplay(x,M);倒放：x,fs,Nbits=wavread('D:1huan.wav');y0=flipud(x);sound(y0,fs);回聲：程序：x,fs,bits=wavread('D:3yang.wav',1 40000);%讀取語音信號n1=0:2000;b=x(:,1); %產生單聲道信號N=3;yy2=filter(1,1,zeros(1,80000/(N+1),0.7,b',zeros(1,40000);figure(3)subplot(2,1,1);plot(yy2); %三次回聲濾波器時域波形title('

11、三次回聲濾波器時域波形');YY2=fft(yy2); %對三次回聲信號做FFT變換subplot(2,1,2);plot(n1(1:1000),YY2(1:1000); %三次回聲濾波器頻譜圖title('三次回聲濾波器頻譜圖');figure(4)subplot(2,1,1);plot(abs(YY2); %經傅里葉變換之后的信號的幅值title('幅值');subplot(2,1,2);plot(angle(YY2); %經傅里葉變換之后的信號的相位title('相位');sound(2*yy2,fs,bits);%經三次回聲濾波器

12、后的語音信號，乘以2是為了加強信號仿真波形： 男女變聲：Voice調用函數：function Y=voice(x,f) %更改采樣率使基頻改變 f>1降低;f<1升高f=round(f*1000);d=resample(x,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 = (

13、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)

14、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 =

15、 xpos+km; Y(ypos+ovix) = (1-xfwin).*Y(ypos+ovix) + (xfwin.*padX(Wsim+xabs+ovix); Y(ypos+newix) = padX(Wsim+xabs+newix); endend變聲程序：y,fs,nbits=wavread('3yang.wav'); %讀取聲音文件x=y(:,1); %讀入的y矩陣有兩列,取第1列sound(voice(x,1.4),fs,nbits); %調整voice()第2個參數轉換音調,>1降調,<1升調4、信號加噪程序：fs=22050;  %語音信號采

16、樣頻率為22050x1=wavread('3yang.wav'); %讀取語音信號數據賦值給x1f=fs*(0:511)/1024;  %將0到511，步長為1的序列的值與fs相乘并除以1024的值賦給ft=0:1/fs:(length(x1)-1)/fs; %將0到x1的長度減1后的值除以fs的值，且步長為1/fs的值，的序列的值，賦予tAu=0.05; %噪聲幅值 d=Au*cos(2*pi*2000*t)'  %干擾信號構建命令函數，構建了一個余弦函數x2=x1(:,1)+d; %取原始語音信號的單聲部信號然后和噪聲信號相加wavwr

17、ite(x2,22050,'jiazaoyang');%生成wav文件sound(x2,22050);  %播放語音信號y1=fft(x1,1024);  %對信號做1024點的FFT變換y2=fft(x2,1024); figure(1);  %創建圖形窗subplot(2,1,1); plot(t,x1);  %做原始語音信號的時域波形title('加噪前的信號'); xlabel('time n');  %x軸的名字是timeylabel('f

18、uzhi n');  %y軸的名字是fuzhisubplot(2,1,2); %創建兩行一列繪圖區間的第1個繪圖區間plot(t,x2) title('加噪后的信號'); xlabel('time n'); ylabel('fuzhi n'); figure(2) subplot(2,1,1); plot(f,abs(y1(1:512); title('原始語音信號頻譜'); xlabel('Hz

19、'); ylabel('fuzhi'); subplot(2,1,2); plot(f,abs(y2(1:512); title('加噪后的信號頻譜'); xlabel('Hz'); ylabel('fuzhi');仿真波形： 5、用窗函數法設計FIR濾波器根據過渡帶寬及阻帶衰減要求，選擇窗函數的類型并估計窗口長度N（或階數M=N-1），窗函數類型可根據最小阻帶衰減As獨立選擇，因為窗口長度N對最小阻帶衰減As沒有影響，在確定窗函數類型以后，可根據過渡帶寬小于給定

20、指標確定所擬用的窗函數的窗口長度N，設待求濾波器的過渡帶寬為w，它與窗口長度N近似成反比，窗函數類型確定后，其計算公式也確定了，不過這些公式是近似的，得出的窗口長度還要在計算中逐步修正，原則是在保證阻帶衰減滿足要求的情況下，盡量選擇較小的N，在N和窗函數類型確定后，即可調用MATLAB中的窗函數求出窗函數wd（n）。根據待求濾波器的理想頻率響應求出理想單位脈沖響應hd(n)，如果給出待求濾波器頻率應為Hd，則理想的單位脈沖響應可以用下面的傅里葉反變換式求出： 在一般情況下，hd(n)是不能用封閉公式表示的，需要采用數值方法表示；從w=0到w=2采樣N點，采用離散傅里葉反變換（IDFT）即可求出

21、。用窗函數wd(n)將hd(n)截斷，并進行加權處理，得到如果要求線性相位特性， 則h(n)還必須滿足：根據上式中的正、 負號和長度N的奇偶性又將線性相位FIR濾波器分成四類。 要根據所設計的濾波特性正確選擇其中一類。 例如， 要設計線性相位低通特性可選擇h(n)=h(N-1-n)一類，而不能選h(n)=-h(N-1-n)一類。 驗算技術指標是否滿足要求，為了計算數字濾波器在頻域中的特性，可調用freqz子程序，如果不滿足要求，可根據具體情況，調整窗函數類型或長度，直到滿足要求為止。FIR低通濾波器：程序：x1,Fs,bits=wavread('jiazapyang');

22、60;derta_Fs = Fs/length(x1);%設置頻譜的間隔，分辨率 ，這里保證了x軸的點數必須和y軸點數一致fs=Fs;  fp1=1000; fs1=1200; As1=100;wp1=2*pi*fp1/fs; %ws1=2*pi*fs1/fs; %BF1=ws1-wp1;wc1=(wp1+ws1)/2;M1=ceil(As1-7.95)/(2.286*BF1)+1;%按凱澤窗計算濾波器階數N1=M1+1;beta1=0.1102*(As1-8.7); Window=(kaise

23、r(N1,beta1);  %求凱澤窗窗函數b1=fir1(M1,wc1/pi,Window); %wc1/pi為歸一化,窗函數法設計函數figure(2); freqz(b1,1,512); %H,w=freqz(B,A,N),(1)中B和A分別為離散系統的系統函數分子、分母多項式的系數向量，返回量H則包含了離散系統頻響在 0pi范圍內N個頻率等分點的值（其中N為正整數），w則包含了范圍內N個頻率等分點。調用默認的N時，其值是512。title('FIR低通濾波器的頻率響應'); x1_low = filt

24、er(b1,1, x1);%對信號進行低通濾波 ,Y = filter(B,A,X) ，輸入X為濾波前序列，Y為濾波結果序列，B/A 提供濾波器系數，B為分子， A為分母 sound(x1_low,Fs,bits); figure(3);subplot(2,1,1);plot(x1_low);title('信號經過FIR低通濾波器(時域)'); subplot(2,1,2);plot(-Fs/2:derta_Fs: Fs/2-derta_Fs,abs(fftshift(

25、fft(x1_low); title('信號經過FIR低通濾波器（頻域）'); 仿真波形： FIR高通濾波器：程序：x1,Fs,bits=wavread('jiazaoyang.wav'); derta_Fs = Fs/length(x1);%設置頻譜的間隔，分辨率 ，這里保證了x軸的點數必須和y軸點數一致fs=Fs; As2=100;fp2=4000; fs2=2000;wp2=2*pi*fp2/fs; ws2=2*pi*fs2/fs; BF2=wp2-ws2;wc2=(

26、wp2+ws2)/2;M2=ceil(As2-7.95)/(2.286*BF2)+1;%按凱澤窗計算濾波器階數N2=M2+1;beta2=0.1102*(As2-8.7); Window=(kaiser(N2,beta2); %求凱澤窗窗函數b2=fir1(M2,wc2/pi,'high',Window);figure(4); freqz(b2,1,512);%數字濾波器頻率響應 title('FIR高通濾波器的頻率響應'); x1_high = filter(b2,1,x1);%對信號進行高通濾波sound(x1_high,Fs,bits);figure(5);subplot(211);plot(x1_high);title('信號經過FIR高通濾波器(時域)');subplot(212);plot(-Fs/2:derta_Fs: Fs/2-derta_Fs,abs(f