G729 之 LPC 线性预测及量化

1. 求LPC 

看了数字语音信号处理中的LPC 自相关求解方式(Levinson-Durbin)，写下自己的理解过程:

 LD 求解是基于最小圴方差  主要是，

     1. 根据最小均方差推出，误差信号与取样信号为正交

     2. 根据正交性质推出： 当前自关值 =   (以前自相关值 * 预测系数 ) 积加和

    

for(int i=1; i<=lpc_order; i++)
{
double k = 0;
double t0 = 0;

for (int j=1; j<i; j++)
{
t0 += a[i-j] * pR[j];
}
t0 += pR[i];

//求出能量比  t0: 实际能量  Alpha 上次预测的差值能量  k 为二个能量比值
k = - t0/Alpha;

//跟据能量比K  重新计算系数 (实际上是把预测的系数修正成为准确的系数, 从而保证，误差为 0)
for (j=1; j<i; j++)
{
     An[j]= a[j] + k*a[i-j];
}
An[i]= k;  

//计算 下一级 预测的能量
Alpha = (1 - k*k)*Alpha;

// 把系数保存
for (j=1; j<=i; j++)
{
      a[j] = An[j];
}

}

for (int j=1; j<=lpc_order; j++)
{
pAz[j] = a[j];
}

 2. LPC 转化为 LSP 

void lpc_Levinson(int *pR, double *pAz, double *pRc, int *pErr, int lpc_order)
{
double Alpha = 0;
double a[LPC_ORDER + 1];
double An[LPC_ORDER + 1];
Alpha = pR[0];

for(int i=1; i<=lpc_order; i++)
{
double k = 0;
double t0 = 0;
for (int j=1; j<i; j++)
{
t0 += a[i-j] * pR[j];
}
t0 += pR[i];

//求出能量比
k = - t0/Alpha;
pRc[i-1] = k;  // 此处为反射系数....

//跟据能量比 重新计算系数 ()
for (j=1; j<i; j++)
{
An[j]= a[j] + k*a[i-j];
}
An[i]= k;

//计算 下一级理论系数
Alpha = (1 - k*k)*Alpha;

for (j=1; j<=i; j++)
{
a[j] = An[j];
}

}

for (int j=1; j<=lpc_order; j++)
{
pAz[j] = a[j];
}

}

3.   lsp 转成 lsf 

void Lsp_lsf(double lsp[],double lsf[], int m )
{

int ind = 63;

for (int i=m-1; i>=0; i--)
{
double tmp = -1.0;
while (1)
{
tmp = lsp[i] - table[ind]/32768.0;
if(tmp < 0)
break;
ind--;
}

double k = (slope[ind]/4096.0);    //查表得,当前的斜率..用来代替弧线....
lsf[i] = ind/128.0 + ( tmp * k ) ;  //  计算归一化频率

// printf("\n %8d %8.4f \n", ind, ind/128.0);
}

}

 4. LSF 的量化 

     1) 计算加权值 

使用公式(并不明白加这个是什么作用，有知道的望告知):

下面为代码实现...

void Get_wegt( double lsf[], double wegt[])
{
int i;
double tmp;

tmp = lsf[1] - PI04 - 1.0;
if (tmp > 0.0)
wegt[0] = 1.0;
else
wegt[0] = tmp * tmp * 10.0 + 1.0;

for ( i=1; i<M-1; i++ )
{
tmp = lsf[i+1] - lsf[i-1] - 1.0;
if (tmp > 0.0)
wegt[i] = 1.0;
else
wegt[i] = tmp * tmp * 10.0 + 1.0;
}

tmp = PI92 - lsf[M-2] - 1.0;

if (tmp > 0.0)
wegt[M-1] = 1.0;
else
wegt[M-1] = tmp * tmp * 10.0 + 1.0;

wegt[4] *= CONST12;
wegt[5] *= CONST12;
}

   2)   MA 预测滤波器及 二级码本量化 

          采用二种不同的MA滤波器 分别计算 LSF的目标矢量

          使用码本一进行一级量化

          使用码本二进行二级量化（这次量化的是一级量化后残差值，并且分前半段，后半段分别 加权 量化 ）

void Lsp_qua_cs( double lsf_in[M], double lsfq_out[M], WORD *code)
{
double wegt[M];

// 计算加权值
Get_wegt(lsf_in, wegt);

{
int k=0;
printf("\n wegt权值 = ");
for (k=0; k<M; k++)
{
printf(" %8.4f ", wegt[k]*2);
}
}

double buff[M];

double rbuf[M] = {0.0};
int cand[MODE] = {0};
int tindex1[MODE] = {0};
int tindex2[MODE] = {0};
double L_tdist[MODE];
int j=0;

for (int mode = 0; mode<MODE; mode++)
{
//使用MA滤波器 -- 取得待量化矢量
Lsp_prev_extract(lsf_in,rbuf, mode);

{
int k=0;

printf("\n rbuf = ");
for (k=0; k<M; k++)
{
printf(" %12.8f ", rbuf[k]);
}
}
//////////////////////////////////////
// 在表一中查找
//////////////////////////////////////
Lsp_pre_select(rbuf, &cand[mode]);
int cand_cur = cand[mode];

//////////////////////////////////////
// 计算前半段
//////////////////////////////////////
Lsp_select_1(rbuf, lspcb1[cand_cur], wegt, &tindex1[mode]);
int index = tindex1[mode];

for(j=0; j<M/2; j++)
{
//前段量化值
buff[j] = lspcb1[cand_cur][j]/8192.0 + lspcb2[index][j]/8192.0;
}
Lsp_expand_1(buff, 0.0012); //-->防止坚锐噪声

//////////////////////////////////////
//计算后半段
//////////////////////////////////////
Lsp_select_2(rbuf, lspcb1[cand_cur], wegt, &tindex2[mode]);
index = tindex2[mode];

for(j=M/2; j<M; j++)
{
//后段量化值
buff[j] = lspcb1[cand_cur][j]/8192.0 + lspcb2[index][j]/8192.0;
}
Lsp_expand_2(buff, 0.0012); //-->防止坚锐噪声
Lsp_expand_1_2(buff, 0.0006); //-->防止坚锐噪声

// printf("\n cand[mode] = %6d index = %6d \n" ,cand[mode], index);

Lsp_get_tdist(wegt, buff, &L_tdist[mode], rbuf, fg_sum[mode]);

printf("\n============= cand = %6d index1 = %6d index2 = %6d [%10.4f] ===========\n" , cand_cur, tindex1[mode], tindex2[mode], L_tdist[mode]);

}

int mode_index = L_tdist[0]>L_tdist[1]?1:0;

int cand_cur = cand[mode_index];
int index1 = tindex1[mode_index];
int index2 = tindex2[mode_index];
for(j=0; j<M/2; j++)
{
buff[j] = lspcb1[cand_cur][j]/8192.0 + lspcb2[index1][j]/8192.0;
}
for(j=M/2; j<M; j++)
{
buff[j] = lspcb1[cand_cur][j]/8192.0 + lspcb2[index2][j]/8192.0;
}

Lsp_expand_1_2(buff, 0.0012);
Lsp_expand_1_2(buff, 0.0006);

// 由矢量转化为LSF
Lsp_prev_compose(buff, lsfq_out, fg[mode_index], freq_prev, fg_sum[mode_index]);

{
int k=0;

printf("\n lspq = ");
for (k=0; k<M; k++)
{
printf(" %12.8f ", lsfq_out[k]);
}
}

// 更新LSP 向量
int k;
for ( k = MA_NP-1 ; k > 0 ; k-- )
{
for(int i=0; i<M; i++) freq_prev[k][i] = freq_prev[k-1][i];
}
for(int i=0; i<M; i++) freq_prev[0][i] = buff[i];

// 排序
for(k=0; k<(M-1); k++)
{
if(lsfq_out[k] > lsfq_out[k+1])
{
double tmp = lsfq_out[k+1];
lsfq_out[k+1] = lsfq_out[k];
lsfq_out[k] = lsfq_out[k + 1];
}
}

}

   3)  计算量化后的LSF 

void Lsf_lsp2( double lsf[], double lsp[], int m )
{

for(int i=0; i<m; i++)
{
double freq = lsf[i]/(PI);
int freq_t = (int)(freq*16384);
int ind = freq_t/256;
int offset = freq_t%256;

lsp[i] = table2[ind]/32768.0 + ( (slope_cos[ind]/524288.0) * (offset/4096.0)) ;

// printf("\n ind = %d offset = %d lsp = %8.4f", ind, offset, lsp[i]);
printf("\n freq_t = %d lsf=%8.4f ind = %6d || offset = %6d lsp = %8.4f", freq_t, lsf[i], ind, offset, lsp[i]);

}

}

    至此 对LPC 的量化，全部完成...................

 问题：

1. Get_wegt 计算加权值，这个有什么用？