南京邮电大学实验报告实验名称语音编码实验课程名称信息处理技术专业综合实验班级学号B13011413姓名陈超开课时间2016/2017学年,第一学期实验二语音编码一、实验目的熟悉语音基本压缩编码的方法,观察语音压缩效果,加深对语音线性预测编码(LPC)的理解。二、实验内容1、编写并调试语音LPC参数提取程序。2、编写并调试语音基音周期提取程序。3、编写并调试语音LPC合成程序。三、实验原理语音信号中含有大量的冗余信息,采用各种信源编码技术减除语音信号的冗余度,并充分利用人耳的听觉掩蔽效应,就可以将其编码速率压缩很多倍,而仍能提供可懂语音。LPC声码器是一种比较简单实用的语音压缩方法,其基本原理是:根据语音生成模型,将语音看作激励源通过一个线性时不变系统产生的输出,利用线性预测分析对声道参数进行估值,将求得的线性预测系数,结合基音周期等少量参数进行传输,就可以在接收端利用合成滤波器重构语音信号。线性预测系数的估计方法为:假设语音的当前样值可以用过去的M个语音样值来进行预测MiiMinxaMnxanxanxanx12121~式中ia即为线性预测系数。实际值和预测值之间的均方误差可表示为nMiininxanxnE212要求均方误差总和最小,将E关于ia的偏导数设置为零,可以得到01Miininxanxknx通过采用自相关法、协方差法或格形法求解该方程,即可得到最优的ia。四、实验方法及程序1.调用xcorr命令计算一帧语音的自相关函数。2.调用toeplitz命令形成该帧语音的自相关矩阵。3.调用durbin命令,采用杜宾递推算法计算该帧语音的线性预测系数。4.编写lpcauto.m函数,求取一句语音信号的线性预测系数及预测残差。选择设当的窗函数对语音信号进行分幀。5.编写lpcpitch函数,由残差信号计算该句语音的基音周期。6.编写lpcgain函数,由预测残差能量,求出该句语音的增益。7.编写lpcsyn函数,由该句语音的基音周期、预测残差能量和增益进行LPC合成。五、实验结果与分析1.如何对全极点模型的线性预测参数进行提取?加窗等预处理对提取结果有何影响?2.如何实现线性预测参数的各种表现方式之间的转换?3.使用预测残差求基音周期有何优点?4.如何根据线性预测系数求得LPC频谱?5.LPC频谱与基于DFT求得的语音对数幅度谱相比有何不同?6.LPC合成语音的质量如何?有何改进措施?实验过程记录与结果分析1、计算短时自相关函数1.1分别计算正弦信号和白噪声的短时自相关函数,估计正弦信号的基音周期。其Matlab代码如下:x=sin(2*pi*0.01*(0:499)');[r,eta]=xcorr(x,100,'unbiased');stem(eta,r);w=randn(500,1);[r,eta]=xcorr(w,100,'unbiased');stem(eta,r);-100-80-60-40-20020406080100-0.200.20.40.60.811.21.2计算正弦信号叠加白噪声的短时自相关函数,试估计正弦信号的基音周期。其Matlab代码如下:x=sin(2*pi*0.01*(0:499)');w=randn(500,1);x1=x+w;[r,eta]=xcorr(x1,100,'unbiased');stem(eta,r);-100-80-60-40-20020406080100-0.6-0.4-0.200.20.40.60.811.21.41.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+1:m;[r,eta]=xcorr(x(n),250,'unbiased');plot(eta,r);-250-200-150-100-50050100150200250-3-2-1012345x104-250-200-150-100-50050100150200250-30-20-100102030405060702、求取LPC系数2.1加载“digits.three1”语音数据。基于自相关法,求出其中一帧数据(采样频率为10kHz,帧长为25.6ms,每帧有256个样点)的14阶LPC系数。画出相应的LPC谱。loaddigits;x=digits.three1;m=2756;N=256;n=m-N+1:m;M=14;[r,eta]=xcorr(x(n),250,'unbiased');Rx=toeplitz(r(M+1:2*M));rx=r(M+2:2*M+1);a=Rx\rx;NFFT=1024;k=1:NFFT/2;X=fft(x(n).*hann(N),NFFT);Theta=1./fft([1:-a],NFFT);plot(k,20*log10(abs([353*Theta(k)’X(k)])))axis([0NFFT/2-infinf])05010015020025030035040045050010203040506070802.2编写一个用以实现杜宾递推算法的函数“durbin.m”,其Matlab代码如下:function[a,xi,kappa]=durbin(r,M)kappa=zeros(M,1);a=zeros(M,1);xi=[r(1);zeros(M,1)];for(j=1:M)kappa(j)=(r(j+1)-a(1:j-1)'*r(j:-1:2))/xi(j);a(j)=kappa(j);a(1:j-1)=a(1:j-1)-kappa(j)*a(j-1:-1:1);xi(j+1)=xi(j)*(1-kappa(j)^2);end2.3利用函数“durbin.m”,计算2.1中语音数据帧的14阶LPC系数,并与2.1中的结果进行比较:x=digits.three1;m=2756;N=256;n=m-N+1:m;M=14;[r,eta]=xcorr(x(n),250,'unbiased');[aLD,xi,kappa]=durbin(r(M+1:2*M+1),M);a,aLD,norm(a-aLD)2.4编写一个用以实现反射系数转换为LPC系数的函数“rf2lpc.m”,其Matlab代码如下:functiona=rf2lpc(kappa)M=length(kappa);a=zeros(M,1);for(j=1:M)a(j)=kappa(j);a(1:j-1)=a(1:j-1)-kappa(j)*a(j-1:-1:1);end2.5编写一个用以实现LPC系数转换为反射系数的函数“lpc2rf.m”,其Matlab代码如下:functionkappa=lpc2rf(a)M=length(a);kappa=zeros(M,1);for(j=M:-1:1)kappa(j)=a(j);a(1:j-1)=(a(1:j-1)+a(j)*a(j-1:-1:1))/(1-kappa(j)^2);end2.6使用函数“rf2lpc.m”和“lpc2rf.m”,检验反射系数和LPC系数相互转换结果。norm(kappa-lpc2rf(aLD))norm(aLD-rf2lpc(kappa))3、语音信号的逆滤波3.1利用2.3中求出的一帧语音数据的LPC系数,构造逆滤波器,并画出该帧语音信号的残差波形。ehat=filter([1;-1],1,x(n));plot([x(n)ehat])050100150200250300-600-400-20002004006003.2编写一个用以实现语音信号LPC分析的函数“lpcauto”,其Matlab代码如下:function[ar,xi,e,m]=lpcauto(x,M,win,Olap)Nx=length(x);N=length(win);if(N==1)N=win;win=ones(N,1);endF=fix((Nx-Olap)/(N-Olap));ar=zeros(M+1,F);xi=zeros(M+1,F);e=zeros(Nx,1);m=zeros(F,1);n=1:N;n1=1:Olap;n2=N-Olap+1:N;n3=Olap+1:N;win1=win(n1)./(win(n1)+win(n2)+eps);win2=win(n2)./(win(n1)+win(n2)+eps);for(f=1:F)[r,eta]=xcorr(x(n).*win,M,'biased');[a,xi(:,f),kappa]=durbin(r(M+1:2*M+1),M);ar(:,f)=[1;-a];ehat=filter(ar(:,f),1,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加载“timit1”语音数据,利用“lpcauto”函数对该句语音进行LPC分析。画出语音信号及其残差波形。loadtimit1;x=timit1;M=14;N=256;[ar,xi,e,m]=lpcauto(x,M,hann(N),N/2);plot([x,e])soundsc(x)soundsc(e)soundsc(x-e)0123456x104-1-0.8-0.6-0.4-0.200.20.40.60.814、LPC谱估计4.1编写一个用以计算和显示LPC谱的函数“lpcplot”,其Matlab代码如下:functionlpcplot(A,Nfft,Fs,m)[M,N]=size(A);if(N==1)[Theta,F]=freqz(1,A,Nfft,Fs);plot(F,20*log10(abs(Theta)));xlabel('Frequency,{\itF}[Hz]');ylabel('Magnitude,|\theta(\omega)|[dB]');elseif(length(m)~=N)error('ThecolumndimensionofAmustbeequaltothelengthofm.')endTheta=zeros(Nfft,N);for(n=1:N)[Theta(:,n),F]=freqz(1,A(:,n),Nfft,Fs);endMeshHndl=meshz(m,F,20*log10(abs(Theta)));axisij;view(-45,45);set(MeshHndl,'MeshStyle','Column');axistight;axis'autoy';axis'autoz';xlabel('SampleNumber,{\itn}');ylabel('Frequency,{\itF}[Hz]');zlabel('Magnitude,|\theta(\omega)|[dB]');end4.2利用“lpcplot”函数,画出3.3中语音数据的几帧(采样频率为16kHz,帧长为32ms,每帧有512个样点)LPC谱。plot(x(m(71):m(147)))010002000300040005000600070008000900010000-1-0.8-0.6-0.4-0.200.20.40.60.81lpcplot(ar(:,71:147),512,16000,m(71:147))11.21.41.61.8x1040200040006000-50050SampleNumber,nFrequency,F[Hz]Magnitude,|()|[dB]lpcplo