南邮通达语音时频参数的提取和应用和语音编码实验.docx

1、南京邮电大学通达学院实验报告实验名称语音时频参数的提取和应用实验语音编码实验课程名称信息处理技术专业综合实验班级学号姓名开课时间2016/2017学年，第二学期实验一语音时频参数的提取和应用一、实验目的熟悉根本的语音时频参数提取方法，观察并比拟各类参数采用不同提取方法的效果，了解这些时频参数在语音处理中的应用。二、实验内容1、编写并调试语音时域参数提取程序，包括短时功率、短时过零率等函数。2、编写并调试语音频域参数提取程序，包括频谱、语谱图、基音频率、共振峰等。3、利用提取出的时频参数对语音信号进展静/清/浊音的划分。三、实验原理语音信号具有短时平稳性，需对其进展分段帧）处理帧长一般取10至3

2、0ms之间）。根据鼓励方式的不同，语音信号可以分成浊音和清音两大类：浊音：鼓励源可模拟成准周期性的脉冲串。该周期称为基音周期，其倒数就是基音频率，基音频率一般在几百赫兹以内。清音：鼓励源可模拟成随机白噪声。简化的语音生成数学模型如下图：对语音时频波形进展绘制与观察，有助于了解语音的清/浊音等特性。对基音周期等声门参数和共振峰等声道参数进展提取与分析，是语音处理与应用的重要内容。对于语音信号，确定了各线性预测系数后，根据H(z)可得其频率响应的估值即LPC谱。5. LPC频谱与基于DFT求得的语音对数幅度谱相比有何不同？LPC谱优点在于可以很好地表示共振峰构造而不出现额外的蜂起和起伏，根本不含有

3、噪声，其频率灵敏度与人耳不相匹配；DET得到的对数幅度谱受基频谐波的影响，最大值只出现在谐波频率上，共振峰测定误差较大。6. LPC合成语音的质量如何？有何改良措施？LPC合成技术的优点是简单直观。其合成过程实质上只是一种简单的解码和拼接过程。另外，由于波形拼接技术的合成基元是语音的波形数据，保存了语音的全部信息，因而对于单个合成基元来说能够获得很高的自然度。但是，由于自然语流中的语音和孤立状况下的语音有着极大的区别，如果只是简单地把各个孤立的语音生硬地拼接在一起，其整个语流的质量势必是不太理想的。而LPC技术从本质上来说只是一种录音+重放，对于合成整个连续语流LPC合成技术的效果是不理想的。

4、因此，LPC合成技术必须和其他技术相结合，才能明显改善LPC合成的质量。实验过程记录与结果分析1、计算短时自相关函数1.1分别计算正弦信号和白噪声的短时自相关函数，估计正弦信号的基音周期。其Matlab代码如下:x=sin(2*pi*0.01*(0:499);r,eta=xcorr(x,l00,unbiased);stem(eta,r);w=randn(500,1);r,eta=xcorr(w,100,unbiased);stcm(cta,r);1.2计算正弦信号叠加白噪声的短时自相关函数，试估计正弦信号的基音周期。其Matlab代码如下:x=sin(2*pi*0.01*(0:499);w=r

5、andn(500,l);x1=x+w;r,eta=xcorr(x1,1()0,unbiased);stem(eta,r);1.3分别画出一帧浊音和一帧清音的语音自相关波形(采样频率为10kHz,帧长为25.6ms,每帧有256个样点)，估计浊音的基音周期。loaddigits;x=digits.three1;m=2756;N=256;n=m-N+1:m;r,eta=xcorr(x(n),250,unbiased);plot(eta,r);m=500;N=256;n=m-N+l:ni;r,eta=xcorr(x(n),250,unbiased);plot(eta,r);2、求取LPC系数2.1加

6、载“digits.threel”语音数据。基于自相关法，求出其中一帧数据采样频率为10kHz,帧长为25.6ms,每帧有256个样点)的14阶LPC系数。画出相应的LPC谱。loaddigits;x=digits.three1;m=2756:N=256;n=m-N+1:m;M=14;|r,cta=xcorr(x(n),250,unbiased);Rx=tocplitz(r(M+1:2*M);rx=r(M+2:2*M+l);a=Rxrx;NFFT=IO24;k=l:NFFT/2;X=fft(x(n).*hann(N),NFFT);Thcta=l./fft(l;-a,NFFT);plot(k,20

7、*logl0(abs(353*Theta(k)X(k)axis(0NFFT/2-infinf)2.2编写一个用以实现杜宾递推算法的函数“durbin.m”,其Matlab代码如下：functiona,xi,kappa=durbin(r,M)kappa=zeros(M,l);a=zeros(M,1);xi=r(l);zeros(MJ);for(j=l:M)kappa(j)=r(j+l)-a(l:j-l)*r(j:-l:2)/xi(j);a(j)=kappa；a(1:j-1)=a(1:j-l)-kappa(j)*a(j-l:-l:1);xi(j+l)=a(l:j-l)-kappa(j)*a(j-l

8、xi(j+l)=xi(j)*(l-kappa(j)A2);end2.3利用函数“durbin.m”,计算2.1中语音数据帧的14阶LPC系数，并与2.1中的结果进展比拟：x=digils.three1;m=2756;N=256;n=m-N+l:m;M=14;Ir,eta=xcorr(x(n),250,unbiased);aLD,xi,kappa=durbin(r(M+l:2*M+1),M);a,aLD,norm(a-aLD)a=0.35120.20720.38080.31630.0522-0.1639-0.2903-0.2898-0.12220.19600.34920.30240.1713-0

9、.8258aLD=0.35120.20720.38080.31630.0522-0.1639-0.2903-0.2898-0.12220.19600.34920.30240.1713-0.8258ans=l. (X)47e-142.4编写一个用以实现反射系数转换为LPC系数的函数“rf21pc.m”,其Matlab代码如下：functiona=rf21pc(kappa)M=length(kappa);a=zcros(MJ);for(j=l:M)a(j)=kappa(j);a(1:j-1)=a(1:j-1)-kappa(j)*a(j-11:1);end2.5编写一个用以实现LPC系数转换为反射系

10、数的函数“lpc2rf.m”，其Matlab代码如下：functionkappa=lpc2rf(a)M=lengih(a);kappa=zeros(M,1);for(j=M:-l:l)kappa=a(j);a(l：j-D=(a(1:j-1)+a(j)*a(j-l:-l:l)/(1-kappa(j)A2);end2.6使用函数“rf21pc.m”和Tpc2rf.m”,检验反射系数和LPC系数相互转换结果。norm(kappa-lpc2rf(aLD)norm(aLD-rf2Ipc(kappa)ans=1.4760e-14ans=3、语音信号的逆滤波3.1利用2.3中求出的一帧语音数据的LPC系数，

11、构造逆滤波器，并画出该帧语音信号的残差波形。ehat=filter(11,1,x(n);plot(x(n)ehal)3.2编写一个用以实现语音信号LPC分析的函数“Ipcaut。，其Matlab代码如下：functionar,xi,e,m=lpcauto(x,M,win,Olap)Nx=lcngth(x);N=length(win);if(N=1)N=win;win=ones(NJ);endF=fix(Nx-Olap)/(N-OIap);ar=zeros(M+l,F);xi=zeros(M+l,F);e=zcros(Nx,l);m=zeros(F,I);n=1:N;nl=1:Olap;n2=N

12、-Olap+1:N;n3=Olap+1:N;win1=win(nl)./(win(n1)+win(n2)+eps);win2=win(n2)./(win(n1)+win(n2)+eps);for(f=l:F)r,eta=xcorr(x(n).*win.M,biased);a,xi(:,f),kappa=durbin(r(M+1:2*M+1),M);ar(:,f)=l;-a;ehat=filtcr(ar(:,f).l,x(n);e(n)=e(n(n1).*win2+ehat(n1).*win1;ehat(n3);%Overlap-add.m(f)=n(N);n=n+(N-Olap);end3.3

13、加载“timitl”语音数据，利用“Ipcauio”函数对该句语音进展LPC分析。画出语音信号及其残差波形。loadlimitI;x=timit1;M=14;N=256;ar,xi,e,m=lpcauto(x,M,hann(N),N/2);subplot(121),plot(x)subplot(122),plot(e)axis(0,60000,-1,1)plot(x,e)soundsc(x)soundsc(c)soundsc(x-e)4、LPC谱估计4.1编写一个用以计算和显示LPC谱的函数Ipcplot”,其Matlab代码如下：functionlpcplot(A,Nfft,Fs,m)M,N

14、=size(A);if(N=l)Theta,F=freqz(l,A,Nfft,Fs);plot(F,20*logl0(abs(Thcta);xlabel(Frequency,itFHz);ylabel(Magnitude,|theta(oniega)|dB);elseif(length(m)=N)error(ThecolumndimensionofAmustbeequaltothelengthofm.)endTheta=zeros(Nfft,N);for(n=l:N)Theta(:,n),F=freqz(1,A(:,n),Nfft,Fs);endMeshHndl=meshz(m,F,20*Io

15、g10(abs(Theta);axisij;view(-45,45);set(MeshHndl,MeshStyle,Column);axistight;axisautoy;axisautoz;xlabcl(SamplcNumber,(itn);ylabelCFrequency,itF)Hz);zlabcl(Magnitudc,|thcta(onicga)|dB);end4.2利用“Ipcplot”函数，画出3.3中语音数据的几帧采样频率为16kHz,帧长为32ms,每帧有512个样点)LPC谱。plot(x(m(71):m(147)lpcplot(ar(:,71:147),512,16000,

16、m(71:147)lpcplot(ar(:,82),512,16000,0)lpcplot(ar(:,ll0),512,16000,0)lpcplot(ar(:,l34),512,16000,0)对于时变语音信号双)，可以计算它的短时功率乙(7)和短时过零率ZA.(/72)o第，”帧帧长为N个样点)的功率计算公式为：过零率计算公式为：对于浊音语音，可以利用其频谱X(C)具有丰富的谐波分量的特点，求出其谐波乘积谱，计算公式为：式中，R般取为5。在谐波乘积谱中，基频分量变得很大，更易于估计基音周期。四、实验方法及程序1.调用mallab中的wavread和wavplay读入并播放语音文件。2.调用

17、fft命令计算语音的频谱。3.调用specgram命令得到语音的语谱图。4.调用plot命令绘制并观察语音的时频域波形。5.根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时功率、短时过零率、基音频率的函数：functionPx=stpowerfunctionZx=stzerocross(x.N)和functionHPSx=hpspectrum(x,N,R)。6.利用短时功率和短时过零率这两个参数可以对语音信号进展浊音/清音分类。编写实现此分类功能的函数voiunvoi.mo五、实验结果与分析1. 清音和浊音的短时功率、短时过零率各有何特点？清音段的能量一般比浊音段的小得多。

18、浊音时能量集中于较低频率段内，具有较低的过零率，而清音时能量集中于较高频率段内，具有较高的过零率。2. 如何利用短时功率和短时过零率这两个参数对语音信号进展浊音/清音的分类？发浊音时，由于声门波引起了谱的高频跌落。所以其语音能量集中在3kHz以下；发清音时，多数能量出现在高频。高频率意味着高的平均过零率，低频率意味着低的平均过零率，-般可认为浊音具有较低的平均过零率，清音具有较高的平均过零率。3. 如何根据时域波形估计各帧元音的基音周期？看图可知四帧语音的周期均为0.01s,第一帧()l()()Hz有10个峰值点，基频约为10Hz,第二帧0l()()Hz有5个峰值点，基频约为2()Hz,第三帧

19、()l()()Hz有5个峰值点，基频约为20Hz,第四帧()l()0Hz有4个峰值点，基频约为25Hzo4. 如何根据基于DFT的对数幅度谱估计某帔浊音的基音周期？01172为一个基因周期。5. 如何根据基于DFT的对数幅度谱，估计出共振峰频率？有7个峰点，因此基频约为167。6. 时域对语音信号进展加窗，反映在频域，其窗谱对基于DFT的对数幅度谱有何影响？如何估计出窗谱的主瓣宽度？因为时间窗幅度的傅立叶变换为脉冲状态，那么频谱偏差很小，当定常过程x(n)的平均值6为零且N很大时ETN(o)=f(a)(l2)所以加窗的周期图TN(3)是频谱密度函数的无偏估计。固有TN(o)=l+顷如果原始信号

20、的频谱成份与FFT中的谱线完全一致，这种情况下采样数据的长度为信号周期的整数倍，频谱中只有主瓣。7. 如何根据窄带、宽带语谱图，提取出语音的基音和共振峰轨迹？在基音提取中，广泛采用语音波形或误差信号波形的低通滤波。提取共振峰特性最简便的手段是使用语谱仪。提取共振峰还有倒谱法、LPC分析法等更为有效、准确的方法。做完试验后，记住要答复这些问题)实验过程记录与结果分析1、进入matlab1.1在Windows环境中，创立一个名为speech的文件夹，存放语音数据和与实验相关的Matlab文件。如：E:speecho用于实验的语音数据(.mat)及相应说明(.txt)包括：digits.matdig

21、its.txt%gliss.mat%两个包含/i/的滑音gliss.txtletters.matletlers.txt%英文字母表中26个字母的发音mal.txtmal-1.matma1-2.matma1-3.matmal-4.mat%一句连续语音%一句连续语音%一句连续语音%一句连续语音timit.txtlimit1.mattimit2.mattimit3.mattimit4.ma(vowels.matvowels.txtwords,matwords.txt%一句连续语音%一句连续语音%一句连续语音%一句连续语音%元音/a/,/i/,/o/,/u/的发音%十个孤立字的发音1.2在Window

22、s桌面上，双击Matlab的图标，进入Matlab的工作环境。在命令窗口提示符后，键入cdE:spcech指令，进入实验文件夹路径。1.3在命令窗口键入clearall指令，去除工作空间的所有变量。2、加载语音数据，熟悉matlab命令2.1加载timitl.mat语音数据，命令为：loadlimitI2.2用命令“who”和“whos列出当前工作空间中的变量。变量“timitl”会以包含55911个元素的矩阵列向量)形式出现Yourvariablesare:timitlNameSizeBytesClassAttributestimitl55911x1447288double2.3利用命令“l

23、ength”和“size”查看语音信号timitl的长度和维数：m=length(timitl)lm,nj=size(timitl)(粘贴实验结果)in=55911m=55911n=I2.4通过以下命令，分别查看语音信号timitl中间的1个、2个和20个数据值:limit1(5001)timitl(5001,5003)timitl(5001:5020)(粘贴实验结果)ans=-0.0646ans=-0.06460.0111ans=-0.06460.03070.0111-0.02000.20210.40120.1869-0.1756-0.19020.08230.23630.17120.0637

24、-0.0209-0.0683-0.03940.00810.01050.03110.06962.5Matlab中，符号/可以实现转置运算，运行以下命令进展观察：size(timill)粘贴实验结果)ans=1559112.6使用函数max和min可以分别得到信号的最大值和最小值：max(timit1)(粘贴实验结果ans=Ians=-0.80683、语音信号的时域分析3.1使用以下命令，画出语音信号timitl(采样频率为16kHz,长度为3.5秒)的时域波形。plot(timitl)title(4imitr)xlabel(TimeIndex,iln)ylabel(Amplitude)axis(

25、0length(tiniitI)-1l|)grid3.2分别画出一帧浊音和一帧清音的语音时域波形(采样频率为16kHz,帧长为25ms,每帧有400个样点)。subplot(211),plot(timitl(14501:14900)%/a/indark.subplot(212),plot(timitl(3550l:35900)%/s/inwash.3.3编写个一用以计算信号短时功率的函数“stpower.m”,其Matlab代码如下：functionPx=stpower(x,N)M=length(x);Px=zeros(M,1);Px(N)=x(1:N)*x(1:N)/N;For(m=(N+l

26、):M)Px(m)=Px(m-1)4-(x(m)A2-x(m-N)A2)/N;3.4编写一个用以计算信号短时过零率的函数“stzerocross.ni,其Matlab代码如下：functionZx=stzerocross(x,N)M=length(x);Zx=zeros(M,l);Zx(N+1)=suni(abs(sign(x(2:N+1)-sign(x(1:N)/(2*N);for(m=(N+2):M)Zx(m)=Zx(m-1)+(abs(sign(x(m)-sign(x(m-l).abs(sign(x(m-N)-sign(x(m-N-1)/(2*N);end3.5加载“digits.mat

27、”语音数据，该数据包含英文数字“0”至U“9”的发音。计算其中的单词“four”的短时功率和过零率(采样频率为10kHz,帧长为30ms,每帧有300个样点)。loaddigits;N=300;x=digits.fourl;Px=stpower(x.N);Zx=stzerocross(x,N);plot(Px*le-5Zxx/20001)3.6加载vowels.mat”语音数据,该数据包含元音/a/,/i/o/,/u/的发音。分别画出一帧/a/、一帧/i/、一帧/o/和一帧/u/的时域波形采样频率为10kHz,帧长为30ms,每帧有300个样点)。loadvowelssubplot(221)p

28、lot(vowels.a_1(2001:2300)subplot(222)plot(vowels.i_1(2001:2300)subplot(223)plot(vowels.o_1(2001:2300)subplot(224)plot(vowels.u_1(2001:23OO)3.7编写一个利用短时功率和短时过零率这两个参数对语音信号进展浊音/清音分类的函数voiunvoi.m”,其Matlab代码如下：functionvoi=voiunvoi(x,N,Pth,Zth)%Short-timepowerPxandzerocrossingZxmeasures.Px=stpower(x,N);Zx=

29、stzerocross(x,N);%Compareestimateswiththresholdvalues.voi=(PxPth*max(Px)&(ZxZth);%Shiftthevoi-flagN/2samplestotheleft(middlesampleinwindow).voi=voi(fix(N/2)+1:length(voi);voi(length(voi)*ones(fix(N/2),1);4、语音信号的频域分析4.1加载“mal_l”语音数据。基于DFT变换，画出其中一帧数据采样频率为8kHz,帧长为37.5ms,每帧有30()个样点)的频域波形(对数幅度谱)。loadmal_

30、l;x=mal_l(4161:4460);subplot(121),plot(x)N=1024;k=0:N/2-l;X=fftshift(fft(x.*hann(length(x).N);subplot。22),plot(k,20*logl0(abs(X(N/2:-1:1),axis(0N/2-1-infinf)4.2编写一个用以计算信号谐波乘积谱的函数“hpspectrum.m,其Matlab代码如下：functionHPSx=hpspcctrum(x,N,R)k=l:R:N/2;K=length(k);X=fft(x.*hann(length(x),N);HPSx=X(k);for(r=H

31、PSx=HPSx.*X(l:r:r*K);end4.3利用函数“hpspectrum.m”,来估计一帧浊音的基音周期。代码如下：X=ma1_1(4161:4460);N=1024;R=5;HPSx=hpspectrum(x,N,R);plot(20*logl0(abs(HPSx)4.4加载“vowels.mat”语音数据，分别画出一*帧/i/和一帧/u/采样频率为10kHz,帧长为30ms,每帧有300个样点)的基于DFT的对数幅度谱。其Matlab代码如下：loadvowelsx=vowels.i_l(2001:2300);N=1024;k=-N/2:N/2-l;X=fftshift(fft

32、(x.*hann(length(x),N);plot(k,20*log!0(abs(X),axis(0fix(N/2)0100J)x=vowels.u_l(2001:2300);N=1024;k=-N/2:N/2-l;X=fftshift(fft(x.*hann(length(x),N);plot(k,20*logl0(abs(X),axis(0fix(N/2)0100)4.5画出一帧清音语音的基于DFT的对数幅度谱。语音数据为mal.l中的第15701-15860个样点采样频率为8kHz,帧长为20ms,每帧有160个样点)。loadma1_1;x=mal-l(4161:4460);plot

33、(x)N=1024;k=-N/2;N/21;X=fftshift(fft(x.*hann(length(x),N);plot(k,20*log!0(abs(X),axis(0fix(N/2)-infinf)4.6加载“timitl”语音数据。利用函数specgram,画出该句语音的语谱图。其Matlab代码如下:loadtimitl;NFFT=256;Fs=16000;Win=256;Noverlap=128;specgram(timitl,NFFT,Fs,Win,Noverlap);4.7加载“gliss语音数据。画出该句语音的语谱图。loadgliss;NFFT=256;Fs=10000;

34、Win=256;Noverlap=128;specgram(gliss.i_2,NFFT,Fs,Win,Noverlap);4.8基于浊音/清音分类函数voiunvoi.m和谐波乘积谱函数hpspectrum.m,编写一个用以计算浊音短时基音周期的函数“slpilch.m”。其Matlab代码如下：functionFp=stpitch(x,N,Pth,Zth,NFFT,R,Fs)M=length(x);N=2*fix(N/2);win=hann(N);Fp=zeros(M,l);voi=voiunvoi(x,N,Pth,Zth);invoi=0;for(m=N:N/2:M)n1-m-N+l:m

35、-N/2;n2=m-N/2+1:m;n=nln2;if(any(*voi(n)if(invoi)Fp(n1)=Fmax;invoi=0;endelseHPSmax,Fmax=max(abs(hpspectrum(x(n),NFFT,R);if(invoi)Fp(n)=Fp(n)+win*Fmax;elseFp(n)=Fmax*ones(N/2,1)win(N/2+l:N);invoi=1;endendendwsave=warning;warning(off);Fp=(Fp*Fs/NFFT).*(Fp./Fp);warning(wsave);1. 调用matlab中的wavread和wavplay读入并播放语音文件。2. 调用fft命令计算语音的频谱。3. 调用speegram命令得到语音的语谱图。4. 调用plot命令绘制并观察语音的时频域波形。5. 根据短时功率、短时过零率、谐波乘枳诺的计算公式，编写分别用以计算短时功率、短时过零率、基音频率的函数：functionPx=stpowcr、functionZx=stzerocross(x,N)和functionHPSx=hpspectr