西工大LPC线性预测编码 实验报告

1、LPC线性预测编码10011206 2012302516邹维韬一、 LPC基本原理基本原理：线性预测编码就是利用过去的样值对新样值进行预测，然后将样值的实际值与其预测值相减得到一个误差信号，显然误差信号的动态范围远小于原始语音信号的动态范围，对误差信号进行量化编码，可大大减少量化所需的比特数，使编码速率降低。设语音信号的样值序列为s(n),n=1,2,n，其中语音信号的当前取样值，即第n时刻的取样值s(n)。而p阶线性预测，是根据信号过去p个取样值的加权和来预测信号当前取样值s(n)，此时的预测器称为p阶预测器。设(n)为s(n)的预测值，则有 (1) 式中，a1,a2,ap，称为线性预测系数

2、，式(1)称为线性预测器，预测器的阶数为p阶。p阶线性预测器的传递函数为 (2)信号s(n)与其线性预测值(n)之差称为线性预测误差，用e(n)表示。则e(n)为 (3)可见，预测误差e(n)是信号s(n)通过具有如下传递函数的系统输出 (4) 图一 LPC误差滤波器如图一所示。称系统A(z)为LPC误差滤波器，设计预测误差滤波器A (z)就是求解预测系数a1,a2,ap，使得预测器的误差e(n)在某个预定的准则下最小，这个过程称为LPC分析。 Levinson-Durbin算法递推求解自相关法求解式据线性预测分析的原理可知，求解p个线性预测系数的依据，是预测误差滤波器的输出方均值或输出功率最

3、小。可称这一最小方均误差为正向预测误差功率Ep，即 (17)已知上式的第二项为0。再将式（3）代入上式可得 （18） 将式（18）与式（12）组合起来可得 （19） 式（19）方程的系数矩阵是对称的，且沿着任一与主对角线平行的斜对角线上的所有元素相等，系数矩阵大小为p×p,这样的矩阵称为Toeplitz(特普利茨)矩阵。式（19）称为Yule-Walker方程，其中R(p)为根据式（9）确定的待分析语音信号s(n)的自相关序列。可见，为了解得线性预测系数，必须先计算出R（k），然后解式（19）即可。经典的线性预测求解方法之自相关法，该方法假定语音信号序列s(n)在间隔0nN-1以外为

4、0；用窗函数从语音序列中截取出选定的序列部分，截取出的序列记为s(0),s(1),s(N-1)。 用对称Toeplitz矩阵的性质，自相关法求解式( 19)可用Levinson-Durbin(莱文森一杜宾)递推算法求解。算法的计算复杂度为O（p2），利用Levinson-Durbin算法递推时，从最低阶预测器开始，由低阶到高阶进行逐阶递推计算。其递推过程如下： （20） （21） （22） （23） （24）由式（20）到式（24）可对i=1,2,p进行递推求解，其最终解为 在上面的一组式子中，i表示预测器阶数，如表示i预测器的第j个预测系数。对p阶预测器，在上述求解预测器系数的过程中，阶数低

5、于p的各阶预测器系数也同时得到。二、 LPC的代码实现 % Levinson's method %初始化矩阵 R=zeros(1,order+1); pc=zeros(1,order); az=zeros(order,order); N=size(frameData,1); %计算自相关系数R的矩阵 for h=1:order+1 R(h)=0; for f=h:N R(h)=R(h)+frameData(f)*frameData(f-h+1); end end %开始计算初值 pc(1)=R(2)/R(1); az(1,1)=pc(1); E=(1-pc(1)2)*R(1); tem

6、p=0; for h=2:order temp=0; for f=1:h-1 temp=temp+az(h-1,f)*R(h-f+1); end %计算反射系数pc与预测系数az pc(h)=(R(h+1)-temp)/E; az(h,h)=pc(h); for f=1:h-1 az(h,f)=az(h-1,f)-pc(h)*az(h-1,h-f); end E=E*(1-pc(h)2); end %得出预测系数矩阵z z=-az(order,:); %计算增益，使用的是pdf中的公式 gain=R(1); for h=1:order gain=gain+z(h)*R(h+1); end ai

7、(:,nframe) = z; % a1,a2,.,a20 G(nframe) =sqrt(6898403*gain); % gain，系数6898403 三、 编码客观评价A、波形与声谱图比较：aslp_zhy_00001.wav与aslp_zhy_00001.wav比较波形图大致吻合，部分细节有所出入，为LPC导致的失真。声谱图吻合。播放效果除了音质有所下降外，无明显失真。B、信噪比与压缩比计算inspeech, Fs1, bits1 = wavread('wavaslp_zhy_00001.wav');outspeech, Fs2, bits2 = wavread('lpcWavaslp_zhy_00001.wav');%计算信噪比duration1 = length(inspeech);duration2 = length(outspeech);up=0;down=0;for i=1:duration1 up=up+inspeech(i)2;endfor i=1:duration2 down=down+(inspeech(i)-outspeech(i)2;endSNR=10*log10(