語音信號(hào)處理實(shí)驗(yàn)報(bào)告 (2)

1、語音信號(hào)處理實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)一1用Matlab讀取一段話音（自己錄制一段，最好其中含有漢語四種聲調(diào)變化， 該段話音作為本課實(shí)驗(yàn)原始材料），繪制原始語音波形圖。2.用Matlab計(jì)算這段語音的短時(shí)平均過零率、短時(shí)平均能量和短時(shí)平均幅度，并將多個(gè)波形同步顯示繪圖。 3.觀察各波形在不同音情況下的參數(shù)特點(diǎn)，并歸納總結(jié)其中的規(guī)律。clc clc;x,fs=wavread('benpao.wav');figureplot(x);axis(0 length(x) min(x) max(x);title('原始語音波形')xlabel('時(shí)間')f=enframe

2、(x,300,100);m,n=size(f);for i=1:menergy(i)=sum(f(i,1:n).2);mn(i)=sum(abs(f(i,1:n);endfigureplot(energy);axis(0 length(energy) min(energy) max(energy);title('短時(shí)能量')figureplot(mn);axis(0 length(mn) min(mn) max(mn);title('短時(shí)幅度')lingd=zeros(m);for x=1:m for y=1:n-1 temp=f(x,y)*f(x,y+1) ;

3、 if temp<= 0 lingd(x)=lingd(x)+1; endend%temp1=num(x,1)/300;%count(x)=temp1;endfigureplot(lingd);%axis(0 length(lingd) min(lingd) max(lingd);title('短時(shí)過零率')子函數(shù)：function f=enframe(x,win,inc)%定義函數(shù)。nx=length(x(:);%x(:)的作用是把x給弄成一個(gè)向量，x為一行，則變成一列，如果為矩陣，則按每一列的順序排成一列。得出的nx為序列的數(shù)據(jù)個(gè)數(shù)。nwin=length(win);

4、if(nwin=1)len=win;%如果win中就一個(gè)數(shù)，則len就=該數(shù)，此例中為256個(gè)點(diǎn)。即每幀長elselen=nwin;%如果有多個(gè)數(shù)，則len=個(gè)數(shù)。endif(nargin<3)%nargin返回的是函數(shù)輸入的個(gè)數(shù)，如果中間有變量，返回的是負(fù)值。inc=len;%也就是說，如果函數(shù)enframe的輸入只有兩個(gè)的話，系統(tǒng)就自動(dòng)賦incendnf=fix(nx-len+inc)/inc);%這個(gè)比較關(guān)鍵，nf為分幀的組數(shù)f=zeros(nf,len);%構(gòu)成以組數(shù)為行，幀長為列的矩陣。indf=inc*(0:(nf-1).'%indf為一列nf個(gè)數(shù)據(jù)，即0到nf-1的

5、inc倍，即分好了每幀起點(diǎn)。inds=(1:len);%構(gòu)成了長度為len的一行。f(:)=x(indf(:,ones(1,len)+inds(ones(nf,1),:);if(nwin>1)%nwin大于1的情況就不說了。w=win(:)'f=f.*w(ones(nf,1),:);end結(jié)論：濁音能量集中與較低頻段中，所以具有較低的過零率；相反，清音能量集中于較高頻率中，所以有較高的過零率。濁音短時(shí)的平均幅度比清音的要高，所以可以從圖中區(qū)別濁音和清音?？梢詮倪吘墔^(qū)別聲母和韻母的分界、無聲和有聲的分界。實(shí)驗(yàn)二要求：1.分別選定濁音段和清音段一幀語音信號(hào)，計(jì)算其自相關(guān)函數(shù)和平均幅度

6、差函數(shù)，并同步繪制兩參數(shù)波形，總結(jié)在不同音情況下，這兩參數(shù)表現(xiàn)的特點(diǎn)。2.利用三電平中心削波自相關(guān)法實(shí)現(xiàn)清、濁音判定，并提取一段話音的基頻參數(shù)，繪制基頻曲線，并將該參數(shù)曲線合并到前面提取的短時(shí)過零率和短時(shí)能量參數(shù)曲線圖中，以示比較。3.加入3點(diǎn)或5點(diǎn)中值平滑處理后再繪制上述基頻曲線與未平滑曲線比較，觀測完整一段話音的基頻包絡(luò)變化情況。Y,Fs,bits=wavread('benpao.wav'); %濁音，取13270-13510個(gè)點(diǎn) %短時(shí)自相關(guān)函數(shù) temp=Y(13271:13510);Rn1=zeros(1,240); for nn=1:240, for ii=1:24

7、0-nn, Rn1(nn) =Rn1(nn)+ temp(ii)*temp(nn+ii); end end figure(1)subplot(4,1,1);jj=1:240;plot(jj, Rn1,'b');title('濁音-自相關(guān)函數(shù)')grid %短時(shí)平均幅度差函數(shù)Yn1=zeros(1,240); for nn=1:240, for ii=1:240-nn, Yn1(nn) =Yn1(nn)+ abs(temp(ii)-temp(nn+ii); endend figure(1) subplot(4,1,2)jj=1:240;plot(jj, Yn1,&

8、#39;b');title('濁音-短時(shí)平均幅度差函數(shù)')grid %清音，取12120-12360個(gè)點(diǎn) %短時(shí)自相關(guān)函數(shù) temp=Y(12121:12360);Rn2=zeros(1,240); for nn=1:240, for ii=1:240-nn, Rn2(nn) =Rn2(nn)+ temp(ii)*temp(nn+ii); end endfigure(1) subplot(4,1,3)jj=1:240;plot(jj, Rn2,'g');title('清音-自相關(guān)函數(shù)')grid %短時(shí)平均幅度差函數(shù)Yn2=zeros(1

9、,240); for nn=1:240, for ii=1:240-nn, Yn2(nn) =Yn2(nn)+ abs(temp(ii)-temp(nn+ii); endend figure(1) subplot(4,1,4)jj=1:240;plot(jj, Yn2,'g');title('清音-短時(shí)平均幅度差函數(shù)')gridfid=fopen('benpao.txt','rt'); a,count=fscanf(fid,'%f',1,inf); L=length(a); m=max(a);for i=1:La(

10、i)=a(i)/m; endm=max(a); n=min(a); ht=(m+n)/2;for i=1:L; a(i)=a(i)-ht;endfigure(1); subplot(2,1,1); plot(a,'k');axis(0,1711,-1,1); title('中心削波前語音波形'); xlabel('樣點(diǎn)數(shù)'); ylabel('幅度'); coeff=0.7; th0=max(a)*coeff; for k=1:L ; if a(k)>=th0 a(k)=a(k)-th0; elseif a(k)<=(

11、-th0); a(k)=a(k)+th0; else a(k)=0; end endm=max(a); for i=1:L; a(i)=a(i)/m;endsubplot(2,1,2); plot(a,'k');axis(0,1711,-1,1); title('中心削波后語音波形'); xlabel('樣點(diǎn)數(shù)'); ylabel('幅度'); fclose(fid); fid=fopen('sara.txt','rt');b,count=fscanf(fid,'%f',1,inf)

12、; fclose(fid);N=320; A=;for k=1:320 sum=0;for m=1:N-(k-1)sum=sum+b(m)*b(m+k-1); endA(k)=sum;endfor k=1:320 B(k)=A(k)/A(1); endfigure(2); subplot(2,1,1); plot(B,'k');title('中心削波前修正自相關(guān)'); xlabel('延時(shí)k'); ylabel('幅度'); axis(0,320,-1,1); N=160; A=;for k=1:160; sum=0;for m=

13、1:N;sum=sum+a(m)*a(m+k-1); endA(k)=sum;endfor k=1:320 C(k)=A(k)/A(1); endsubplot(2,1,2); plot(C,'k');title('中心削波后修正自相關(guān)'); xlabel('延時(shí)k'); ylabel('幅度'); axis(0,320,-1,1); 結(jié)果：主函數(shù)x,fs=wavread('benpao.wav');zhouqi0=medfilt1(x,5);zhouqi2=linsmooth(zhouqi0,5);w=;w=x;

14、w1=w-zhouqi2;w1=medfilt1(w1,5);w1=linsmooth(w1,5);plot(zhouqi2);title('5點(diǎn)中值平滑處理后波形');axis(0,length(x),min(zhouqi2),max(zhouqi2);子函數(shù)function y=linsmooth(x,n,wintype)if nargin<3 wintype='hann'endif nargin<2 n=3;endwin=hann(n);win=win/sum(win);r,c=size(x);if min(r,c)=1; error('

15、;sorry,no matrix here!:(');endif r=1 len=c;else len=r; x=x.'endy=zeros(len,1);if mod(n,2)=0 l=n/2; x=ones(1,l)*x(l) x ones(1,l)*x(len)'else l=(n-1)/2; x=ones(1,l)*x(l) x ones(1,l+1)*x(len)'endfor k=1:len y(k)=win'*x(k:k+n-1);end實(shí)驗(yàn)三要求：1.繪制并觀查清音與濁音幀信號(hào)短時(shí)頻譜的不同。2.繪制并觀察清音和濁音幀信號(hào)倒譜的不同。3.

16、繪制并觀察加不同窗函數(shù)時(shí)上述頻譜的不同。加漢明窗時(shí)的短時(shí)譜cleara=wavread('benpao.wav');subplot(2,1,1);plot(a);title('original signal');gridN=256;k=hamming(N);for m=1:N b(m)=a(m)*k(m);endy=20*log(abs(fft(b,1024);pinlv=(0:1:255)*8000/512;%點(diǎn)和頻率的對應(yīng)關(guān)系subplot(2,1,2);y1(1:256)=y(1:256);plot(pinlv,y1);title('短時(shí)譜'

17、;);xlabel('頻率/Hz')ylabel('對數(shù)幅度/dB')grid清音的發(fā)聲類似于隨機(jī)噪聲，其頻譜具有慢速變化的趨勢，有頻繁的尖峰和谷加矩形窗時(shí)的短時(shí)譜cleara=wavread('benpao.wav');subplot(2,1,1);plot(a);title('original signal');gridN=256;k=linspace(1,1,512);%形成一個(gè)矩形窗，長度為512for m=1:N b(m)=a(m)*k(m);endy=20*log(abs(fft(b,1024);pinlv=(0:1:

18、255)*8000/512;%點(diǎn)和頻率的對應(yīng)關(guān)系subplot(2,1,2);y1(1:256)=y(1:256);plot(pinlv,y1);title('短時(shí)譜');xlabel('頻率/Hz')ylabel('對數(shù)幅度/dB')grid cleara=wavread('benpao');N=300;h=hamming(N);for m=1:N b(m)=a(m)*h(m);endd=rceps(b);d=fftshift(d);plot(d);title('h 加漢明窗時(shí)的倒譜')cleara=wavrea

19、d('benpao');N=300;h=linspace(1,1,512);%形成一個(gè)矩形窗，長度為512for m=1:N b(m)=a(m)*h(m);endd=rceps(b);%MATLAB提供的倒頻譜函數(shù)d=fftshift(d);% fftshift的作用正是讓正半軸部分和負(fù)半軸部分的圖像分別關(guān)于各自的中心對稱plot(d);title('h 加矩形窗時(shí)的倒譜')實(shí)驗(yàn)四要求：1.利用Durbin算法求解LPC線性預(yù)測系數(shù)，進(jìn)而求得反射系數(shù)和預(yù)測器殘差能量，列表輸出相應(yīng)數(shù)值（可求解濁音一幀數(shù)值即可）。2.繪制濁音段和清音段的一幀LPC功率譜，并與對應(yīng)短

20、時(shí)譜圖形比較。3.從LPC系數(shù)求解LPCC參數(shù)，并輸出其值。 close allclearclcx,fs=wavread('benpao.wav');% 預(yù)加重濾波器xx=double(x);y=filter(1 -0.9495,1,xx);N=160;y1=y(1:N);w1=hamming(N);y2=(y1.*w1)'%加窗 取一幀數(shù)據(jù)p=30;%預(yù)測階數(shù)%首先求自相關(guān)函數(shù)r=zeros(1,p+1);for k=1:p+1sum=0;for m=1:N+1-ksum=sum+y2(m).*y2(m-1+k)'endr(k)=sum;end%根據(jù)durbi

21、n算法求線性預(yù)測系數(shù)k=zeros(1,p);k(1)=r(2)/r(1);a=zeros(p,p);a(1,1)=k(1);e=zeros(1,p);e(1)=(1-k(1)2)*r(1);%遞推過程for i=2:pc=zeros(1,i);sum=0;for j=1:i-1sum=sum+(a(i-1,j).*r(i+1-j);endc(i)=sum;k(i)=(r(i+1)-c(i)/e(i-1);%求反射系數(shù)if find(abs(k)>1) disp('default')elsesubplot(413);plot(abs(k);title('|k(i)

22、|<=1')enda(i,i)=k(i);for j=1:i-1a(i,j)=a(i-1,j)-k(i).*a(i-1,i-j);ende(i)=(1-k(i)2)*e(i-1);%預(yù)測器殘差能量subplot(414);plot(e);title('預(yù)測器殘差能量E(i)')end%遞推結(jié)束后提取預(yù)測系數(shù)d=zeros(1,p);for t=1:pd(t)=a(p,t);endz=zeros(1,N);for i=1:pz(i)=y2(i);endfigure(1);subplot(411);plot(y2);title('原始數(shù)據(jù)')subpl

23、ot(412);plot(z);title('durbin算法求線性預(yù)測系數(shù)')waveFile = 'benpao.wav 'startIndex=round(15000);y, fs, nbits = wavread(waveFile);time=(1:length(y)/fs;frameSize=floor(44*fs/1000);endIndex=startIndex+frameSize-1;frame = y(startIndex:endIndex);frameSize=length(frame);frame2=frame.*hamming(lengt

24、h(frame); % ¼Óhamming window%短時(shí)譜spectral=20*log10(abs(fft(frame2);len=length(spectral);% 得到線性預(yù)測濾波器ncoeff=round(2+fs/1000);a=lpc(frame2,ncoeff);% 頻率響應(yīng)h,f=freqz(1,a,512,fs);% 畫圖subplot(3,1,1);plot(time,y); title(waveFile); axis tightylim=get(gca, 'ylim');line(time(startIndex), time(startIndex), ylim, 'color', 'r');line(time(endIndex), time(endIndex), ylim, 'color', 'r');axis xyxlabel('時(shí)間/ s');ylabel('幅度');subplot(3,1,2);xj=(1:len/2)*fs/len;plot(xj,spectral(1:len/2);title('短時(shí)譜')