基于matlab的語音信號濾波處理——數(shù)字信號處理課程設(shè)計

1、通信與信息工程學(xué)院數(shù)字信號處理課程設(shè)計(綜合實驗)班 級：姓 名：學(xué) 號：指導(dǎo)教師：設(shè)計時間：成 績：評 語：通信與信息工程學(xué)院二一三年設(shè)計一 正余弦信號的譜分析【一】 設(shè)計目的1. 用DFT實現(xiàn)對正余弦信號的譜分析；2. 觀察DFT長度和窗函數(shù)長度對頻譜的影響；3. 對DFT進(jìn)行譜分析中的誤差現(xiàn)象獲得感性認(rèn)識?！径?設(shè)計原理一、譜分析原理數(shù)字信號處理方法的一個重要用途是在離散時間域中確定一個連續(xù)時間信號的頻譜，通常稱為頻譜分析，更具體地說，它也包括確定能量譜和功率譜。數(shù)字頻譜分析可以應(yīng)用在很廣闊領(lǐng)域，頻譜分析方法是基于以下的觀測：如果連續(xù)時間信號是頻帶有限的，那么對其離散時間等效信號的DF

2、T進(jìn)行譜分析。然而，在大多數(shù)情況下，是在范圍內(nèi)定義的，因此也就定義在的無線范圍內(nèi)，要估計一個無限長信號的頻譜是不可能的。實用的方法是：先讓模擬連續(xù)信號通過一個抗混疊的模擬濾波器，然后把它采樣成一個離散序列。假定反混疊濾波器的設(shè)計是正確的，則混疊效應(yīng)可以忽略，又假設(shè)A/D變換器的字長足夠長，則A/D變換中的量化噪聲也可忽略。假定表征正余弦信號的基本參數(shù)，如振幅、頻率和相位不隨時間改變，則此信號的傅立葉變換可以用計算它的DTFT得到 （1.1）實際上無限長序列首先乘以一個長度為M的窗函數(shù)，使它變成一個長為M的有限長序列，對求出的DTFT 應(yīng)該可以作為原連續(xù)模擬信號的頻譜估計，然后求出在區(qū)間等分為N

3、點(diǎn)的離散傅立葉變換DFT。為保證足夠的分辨率，DFT的長度N選的比窗長度M大，其方法是在截斷了的序列后面補(bǔ)上NM個零。計算采用FFT算法?！救?設(shè)計結(jié)果1.（1）程序：>> N=input('輸入譜分析的長度')n=0:N-1;T=1/64;x1=cos(20*pi*n*T);subplot(2,2,1);stem(n,x1); xlabel('n');ylabel('x1(n)'); %繪制序列的波形title('余弦序列');X1=abs(fft(x1,N); %求余弦序列的N點(diǎn)FFTsubplot(2,2,2)

4、k=0:N-1;stem(k,X1) %繪制序列的幅頻特性曲線xlabel('k');ylabel('X(k)'); string=num2str(N),'點(diǎn)FFT幅頻曲線'title(string);x2=cos(22*pi*n*T);subplot(2,2,3);stem(n,x2);xlabel('n');ylabel('x1(n)'); %繪制序列的波形title('余弦序列');X2=abs(fft(x2,N); %求余弦序列的N點(diǎn)FFTsubplot(2,2,4)k=0:N-1;stem

5、(k,X2) %繪制序列的幅頻特性曲線xlabel('k');ylabel('X(k)'); string=num2str(N),'點(diǎn)FFT幅頻曲線'title(string);（2）調(diào)試結(jié)果：2.（1）程序：f1=0.22;f2=0.34;n=0:16;x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n);X1=abs(fft(x,16); subplot(2,2,1),stem(n,x), xlabel('n');ylabel('x1(n)'); k=0:15;stem(k,X1) %繪制序列

6、的幅頻特性曲線xlabel('k');ylabel('X(k)'); string=num2str(16),'點(diǎn)FFT幅頻曲線'title(string);X2=abs(fft(x,32); subplot(2,2,2),stem(n,x), xlabel('n');ylabel('x1(n)'); k=0:31;stem(k,X2) %繪制序列的幅頻特性曲線xlabel('k');ylabel('X(k)'); string=num2str(32),'點(diǎn)FFT幅頻曲線

7、9;title(string);X3=abs(fft(x,64); subplot(2,2,3),stem(n,x), xlabel('n');ylabel('x1(n)'); k=0:63;stem(k,X3) %繪制序列的幅頻特性曲線xlabel('k');ylabel('X(k)'); string=num2str(64),'點(diǎn)FFT幅頻曲線'title(string);X4=abs(fft(x,128); subplot(2,2,4),stem(n,x), xlabel('n');ylabe

8、l('x1(n)'); k=0:127;stem(k,X4) %繪制序列的幅頻特性曲線xlabel('k');ylabel('X(k)'); string=num2str(128),'點(diǎn)FFT幅頻曲線'title(string);（2）調(diào)試結(jié)果：3.【四】 設(shè)計心得設(shè)計二 數(shù)字濾波器的設(shè)計及實現(xiàn)【一】設(shè)計目的1. 熟悉IIR數(shù)字濾波器和FIR數(shù)字濾波器的設(shè)計原理和方法；2. 學(xué)會調(diào)用MATLAB信號處理工具箱中的濾波器設(shè)計函數(shù)設(shè)計各種IIR和FIR數(shù)字濾波器，學(xué)會根據(jù)濾波要求確定濾波器指標(biāo)參數(shù)；3. 掌握用IIR和FIR數(shù)字濾波器

9、的MATLAB實現(xiàn)方法，并能繪制濾波器的幅頻特性、相頻特性；4. 通過觀察濾波器的輸入、輸出信號的時域波形及其頻譜，建立數(shù)字濾波的概念?！径吭O(shè)計原理抑制載波單頻調(diào)幅信號的數(shù)學(xué)表達(dá)式為 其中，稱為載波，為載波頻率，稱為單頻調(diào)制信號，為調(diào)制正弦波信號頻率，且滿足>。由上式可見，所謂抑制載波單頻調(diào)制信號，就是兩個正弦信號相乘，它有2個頻率成分：和頻+，差頻-，這兩個頻率成分關(guān)于載波頻率對稱。所以，1路抑制載波單頻調(diào)幅信號的頻譜圖是關(guān)于載波頻率對稱的兩根譜線。顯然，當(dāng)調(diào)制頻率和（或）載波頻率不同時，可以得到包含不同頻率成分的單頻調(diào)幅信號，將幾路不同頻率成分的單頻調(diào)幅信號相加后形成混合信號?！救?/p>

10、】設(shè)計結(jié)果1.（1）程序：>>%產(chǎn)生信號序列st，并顯示st的時域波形和頻譜>>%st=mstg返回三路調(diào)幅信號相加形成的混合信號，長度N=800>>N=800; %信號長度N為800>>Fs=10000;T=1/Fs;Tp=N*T; %采樣頻率Fs=10kHz，Tp為采樣時間>>t=0:T:(N-1)*T;k=0:N-1;f=k/Tp; >>fc1=Fs/10; %第1路調(diào)幅信號載波頻率fc1=1000Hz>>fm1=fc1/10; %第1路調(diào)幅信號的調(diào)制信號頻率fm1=100Hz>>fc2=Fs

11、/20; %第2路調(diào)幅信號載波頻率fc2=500Hz>>fm2=fc2/10; %第2路調(diào)幅信號的調(diào)制信號頻率fm2=50Hz>>fc3=Fs/40; %第3路調(diào)幅信號載波頻率fc3=250Hz>>fm3=fc3/10; %第3路調(diào)幅信號的調(diào)制信號頻率fm3=25Hz>>xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %產(chǎn)生第1路調(diào)幅信號>>xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %產(chǎn)生第2路調(diào)幅信號>>xt3=cos(2*pi*fm3*t).*cos(2*p

12、i*fc3*t); %產(chǎn)生第3路調(diào)幅信號>>st=xt1+xt2+xt3; %三路信號相加，得到復(fù)合信號>>fxt=fft(st,N); %計算信號st的頻譜>>%以下為繪圖命令>>figure(1)>>subplot(4,2,1);>>plot(t,st);grid;xlabel('t/s');ylabel('s(t)');>>axis(0,Tp,min(st),max(st);title('(a)s(t)的波形')>>figure(1)>&g

13、t;subplot(4,2,2);>>stem(f,abs(fxt)/max(abs(fxt),'.');grid;title('(b)s(t)的頻譜')>>axis(0,Fs/8,0,1.2);>>xlabel('f/Hz');ylabel('幅度');低通濾波器：>>%低通濾波器>>fp=300;fs=320;>>wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指標(biāo)（低通濾波器的通、阻帶邊界頻）>> N,wp=el

14、lipord(wp,ws,rp,rs); %調(diào)用ellipord計算橢圓DF階數(shù)N和通帶截止頻率wp>> B,A=ellip(N,rp,rs,wp); %調(diào)用ellip計算橢圓帶通DF系統(tǒng)函數(shù)系數(shù)向量B和A>>y1t=filter(B,A,st); %濾波器軟件實現(xiàn)>>y1=fft(y1t);% 繪圖部分>>figure(2)>>subplot(311)>>H,W=freqz(B,A,800);>>plot(W*Fs/2/pi,abs(H);>>xlabel('Hz');>&g

15、t;ylabel('H(w)');>>title('低通濾波器');>>axis(0,2000,0,1.2);>>grid;>>figure(1)>>subplot(4,2,3);>>plot(t,y1t);>>xlabel('t');>>ylabel('y(t)');>>title('分離出的250Hz的波形');>>figure(1)>>subplot(4,2,4);>>

16、;stem(f,abs(y1)/max(abs(y1),'.');>>xlabel('Hz');>>ylabel('|H|');>>title('250Hz的頻譜');>>vaxis(0,1200,0,1);>>grid;帶通濾波器：%帶通濾波器>>fpl=400;fpu=580;fsl=300;fsu=700;>>wp=2*fpl/Fs,2*fpu/Fs;ws=2*fsl/Fs,2*fsu/Fs;rp=0.1;rs=60; >> N,

17、wp=ellipord(wp,ws,rp,rs); %調(diào)用ellipord計算橢圓DF階數(shù)N和通帶截止頻率wp>> B,A=ellip(N,rp,rs,wp); %調(diào)用ellip計算橢圓帶通DF系統(tǒng)函數(shù)系數(shù)向量B和A>>y2t=filter(B,A,st); %濾波器軟件實現(xiàn)>>y2=fft(y2t);%繪圖部分>>figure(2)>>subplot(312)>> H,W=freqz(B,A,800);>>plot(W*Fs/2/pi,abs(H);>>xlabel('Hz');&

18、gt;>ylabel('H(w)');>>title('帶通濾波器');>>axis(0,2000,0,1.2);>>grid;>>figure(1)>>subplot(4,2,5);>>plot(t,y2t);>>xlabel('t');>>ylabel('y(t)');>>title('分離出的500Hz的波形');>>figure(1)>>subplot(4,2,6);&g

19、t;>stem(f,abs(y2)/max(abs(y2),'.');>>xlabel('Hz');>>ylabel('|H|');>>title('頻譜');>>axis(0,1200,0,1);>>grid;高通濾波器：%高通濾波器>>fp=800;fs=780;>>wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指標(biāo)（低通濾波器的通、阻帶邊界頻）>> N,wp=ellipord(wp,ws,rp

20、,rs); %調(diào)用ellipord計算橢圓DF階數(shù)N和通帶截止頻率wp>> B,A=ellip(N,rp,rs,wp,'high'); %調(diào)用ellip計算橢圓帶通DF系統(tǒng)函數(shù)系數(shù)向量B和A>>y3t=filter(B,A,st); %濾波器軟件實現(xiàn)>>y3=fft(y3t);%繪圖部分>>figure(2)>>subplot(313)>> H,W=freqz(B,A,800);>>plot(W*Fs/2/pi,abs(H);>>xlabel('Hz');>&g

21、t;ylabel('H(w)');>>title('高通濾波器');濾波分離：>>axis(0,2000,0,1.2);>>grid;>>figure(1)>>subplot(4,2,7);>>plot(t,y3t);>>xlabel('t');ylabel('y(t)');>>title('分離出的1000Hz的波形');>>figure(1)>>subplot(4,2,8);>>s

22、tem(f,abs(y3)/max(abs(y3),'.');>>xlabel('Hz');>>ylabel('|H|');>>title('頻譜');>>axis(0,1200,0,1); grid;（2）調(diào)試結(jié)果：三種濾波器：【四】設(shè)計體會設(shè)計三 語音信號濾波處理【一】設(shè)計目的1. 了解語音信號的產(chǎn)生、采集，能繪制語音信號的頻率響應(yīng)曲線及頻譜圖；2. 學(xué)會用MATLAB對語音信號進(jìn)行分析和處理；3. 掌握用濾波器去除語音信號噪聲的方法，觀察去噪前后的語音信號。【二】設(shè)計原理1.

23、語音信號的采集在MATLAB軟件平臺下，利用函數(shù)wavread（ ）對語音信號采集，并記錄采樣頻率和采樣點(diǎn)數(shù)。將語音信號轉(zhuǎn)換成計算機(jī)能夠運(yùn)算的有限長序列。wavread函數(shù)的調(diào)用格式如下：y=wavread(file)讀取file所規(guī)定的wav文件，返回采樣值放在向量y中。y,fs=wavread(file)采樣值放在向量y中，fs表示采樣頻率（Hz）。 2. 用FFT作譜分析FFT即快速傅立葉變換，它是從DFT運(yùn)算中發(fā)展起來的，利用系數(shù) 的對稱性和周期性減少運(yùn)算量。長度為N的序列直接計算DFT需要次復(fù)乘和次復(fù)加，而用FFT進(jìn)行運(yùn)算一般需要次復(fù)乘和次復(fù)加，從而使DFT的運(yùn)算大大簡化。用FFT對

24、連續(xù)信號進(jìn)行譜分析的步驟如下：Sc(t)LPFA/DDFT 圖3.1 連續(xù)信號譜分析過程引入前置低通濾波器LPF是為了消除或減少時域連續(xù)信號轉(zhuǎn)換成序列時可能出現(xiàn)的頻譜混疊現(xiàn)象。表示時域有限的窗函數(shù)。3.設(shè)計濾波器去除語音信號的噪聲 通過wavread（）函數(shù)將語音信號讀入，通過頻率采樣及fft ( )產(chǎn)生信號，并對之加噪，通過窗函數(shù)法設(shè)計濾波器濾掉該語音信號的噪聲，對比濾波前后的語音波形和頻譜?！救吭O(shè)計結(jié)果1. 采樣：>> fs=20000;>> x,Fs=wavread('C:Records3.wav');>> figure(1)>

25、> stem(x)采樣后的圖：2（1）程序：>> Fs =20000; %采樣頻率>> x,Fs=wavread('C:Records3.wav');>> figure(1)>> subplot(2,2,1)>> plot(x,'r'); %做原始語音信號的時域圖形>> title('原語音信號');grid on;>> xlabel('時間 n');ylabel('音量 n');>> y=fft(x); %做le

26、ngth(x)點(diǎn)的FFT>> subplot(2,2,2)>> y=fft(x); %做length(x)點(diǎn)的FFT>> y=fftshift(y);%平移，是頻率中心為0>> derta_Fs = Fs/length(x);%設(shè)置頻譜的間隔，分辨率>> plot(-Fs/2:derta_Fs: Fs/2-derta_Fs,abs(y);%畫出原始語音信號的頻譜圖>> title('原語音信號的頻譜');grid on;%加正弦噪聲>> y=sin(x);>> y=x+y;>&

27、gt; figure(2)>> plot(y); %做正弦噪聲語音信號的時域圖形>> subplot(2,2,1);>> plot(y); %做正弦噪聲語音信號的時域圖形>> title('加正弦噪聲后語音信號');>> xlabel('時間 n');ylabel('音量 n');>> subplot(2,2,2);>> y1=fft(x); %做length(x)點(diǎn)的FFT>> y1=fftshift(y1);%平移，是頻率中心為0>>

28、derta_Fs = Fs/length(x);%設(shè)置頻譜的間隔，分辨率>> plot(-Fs/2:derta_Fs: Fs/2-derta_Fs,abs(y1);%畫出加正弦噪聲語音信號的頻譜圖>> title('加正弦噪聲后語音信號的頻譜');grid on; %加白噪聲>> y=rand(size(x)/5;>> y=x+y;>> subplot(2,2,3);>> plot(y); %做白噪聲語音信號的時域圖形>> title('加白噪聲后語音信號');>> xlabel('時間 n');ylabel('音量 n');>> subplot(2,2,4);>> y2=fft(x); %做length