实验4 基于MATLAB的语音信号LPC分析

实验四基于MATLAB的语音信号LPC分析一、实验目的掌握LPC原理，会利用已学的知识，编写程序估计线性预测系数以及LPC的推演参数。能利用所求的相关参数估计语音的端点、清浊音判断、基因周期、共振峰等。二、实验原理LPC分析基本原理LPC分析为线性时不变因果稳定系统V（z）建立一个全极点模型，并利用均方误差准则，对已知的语音信号s(n)进行模型参数估计。1,2,SnSnSnpSnSn如果利用P个取样值来进行预测，则称为P阶线性预测。假设用过去P个取样值的加权之和来预测信号当前取样值，则预测信号为：1pkkSnank(1)其中加权系数用表示，称为预测系数，则预测误差为:ka1pkkensnSnsnank(2)要使预测最佳，则要使短时平均预测误差最小有：2minEen(3)20,(1)kenkpa(4)令,,ikEsniSnk(5)最小的可表示成：min10,00,pkkak(6)显然，误差越接近于零，线性预测的准确度在均方误差最小的意义上为最佳，由此可以计算出预测系数。通过LPC分析，由若干帧语音可以得到若干组LPC参数，每组参数形成一个描绘该帧语音特征的矢量，即LPC特征矢量。由LPC特征矢量可以进一步得到很多种派生特征矢量，例如线性预测倒谱系数、线谱对特征、部分相关系数、对数面积比等等。不同的特征矢量具有不同的特点，它们在语音编码和识别领域有着不同的应用价值。自相关法在最佳线性预测中，若用下式定义的时间平均最小均方准则代替(3)式的集合平均最小均方准则，即令1201minNpnenN（7）事实上就是短时自相关函数，因而,Rikik,RkESnSnk（8）（9）根据平稳随机信号的自相关性质，可得,,1,2;0,1ikRikipkp（10）由(6)式，可得：min10pkkRaRk（11）011102120RRRPRRRPRPRPR综上所述，可以得到如下矩阵形式：（12）123123nRaRaaRaRp协方差法如果在最佳线性预测中，用下式定义的时间平均最小均方准则代替(3)式的集合平均最小均方准则，则可得到类似的方程：121minNnpenN（13）可以看出，这里的数据段两端不需要添加零取样值。在理论上，协方差法计算出来的预测系数有可能造成预测误差滤波器的不稳定，但在实际上当每帧信号取样足够多时，其计算结果将与自相关法的结果很接近，因而稳定性一般是能够保证的(当然这种方法也有量化效应可能引起不稳定的缺点)。协方差解法的最大优点在于不存在自相关法中两端出现很大预测误差的情况，在N和P相差不大时，其参数估值比自相关法要精确的多。但是在语音信号处理时，往往取N在200左右。此时，自相关法具有较大误差的段落在整个语音段中所占的比例很小，参数估值也是比较准确的。在这种情况下，协方差法误差较小的优点就不再突出，其缺乏高效递推算法的缺点成为了制约因素。所以，在语音信号处理中往往使用高效的自相关法。全极点声道模型能量分析是基于语音信号能量随时间有相当大的变化，将线性预测分析应用于语音信号处理，不仅是为了利用其预测功能，更因为它提供了一个非常好的声道模型。将式(2)所示的方程看成是滤波器在语音信号激励下的输入输出方程，则该滤波器称为预测误差滤波器，其e(n)是输出误差。变换到z域，P阶预测误差滤波器的系统函数为11piiiHzaz（14）可以看出，如果将预测误差e(n)作为激励信号，使其通过预测误差滤波器的逆滤波器H(Z)，即1111piiiHzAZaz（15）则H(Z)的输出为语音信号s(n)，也就是说，H(Z)在预测误差e(n)的激励下可以合成语音。因此，H(Z)被称为语音信号的全极点模型，也称为语音合成器。该模型的参数就是P阶线性预测的预测系数。因为预测误差含有语音信号的基音信息，所以对于浊音，模型的激励信号源是以基音周期重复的单位脉冲；对于清音，激励信号源e(n)是自噪声。LPC如果声道特性H(Z)用式(14)所示的全极点模型表示，有（16）式中，S(z)和I(z)分别为语音信号和激励源的Z变换。对人的听觉来说，浊音是最重要的语音信号。对于浊音，模型的激励信号源e(n)是以基音周期重复的单位脉冲，此时有可得的Z变换S(z)为（17）式中，为P阶线性预测系数。根据倒谱的定义，对具有最小相位特征的语音信号，有。式中，为语音信号的倒谱。将式(16)代入式(17)，并对两边求导，得111pnnnSzHzIzaz1Iz111pnnnSzaz1,2,,iaip1lnnnnSzCzcz（18）nc11111,1nnnknkkcakcaacnpn（19）根据上式即可由线性预测系数通过递推得到倒谱系数，将这样得到的倒谱称为线性预测倒谱系数。结合语音帧能量构成LPC组合参数实验证明，组合参数可以提高系统的识别性能。组合参数虽然可以提高系统的性能，但很显然，无论是在特征参数提取环节，还是在模型训练和模型匹配环节都使运算量有所增加。在特征参数提取环节，要计算一种以上的特征参数。在模型训练和模型匹配环节，由于组合参数特征矢量的维数较多，使运算复杂度有所增加。运算量的增加会使系统的识别速度受到影响。为使运算量问题得到较好的解决，所以可以由LPC参数与语音帧能量构成组合参数，能够在运算量增加不明显的情况下改进系统的性能。语音帧能量是指一帧语音信号的能量，它等于该帧语音样值的平方和。选取与语音帧能量构成组合参数主要有以下考虑：1）语音帧能量是语音信号最基本的短时参数之一，它表征一帧语音信号能量的大小，是语音信号一个重要的时域特征；2）由一帧语音求出的语音帧能量是一个标量值，与其它参量构成组合参数不会使原特征矢量的维数明显增加，特征矢量的维数越少，则需要的运算复杂度越小，另外，获取语音帧能量的运算并不复杂；3）语音帧能量与LPC参数之间的相关性不大，它们反映的是语音信号的不同特征，应该有较好的效果。模型增益GGen模型的激励信号表示为：1piiGensnasni（20）预测误差e(n)如式（2），这样当实际的预测系数与模型系数相等时，有nGen（21）这说明激励信号正比于误差信号，其比例常数等于模型增益G。通常假设误差信号的能量等于输入激励信号的能量，因此可以得到：1122200NNnmmGemmE（22）对于式中的激励信号额e(n)，主要分为浊音和清音两种情况。其中为浊音时，考虑到此时实际的激励信号为声门脉冲，因此可以将激励信号表示为n=0时的单位抽样。为了保证这个假设成立，要求分析的区间应该大致和语音基因周期的长度相等。当语音为清音时，我们假定激励信号e(n)为一个零均值、单位方差的平稳白噪声过程。采用自相关解法时，浊音的模型增益为210pnniniERaRiG(23)清音计算模型增益的公式和浊音相同。实验结果(参考)我们使用的原始语音为“kdt_070”，采样频率为11000Hz,运行程序见附录。在这里我们取第30帧进行观察，线性预测阶数为12，看到图3.1所示的原始语音帧的波形，预测语音帧波形和它们之间预测误差的波形。图3.2为原始语音帧和预测语音帧的短时谱和LPC谱的波形图3.1原始语音帧、预测语音帧和预测误差的波形图3.2原始语音帧和预测语音帧的短时谱和LPC谱的波形这里我们可以改变线性误差的阶数来观察语音帧的短时谱和LP谱的变化情况，如图3.3。图3.3预测阶数对语音帧短时谱和LPC谱的影响从图中可以看出，P越大，LPC谱越能反映出语音短时谱的细节部分，但LPC谱的光滑度随之下降。由于我们的目的只是用LPC谱反映声道综合效应的谱的表示式，而具体的谐波形状是通过激励谱来控制的，因此LPC谱只要能够体现出语音的共振峰的结构和谱包络就可以，因此从计算复杂性的角度分析，预测阶数P应该适中。图3.4是原始语音和预测误差的倒谱波形，我们可以从中计算出原始语音的基音周期。从图中看出两峰值之间的间隔为40点左右，基音周期为40/11000=3.6ms，频率为278Hz左右。图3.4原始语音和预测误差的倒谱波形图3.5给出了原始语音的语谱图和预测语音的语谱图，通过比较发现，预测语音的预测效果还可以，基音频率相差无几。TimeFrequency原始语音语谱图01002003004005006007008009000204060TimeFrequency预测语音语谱图01002003004005006007008009000204060图3.5原始语音的语谱图和预测语音的语谱图MusicSource=wavread('kdt_070');Music_source=MusicSource';N=256;%windowlength，N=100--1000;Hamm=hamming(N);%createHammingwindowframe=input('请键入想要处理的帧位置=');%originiscurrentframeorigin=Music_source(((frame-1)*(N/2)+1):((frame-1)*(N/2)+N));Frame=origin.*Hamm';%ShortTimeFourierTransform[s1,f1,t1]=specgram(MusicSource,N,N/2,N);[Xs1,Ys1]=size(s1);fori=1:Xs1FTframe1(i)=s1(i,frame);endN1=input('请键入预测器阶数=');%N1ispredictor'sorder[coef,gain]=lpc(Frame,N1);%LPCanalysisusingLevinson-Durbinrecursionest_Frame=filter([0-coef(2:end)],1,Frame);%estimateframe(LP)FFT_est=fft(est_Frame);err=Frame-est_Frame;%error%FFT_err=fft(err);subplot(3,1,1);plot(MusicSource)%做原始语音信号的时域图形title('原始语音信号');subplot(3,1,2),plot(1:N,Frame,1:N,est_Frame,'-r');grid;title('原始语音帧vs.预测后语音帧')subplot(3,1,3),plot(err);grid;title('误差');pausefLength(1:2*N)=[origin,zeros(1,N)];Xm=fft(fLength,2*N);X=Xm.*conj(Xm);Y=fft(X,2*N);Rk=Y(1:N);PART=sum(coef(2:N1+1).*Rk(1:N1));G=sqrt(sum(Frame.^2)-PART);A=(FTframe1-FFT_est(1:length(f1')))./FTframe1;%inversefilterA(Z)subplot(2,1,1),plot(f1',20*log(abs(FTframe1)),f1',(20*log(abs(1./A))),'-r');grid;title('短时谱');subplot(2,1,2),plot(f1',