LPC系数预测
实验目的
语音线性预测的基本思想是：一个语音信号的抽样值可以用过去若干个取样值的线性组合来逼近。通过使实际语音抽样值与线性预测抽样值的均方误差达到最小，可以确定唯一的一组线性预测系数。
本实验要求掌握LPC原理，利用自相关法，将语音序列加窗，然后对加窗语音进行LP分析，编写程序求12阶线性预测系数。
实验原理
1、线性预测编码LPC算法
由于语音样点之间存在相关性，所以可以用过去的样点值来预测现在或未来的样点值，
从而可以通过使实际语音和线性预测结果之间的误差在某个准则下达到最小值来决定唯一的一组预测系数。而这组系数就能反映语音信号的特性，可以作为语音信号特征参数来用于语音编码、语音合成和语音识别等应用中去。
假设y(n)是一实数据列，???n???，我们可以用过去时刻的N个数据来预测当前时刻的数据y(n)， 即：
?(n)???aN(k)y(n?k) , -??n??? yk?1N （1)
这里aN(k)即为预测系数。定义预测误差e(n)为
?(n) e(n)?y(n)?y (2)
我们将采用最小均方误差准则来选择aN(k)的值，使得式(3)总误差EN最小。
N??2EN??e(n)???y(n)??aN(k)y(n?k)?
n???n????k?1???2 (3)
这种优化参数aN(k)的方法导致了求解如下的正则方程组
?ak?1NN(k)Ryy(k?l)??Ryy(l) , l?,2,?,N
(4)
这里的Ryy(k)是序列y(n)的自相关系数。式（4）可写成如下的矩阵形式：
RN?aN??rN
其中
(5)
aN??aN(1),aN(2),?,aN(N)?
T
(6) (7)
TrN?Ryy(1),Ryy(2),?,Ryy(N)
??注意到Ryy(k)具有Ryy(i?j)?Ryy(i?j)的性质，式(5)中的RN可写成如下形式
页 第 1
Ryy(1)?Ryy(N?1)??Ryy(0)?R(1)?R(0)?R(N?2)yyyyyy? (8) RN??????????R(N?1)R(N?2)?R(0)??yyyy?yy?2、Levinson-Durbin算法
Levinson-Durbin算法是求解正则方程组中的预测系数aN的有效算法。这种算法利用了
自相关矩阵中特殊的对称性。注意到RN(i,j)?Ryy(i?j)，即对角线上的元素都相等，所以这个自相关矩阵是Toeplitz矩阵。
Levinson-Durbin 算法利用了Toeplitz矩阵的特点来进行迭代计算。即首先由一阶预测器(N?1)开始，计算预测系数a1(1)。然后增加阶数，利用低阶的结果得到下一个高阶的计算结果。根据(4)式求解得到的一阶预测器的预测系数a1(1)是： a1(1)??Ryy(1)Ryy(0) (9)
其最小均方误差是：
E1?Ryy(0)?2a1(1)?Ryy(1)?a1(1)Ryy(0)?Ryy(0)1?a1(1)2
这里a1(1)?2??(10)
?1是格形滤波器的第一反射系数。
下一步是求解二阶预测器的系数a2(1)和a2(2)，并将结果用a1(1)表示，根据式(5)得到的两个方程是:
a2(1)Ryy(0)?a2(2)Ryy(1)??Ryy(1)
a2(1)Ryy(1)?a2(2)Ryy(0)??Ryy(2)
通过用(9)的解来消去Ryy(1)，我们得到解： a2(2)??
（11）
Ryy(2)?a1(1)?Ryy(1)Ryy(0)1?a12(1)??
?
Ryy(2)?a1(1)?Ryy(1)E1
a2(1)?a1(1)?a2(2)a1(1)
（12）
这样我们得到了二阶预测器的预测系数，我们再次注意到a2(2)??2是格形滤波器中的第二反射系数。
据此类推，我们可以用(m?1)阶预测器的预测系数来表示m阶预测器的系数。这样，我们可将m阶预测系数矢量am写成两矢量的和，也就是
?am(1)??a(2)?a?m?1??dm?1?m??am??????? ????0???m????am(m)?的。我们将m?m自相关矩阵Rm分区如下：
(13)
这里am?1矢量是第(m?1)阶预测器的预测系数，(m?1)维的矢量dm?1和标量?m是待定
页 第 2
?Rm?1Rm??bt?rm?1?b?rm?1? Ryy(0)?? (14)
b
表示
bbtbT这里rm，rm),Ryy(m?2),?,Ryy(1)?(rm?1的上标?1?Ryy(m?1?1)??Trm?1?Ryy(1),Ryy(2),?,Ryy(m?1)元素的倒序排列。
??根据式(13)和(14)，式(5)可以写成如下形式
?Rm?1?bt??rm?1b???am?1??dm?1??rm?rm?1??1?????????R(m)? ?????Ryy(0)??yy????0??m?? (15)
这是Levinson-Durbin 算法中的关键一步，从式(15)中我们得到两个方程
bRm?1am?1?Rm?1dm?1??mrm?1??rm?1
(16) (17)
btbtrma?r?1m?1m?1dm?1??mRyy(0)??Ryy(m)
由于Rm?1am?1??rm?1， 由式(16)得到
1bdm?1???mR?m?1rm?1
(18)
b又由于rm?1仅是rm?1倒序排列，且Rm?1是Toeplitz，因此可得
b Rm?1abm??rm
(19)
即：
?1babm??Rm?1?rm
(20)
因此式(18)可写成
dm?1??mabm?1?am?1(m?1)??a(m?2)?? ??m?m?1?????a(1)m?1?? (21)
现在可用式(17)这个方程来求解?m，如果我们用式(21)来消去式(17)中的dm?1，可得
btbbt?m?Ryy(0)?rm?1am?1??rm?1am?1??Ryy(m)
btbt注意到rm?1?am?1?rm?1?am?1，由式(22)可得
（22）
?m??btRyy(m)?rm?1am?1tRyy(0)?rm?1am?1 (23)
因此，通过用式(21)和(23)的结果，替换式(13)中的值，我们得到了求解预测器系数的
Levinson-Durbin 迭代算法。
实验内容
有一段语音信号，采样率为8kHz；加长度为120样点的汉宁窗。编写程序，求12阶LPC系数。
实验方法及结果
页 第 3
本实验开发工具为Microsoft Visual Studio 2010，采用语音信号为es01_8k_mono.snd,其信号的采样率为8kHz，对其信号数据进行预处理，采用120个采样点长度的汉宁窗进行窗处理，并取整个语音文件中的其中一帧作为实验对象，用C语言编写程序，实现Lenvinson-Durbin算法，从而求解语音序列的12阶LPC系数。
VS2010程序运行结构如下：
实验代码
#include
float *sample,sample_H[N]; float R0,*R; float E0,EP[M]; float K[M]; float A[M][M]; float Coefficient[M+1];
float sum;
页 第 4
}
int i,j; FILE *fp;
fp=fopen(\C:\\\\Users\\\\Sun\\\\Desktop\\\\新建文件夹(2)\\\\es01_8k_mono.snd\,\); sample=Read_Data(fp); for(i=0;i R0=*correlation(sample_H); R=(correlation(sample_H)+1); E0=R0; K[0]=R[0]/R0; A[0][0]=K[0]; EP[0]=(1-K[0]*K[0])*E0; for(i=1;i Coefficient[0]=1.0; printf(\系数:\\n\); printf(\,Coefficient[0]); for(j=1;j<=M;j++) { } printf(\,__FILE__); system(\); return 0 ; Coefficient[j]=-A[j-1][M-1]; printf(\,j,Coefficient[j]); if((j+1)%3==0) printf(\); sum=0; for(j=0;j<=i-1;j++) { } K[i]=(R[i]-sum)/EP[i-1]; A[i][i]=K[i]; EP[i]=(1-K[i]*K[i])*EP[i-1]; for(j=0;j<=i-1;j++) { } A[j][i]=A[j][i-1]-K[i]*A[i-j-1][i-1]; sum+=A[j][i-1]*R[i-j-1]; sample_H[i] = sample[i]*(float)( 0.5-0.5*cos(2*PI*i/(N-1)) ); 页 第 5 Sample.h #define PI 3.1415926 #define M 12 #define N 120 float *correlation(float a[]) { } float *Read_Data(FILE *fp_speech) { } short data; float static a[N]; for(int i=0;i<=N;i++) { } return a+1; fread(&data,sizeof(short),1,fp_speech); a[i]=(float)data/(float)(1024*32); float static relate[M+1]; int i,k; float sum; for(k=0;k<=M;k++) { } return relate; sum=0; for(i=k;i relate[k]=sum; sum+=a[i]*a[i-k]; 页 第 6