线性预测编码LPC

LPC 系数预测
实验目的
语音线性预测的基本思想是：一个语音信号的抽样值可以用过去若干个取样值的线性组合来逼近。

通过使实际语音抽样值与线性预测抽样值的均方误差达到最小，可以确定唯一的一组线性预测系数。

本实验要求掌握LPC 原理，利用自相关法，将语音序列加窗，然后对加窗语音进行LP 分析，编写程序求12阶线性预测系数。

实验原理
1、线性预测编码LPC 算法
由于语音样点之间存在相关性，所以可以用过去的样点值来预测现在或未来的样点值，
从而可以通过使实际语音和线性预测结果之间的误差在某个准则下达到最小值来决定唯一的一组预测系数。

而这组系数就能反映语音信号的特性，可以作为语音信号特征参数来用于语音编码、语音合成和语音识别等应用中去。

假设)(n y 是一实数据列，+∞<<∞-n ，我们可以用过去时刻的N 个数据来预测当前时刻的数据)(n y ， 即：
+∞<<∞--=∑=n - , )()()(ˆ1
N
k N k n y k a n y
（1)
这里)(k a N 即为预测系数。

定义预测误差)(n e 为
)(ˆ)()(n y
n y n e -= (2)
我们将采用最小均方误差准则来选择)(k a N 的值，使得式(3)总误差N E 最小。

∑∑∑∞
-∞==∞
-∞=⎥⎦
⎤⎢⎣⎡-+==n N
k N n N k n y k a n y n e E 2
12
)()()()(
(3)
这种优化参数)(k a N 的方法导致了求解如下的正则方程组
N l l R l k R k a
yy N
k yy N
,,2, , )()()(1
=-=-∑=
(4)
这里的)(k R yy 是序列)(n y 的自相关系数。

式（4）可写成如下的矩阵形式：
N N N r a R -=⋅
(5)
其中
[]T
N N N N N a a a )(,),2(),1( =a
(6) []
T
yy yy yy N N R R R )(,),2(),1( =r
(7)
注意到)(k R yy 具有)()(j i R j i R yy yy -=-的性质，式(5)中的N R 可写成如下形式
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡----=)0()2()1()2()0()1()1()1()
0(yy yy yy yy yy yy yy yy yy N R N R N R N R R R N R R R R (8) 2、Levinson-Durbin 算法
Levinson-Durbin 算法是求解正则方程组中的预测系数N a 的有效算法。

这种算法利用了
自相关矩阵中特殊的对称性。

注意到)(),(j i R j i R yy N -=，即对角线上的元素都相等，所以这个自相关矩阵是Toeplitz 矩阵。

Levinson-Durbin 算法利用了Toeplitz 矩阵的特点来进行迭代计算。

即首先由一阶预测器(1=N )开始，计算预测系数)1(1a 。

然后增加阶数，利用低阶的结果得到下一个高阶的计算结果。

根据(4)式求解得到的一阶预测器的预测系数)1(1a 是： )
0()1()1(1yy yy R R a -
= (9)
其最小均方误差是：
[]
212
111)1(1)0()0()1()1()1(2)0(a R R a R a R E yy yy yy yy -=+⋅+=
(10)
这里11)1(ρ=
a 是格形滤波器的第一反射系数。

下一步是求解二阶预测器的系数)1(2a 和)2(2a ，并将结果用)1(1a 表示，根据式(5)得到的两个方程是:
)1()1()2()0()1(22yy yy yy R R a R a -=+
)2()0()2()1()1(22yy yy yy R R a R a -=+
（11）
通过用(9)的解来消去)1(yy R ，我们得到解： []
1
12112)
1()1()2()
1(1)0()1()1()2()2(E R a R a R R a R a yy yy yy yy yy ⋅+=
-⋅+-
=
)1()2()1()1(1212a a a a +=
（12）
这样我们得到了二阶预测器的预测系数，我们再次注意到22)2(ρ=a 是格形滤波器中的第二反射系数。

据此类推，我们可以用)1(-m 阶预测器的预测系数来表示m 阶预测器的系数。

这样，我们可将m 阶预测系数矢量m a 写成两矢量的和，也就是
⎥⎦
⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=--m m m m m m m m a a a ρ110)()2()1(d a
a (13)
这里1-m a 矢量是第)1(-m 阶预测器的预测系数，)1(-m 维的矢量1-m d 和标量m ρ是待定
的。

我们将m m ⨯自相关矩阵m R 分区如下：
⎥
⎥⎦⎤⎢⎢⎣⎡=---)0(11
1
yy bt
m b
m m m R r r R R (14)
这里[]
T
b m yy yy yy bt m R m R m R )
()1(,),2(),1(11--=--=r r ，b m 1-r 的上标b
表示
[]
T
yy yy yy m m R R R )1(,),2(),1(1-=- r 元素的倒序排列。

根据式(13)和(14)，式(5)可以写成如下形式
⎥⎦⎤⎢⎣⎡-=⎪⎪⎭
⎫ ⎝⎛⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣⎡------)(0)0(11111
1m R R yy m m m m yy bt m b
m m r d a r r R ρ (15)
这是Levinson-Durbin 算法中的关键一步，从式(15)中我们得到两个方程
111111-------=++m b
m m m m m m r r d R a R ρ
(16) )()0(1111m R R yy yy m m bt m m bt m -=++----ρd r a r
(17)
由于111----=m m m r a R ， 由式(16)得到
b m m m m 1111-----=r R d ρ
(18)
又由于b
m 1-r 仅是1-m r 倒序排列，且1-m R 是Toeplitz ，因此可得
b
m b m m r a R -=-1
(19)
即：
b
m m b m r R a ⋅-=--11
(20)
因此式(18)可写成
⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎣⎡--==-----)1()2()1(11
11
1m m m m b m m m a m a m a ρρa d (21)
现在可用式(17)这个方程来求解m ρ，如果我们用式(21)来消去式(17)中的1-m d ，可得
[])()0(1111m R R yy m bt
m b m bt m yy m -=++----a r a r ρ
（22）
注意到1111----⋅=⋅m t
m b m bt m a r a r ，由式(22)可得
1
111)0()(----++-
=m t m yy m bt
m yy m R m R a r a r ρ (23)
因此，通过用式(21)和(23)的结果，替换式(13)中的值，我们得到了求解预测器系数的
Levinson-Durbin 迭代算法。

实验内容
有一段语音信号，采样率为8kHz ；加长度为120样点的汉宁窗。

编写程序，求12阶LPC 系数。

实验方法及结果
本实验开发工具为Microsoft Visual Studio 2010，采用语音信号为es01_8k_mono.snd,其信号的采样率为8kHz，对其信号数据进行预处理，采用120个采样点长度的汉宁窗进行窗处理，并取整个语音文件中的其中一帧作为实验对象，用C语言编写程序，实现Lenvinson-Durbin算法，从而求解语音序列的12阶LPC系数。

VS2010程序运行结构如下：
实验代码
#include<stdio.h>
#include<iostream>
#include<math.h>
#include"sample.h"
int main()
{
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;
FILE *fp;
fp=fopen("C:\\Users\\Sun\\Desktop\\新建文件夹(2)\\es01_8k_mono.snd","rb");
sample=Read_Data(fp);
for(i=0;i<N;i++)
{
sample_H[i] = sample[i]*(float)( 0.5-0.5*cos(2*PI*i/(N-1)) );
}
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<M;i++)
{
sum=0;
for(j=0;j<=i-1;j++)
{
sum+=A[j][i-1]*R[i-j-1];
}
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];
}
}
Coefficient[0]=1.0;
printf("LPC系数:\n");
printf("a0 = %f\t",Coefficient[0]);
for(j=1;j<=M;j++)
{
Coefficient[j]=-A[j-1][M-1];
printf("a%d = %f\t",j,Coefficient[j]);
if((j+1)%3==0)
printf("\n\n");
}
printf("\n\n%s\n\n",__FILE__);
system("pause");
return 0 ;
}
#define PI 3.1415926
#define M 12
#define N 120
float *correlation(float a[])
{
float static relate[M+1];
int i,k;
float sum;
for(k=0;k<=M;k++)
{
sum=0;
for(i=k;i<N;i++)
{
sum+=a[i]*a[i-k];
}
relate[k]=sum;
}
return relate;
}
float *Read_Data(FILE *fp_speech)
{
short data;
float static a[N];
for(int i=0;i<=N;i++)
{
fread(&data,sizeof(short),1,fp_speech);
a[i]=(float)data/(float)(1024*32);
}
return a+1;
}。