倒譜計算與分析倒譜計算與分析倒譜計算與分析xxx公司倒譜計算與分析文件編號：文件日期：修訂次數：第1.0次更改批準審核制定方案設計，管理制度《視頻語音處理技術》倒譜計算與分析學院名稱：計算機與信息工程學院專業名稱：計算機科學與技術年級班級：姓名：學號：計算機與信息技術學院綜合性、設計性實驗報告專業：計算機科學技術年級/班級：2011級2012—2013學年第一學期課程名稱視頻語音處理技術指導教師張新明本組成員學號姓名53王彥秋實驗地點計科樓324實驗時間項目名稱倒譜計算與分析實驗類型設計性實驗目的：對語音信號進行同態分析可得到語音信號的倒譜參數。語音的倒譜是將語音的短時譜取對數后再進行IDFT得到的，所以濁音信號的激勵反映在倒譜上是同樣周期的沖激，借此，可從倒譜波形中估計出基音周期。對倒譜進行低時窗選，通過語音倒譜分析的最后一級，進行DFT后的輸出即為平滑后的對數模函數，這個平滑的對數譜顯示了特定輸入語音段的諧振結構，即譜的峰值基本上對應于共振峰頻率，對于平滑過的對數譜中的峰值進行定位，即可估計共振峰。對于倒譜計算與分析的設計實驗可作如下訓練：1、復倒譜的幾種計算方法：2、最小相位信號法和遞歸法；3、基音檢測；4、共振峰檢測。二、實驗儀器或設備：windowsXP下的Matlab編程環境三、總體設計（設計原理、設計方案及流程等）1．復倒譜的幾種計算方法：在復倒譜分析中，z變換后得到的是復數，所以取對數時要進行復對數運算。這時存在相位的多值性問題，稱為“相位卷繞”。設信號為則其傅里葉變換為對上式取復對數為則其幅度和相位分別為：上式中，雖然，的范圍均在內，但的值可能超過范圍。計算機處理時總相位值只能用其主值表示，然后把這個相位主值“展開”，得到連續相位。所以存在下面的情況：（K為整數）此時即產生了相位卷繞。下面介紹幾種避免相位卷繞求復倒譜的方法。最小相位信號法這是解決相位卷繞的一種較好的方法。但它有一個限制條件：被處理的信號想x(n)必須是最小相位信號。實際上許多信號就是最小相位信號，或可以看作是最小相位信號。語音信號的模型就是極點都在z平面單位圓內的全極點模型，或者極零點都在z平面單位圓內的極零點模型。設信號x(n)的z變換為X(z)=N(z)/D(z)，則有根據z變換的微分特性有若x(n)是最小相位信號，則必然是穩定的因果序列。由Hilbert變換的性質可知，任一因果復倒譜序列都可分解為偶對稱分量和奇對稱分量之和:其中這兩個分量的傅里葉變換分別為的傅里葉變換的實部和虛部。從而可得此即復倒譜的性質3，也就是說一個因果序列可由其偶對稱分量來恢復。如果引入一個輔助因子g(n)，上式可寫作其中最小相位信號法求復倒譜原理框圖如下遞歸法這種方法僅限于是最小相位信號的情況。根據z變換的微分特性得對上式求逆z變換，根據z變換的微分特性，有所以設x(n)是最小相位序列，而最小相位信號序列一定為因果序列，所以有由于及可得遞推公式遞歸運算后由復倒譜定義可知如果x(n)是最大相位序列，則變為其中2、基音檢測；語音的倒譜是將語音的短時譜取對數后再進行IDFT得到的，所以濁音信號的周期性激勵反映在倒譜上是同樣周期的沖激。借此，可從倒譜波形中估計出基音周期。一般把倒譜波形中第二個沖激，認為是對應激勵源的基頻。下面給出一種倒譜法求基音周期的框圖及流程圖如下3.共振峰檢測倒譜將基音諧波和聲道的頻譜包絡分離開來。對倒譜進行低時窗選，通過語音倒譜分析系統的最后一級，進行DFT后的輸出即為平滑后的對數模函數，這個平滑的對數譜顯示了特定輸入語音段的諧振結構，即譜的峰值基本上對應于共振蜂頻率，對平滑過的對數譜中的峰值進行定位，即可估計共振峰。原理框圖及流程圖如下。四、實驗步驟（包括主要步驟、代碼分析等）1．倒譜MATLAB實現代碼段clearall;%倒譜[s,fs,nbit]=wavread('');%讀入一段語音b=s';%將s轉置x=b(5000:5399);%取400點語音N=length(x);%讀入語音的長度S=fft(x);%對x進行傅里葉變換Sa=log(abs(S));%log為以e為底的對數sa=ifft(Sa);%對Sa進行傅里葉逆變換ylen=length(sa);fori=1:ylen/2;sal(i)=sa(ylen/2+1-i);endfori=(ylen/2+1):ylen;sal(i)=sa(i+1-ylen/2);end%繪圖figure(1);subplot(2,1,1);plot(x);%axis([0,400,,])title('截取的語音段');xlabel('樣點數');ylabel('幅度');subplot(2,1,2);time2=[-199:1:-1,0:1:200];plot(time2,sa1);%axis([-200,200,,])title('截取語音的倒譜');xlabel('樣點數');ylabel('幅度');運行結果如圖：2.倒譜法求濁音、清音的基音周期functions=p5_2pitchdetectwaveFile='';[y,fs,nbits]=wavread(waveFile);time1=1:length(y);time=(1:length(y))/fs;frameSize=floor(50*fs/1000);%幀長startIndex=round(5000);%起始序號endIndex=startIndex+frameSize-1;%結束序號frame=y(startIndex:endIndex);%取出該幀frameSize=length(frame);frame2=frame.*hamming(length(frame));%加漢明窗rwy=rceps(frame2);%求倒譜ylen=length(rwy);cepstrum=rwy(1:ylen/2);fori=1:ylen/2;cepstrum1(i)=rwy(ylen/2+1-i);endfori=(ylen/2+1):ylen;cepstrum1(i)=rwy(i+1-ylen/2);end%基因檢測LF=floor(fs/500);%基因周期的范圍是70~500HzHF=floor(fs/70);cn=cepstrum(LF:HF);[mx_cepind]=max(cn);ifmx_cep>&ind>LF;a=fs/(LF+ind);elsea=0;endpitch=a%畫圖figure(1);subplot(3,1,1);plot(time1,y);title('語音波形');%axistightylim=get(gca,'ylim');line([time1(startIndex),time1(startIndex)],ylim,'color','r');line([time1(endIndex),time1(endIndex)],ylim,'color','r');xlabel('樣點數');ylabel('幅度');subplot(3,1,2);plot(frame);%axis([0,400,,])title('一幀語音');xlabel('樣點數');ylabel('幅度')subplot(3,1,3);time2=[-199:1:-1,0:1:200];plot(time2,cepstrum1);%axis([-200,200,,])title('一幀語音的倒譜');xlabel('樣點數');ylabel('幅度');運行結果如下圖：倒譜法求濁音的基音周期清音的倒譜3.共振峰檢測程序waveFile='';[y,fs,nbits]=wavread(waveFile);time=(1:length(y))/fs;frameSize=floor(40*fs/1000);%幀長startIndex=round(15000);%起始序號endIndex=startIndex+frameSize-1;%結束序號frame=y(startIndex:endIndex);%取出該幀frameSize=length(frame);frame2=frame.*hamming(length(frame));%加漢明窗rwy=rceps(frame2);%倒譜求ylen=length(rwy);cepstrum=rwy(1:ylen/2);%基音檢測LF=floor(fs/500);HF=floor(fs/70);cn=cepstrum(LF:HF);[mx_cepind]=max(cn);%找到最大的突起的位置%共振峰檢測核心代碼NN=ind+LF;ham=hamming(NN);cep=cepstrum(1:NN);ceps=cep.*ham;%漢明窗formant1=20*log(abs(fft(ceps)));formant(1:2)=formant1(1:2);fort=3:NN%dosomemedianfilteringz=formant1(t-2:t);md=median(z);formant2(t)=md;endfort=1:NN-1ift<=2formant(t)=formant1(t);elseformant(t)=formant2(t-1)*+formant2(t)*+formant2(t+1)*;endendsubplot(3,1,1);plot(cepstrum);title('倒譜');xlabel('樣點數');ylabel('幅度')%axis([0,220,,])spectral=20*log(abs(fft(frame2)));subplot(3,1,2);xj=(1:length(spectral)/2)*fs/length(spectral);plot(xj,spectral(1:length(spectral)/2));title('頻譜');xlabel('頻率/Hz');ylabel('幅度/dB')%axis([0,500,-100,50])subplot(3,1,3);xi=(1:NN/2)*fs/NN;plot(xi,formant(1:NN/2));title('平滑對數幅度譜');xlabel('頻率/Hz');ylabel('幅度/dB');%a