g723源码详细分析-7-lsp反量化与插值

Lsp_Inq Lsp_Int lsp反量化与插值

这两个函数完成的功能是将量化的lsp系数反量化
由于参与量化的lsp系数是4个子帧中的最后一帧,
另外3个子帧的lsp系数是由之前的帧的lsp系数(PrevLsp)
与当前的lsp经过加权比例混合得到的.

这4组lsp被转化成为4组lpc系数后,构成一个滤波器,利用此滤波器
得到残差信号,这样就可以对残差信号进行自适应编码与固定码本编码了

先来看Lsp_Inq,
这个函数对量化后的lsp进行反量化

首先是一个判断
传入的参数Crc(这个是由于上层应用程序传入的)值为1时,认为是出现了丢包,
编码不存在丢包这一说,所以此处认为是不丢包的

然后将lspid拆开,每8个bit代表一个量化的下标索引,然后根据数组下标索引
从码本表里查相应的值,代码片段如下:

/*
* Inverse quantize the 10th-order LSP vector. Each band is done
* separately.
*/ //查表反量化
for ( i = LspQntBands-1; i >= 0 ; i -- ) {

/*
* Get the VQ table entry corresponding to the transmitted index
*/
Tmp = (Word16) ( LspId & (Word32) 0x000000ff ) ;
LspId >>= 8 ;

LspQntPnt = BandQntTable[i] ;

for ( j = 0 ; j < BandInfoTable[i][1] ; j ++ )
Lsp[BandInfoTable[i][0] + j] =
LspQntPnt[Tmp*BandInfoTable[i][1] + j] ;
}

回忆一下编码过程,采用的时类似差值方式来量化的,
所以从码本表里得到的值,要加上直流分量以及
过程与量化相反

首先得到之前被去除直分量的lsp系数
/*
* Subtract the DC component from the previous frame's quantized
* vector
*/ //减去直流分量
for ( j = 0 ; j < LpcOrder ; j ++ )
PrevLsp[j] = sub(PrevLsp[j], LspDcTable[j] ) ;

然后与码本里的值相加
/*
* Generate the prediction vector using a fixed first-order
* predictor based on the previous frame's (DC-free) quantized
* vector
*/ //形成预测向量
for ( j = 0 ; j < LpcOrder ; j ++ ) {
Tmp = mult_r( PrevLsp[j], Lprd ) ;
Lsp[j] = add( Lsp[j], Tmp ) ;
}

再加上直流分量
/*
* Add the DC component back to the previous quantized vector,
* which is needed in later routines
*/ //因为是用减去直流分量的lsp向量量化的,所以加回来
for ( j = 0 ; j < LpcOrder ; j ++ ) {
PrevLsp[j] = add( PrevLsp[j], LspDcTable[j] ) ;
Lsp[j] = add( Lsp[j], LspDcTable[j] ) ;
}

可以看到,这样就"恢复"了lsp，即反量化了,当然反量化值与原值是略有不同的

接下来进行lsp系数稳定性检测调整,
其实就是保证每个lsp系数分得足够开,如果两个lsp系数太过靠紧,就调整它们的值,让它
们离得开一些

从itu的文档,我们看出两个lsp之间距离不得小于31.25hz,
即8000/31.25=256,换算到量化的间隔(2pi被量化成512等于)
也就是量化后的lsp量必须相差2,注意到量化的lsp是扩大2^7倍的,
2即 0x100
如果相邻lsp系数之间的差值小于2,就取平均值，然后各减一
与itu文档对应的就是取平均,各减31.25/2hz
代码比较简单,如下:
/* Check the first frequency */
if ( Lsp[0] < (Word16) 0x180 )//lsc 这个是3,最小值限制
Lsp[0] = (Word16) 0x180 ;

/* Check the last frequency */
if ( Lsp[LpcOrder-1] > (Word16) 0x7e00 )//lsc 这个是252,最大值限制
Lsp[LpcOrder-1] = (Word16) 0x7e00 ;

/* Perform the modification */
for ( j = 1 ; j < LpcOrder ; j ++ ) {

Tmp = add( Scon, Lsp[j-1] ) ;//lsc 在编码时 Scon实际上就是2 对应的就是31.25hz, 8000/31.25=256
Tmp = sub( Tmp, Lsp[j] ) ;
if ( Tmp > (Word16) 0 ) {//lsc相互距离小于31.25hz,取平均,然后一个加1,一个减1(也就是31.25/2)
Tmp = shr( Tmp, (Word16) 1 ) ;
Lsp[j-1] = sub( Lsp[j-1], Tmp ) ;
Lsp[j] = add( Lsp[j], Tmp ) ;
}
}

然后检测:
for ( j = 1 ; j < LpcOrder ; j ++ ) {//lsc 检查是否满足距离都大于等于31.25hz即2
Tmp = add( Lsp[j-1], Scon ) ;
Tmp = sub( Tmp, (Word16) 4 ) ;
Tmp = sub( Tmp, Lsp[j] ) ;
if ( Tmp > (Word16) 0 )
Test = True ;
}

如果循环了10次仍然不满足每个lsp的值间隔大于2,就取之前的lsp值

Lsp_Int lsp插值
将当前解码出来的lsp与之前帧的lsp系数进行一定比例的混合,得到4个子帧的lsp系数
因为语音数据的相关性,itu采用了这种作法来减少传输的数据量
这个函数比较简单,对着itu的文档自然一清二楚,笔者不再详细说明了

这里稍解释一些LsptoA这个函数
顾名思义,它的功能就是完成lsp转成lpc系数,真正构造用于计算残差的滤波器的还是lpc系数
每完成一组lsp插值,就计算出相应的lpc系数
算法比较直观,由于采用的是定点数运算,所以存在很多数值放大缩小的情况,显得比较绕

首先,LsptoA做的是查cos表,得到扩大2^15倍的cos值,因为此时已经得到角度的量化值了,2pi为量化成512等分(15:7)为整数部分(6:0)为小数部分
代码片段如下:
for ( i = 0 ; i < LpcOrder ; i ++ ) {

/*
* Do the table lookup using bits [15:7] of the LSP frequency
*/
j = (int) shr( Lsp[i], (Word16) 7 ) ;// lsc lsp的低7位是小数点之后的,剩余的高位才直接从cos表里面查
Acc0 = L_deposit_h( CosineTable[j] ) ;//扩大2^16倍,总共2^30

/*
* Do the linear interpolations using bits [6:0] of the LSP
* frequency
*/ //lsc 用直线段代替曲线,量化小数点后面的值,并与整数部分的cos值相加
Tmp = sub(CosineTable[j+1], CosineTable[j] ) ;
Acc0 = L_mac( Acc0, Tmp, add( shl( (Word16)(Lsp[i] & 0x007f) ,
(Word16)8 ), (Word16) 0x0080 ) ) ;
//lsc shl 8,由于小数部分本身是放大2^7, 8+7 + 14(tmp的放大倍数) + 1(L_mac) = 30,这样就与Acc0是一个数量级了
Acc0 = L_shl( Acc0, (Word16) 1 ) ;
Lsp[i] = negate( round( Acc0 ) ) ;//lsc 因为- 2cos(w_i),有个负号,所以这里直接取反 最终是扩大了2^15
}

得到cos值之后,就可以根据根的性质进行多项式乘法得到lpc系数了,
P_i(z) = 1 - 2cos(w_i) z^{-1} + z^{-2}这个方程的根 为 cos(w_i)+sin(w_i)j cos(w_i)-jsin(w_i),所以P(z)可以分解为5个这样的二次多项式
执行多项式相乘,即可得到对应的系数

笔者这里简要地描述一下多项式乘法的过程
1 c 1 c = 2cos(w_i)
1 b 1 b = 2cos(w_i+1)
根据多项式乘法则有系数如下
z^0 z^1 z^2 z^3 z^4
1 b 1
c bc c
1 b 1
从这里,可以看出系数存在对称性
当然由于对称性itu最代p[0]是最高次(完全可以这么假设),这在之后的计算要加以注意,否则就会无法理解代码,

根有10个,由于对称性,只要做5次循环,求出前6个系数,后面的5个系数由对称性可以直接推出来
先计算p[0] p[1] p[2]然后,再做3次循环,迭代更新,求出前6个系数,每次循环会多求出一个系数,
并在之后的循环中更新所有之前计算出来的值

通过这样一个实例,可以更好地理解以下的代码,因为系数总是对称的,只需要计算"前一半"就可以了(有奇数个项)
1 2c 1
p[0] p[1] p[2] p[1] p[0]

执行多项式乘法后,则如下
p[0] p[1] p[2] p[1] p[0]
2c*p[0] 2c*p[1] 2c*p[2] 2c*p[1] 2c*p[0]
p[0] p[1] p[2] p[1] p[0]

从这里我们看到 新增加的 p[3] = p[1] + 2c*p[2] + p[1]
而p[2]要更新为 p[2] += (2c*p[1] + p[0])
p[1]更新为 p[1] += (2c*p[0])
p[0]不变

在计算的过程中,每循环一次,缩小1/2(总共三次,所以缩小了1/8),这样做一方面也能防止溢出,所以我们看到循环过程中,p[0]也在不断地被更新,
因为我们看到,是<相乘迭代>的过程,每一环缩小1/2,最终就会缩小1/8

相应的代码片段:
求p[0] p[1] p[2]
P[0] = (Word32) 0x10000000L ;//lsc这个表示1 2^28
P[1] = L_mult( Lsp[0], (Word16) 0x2000 ) ;//lsc lsp[0]是放大2^15,0x2000即2 2^12 总共2^28,因为有个乘2因子
P[1] = L_mac( P[1], Lsp[2], (Word16) 0x2000 ) ;//两个一次项系相加,得到最最终的一次项系数
P[2] = L_mult( Lsp[0], Lsp[2] ) ;//这里是算bc了
P[2] = L_shr( P[2], (Word16) 1 ) ;//正好放大了4倍(L_mult贡献了乘2)
P[2] = L_add( P[2], (Word32) 0x20000000L ) ;//这里做+2 即 2+bc,这里得到的值,是原值的2^28倍

//lsc 这段代码含义同上,不多做解释了
Q[0] = (Word32) 0x10000000L ;
Q[1] = L_mult( Lsp[1], (Word16) 0x2000 ) ;
Q[1] = L_mac( Q[1], Lsp[3], (Word16) 0x2000 ) ;
Q[2] = L_mult( Lsp[1], Lsp[3] ) ;
Q[2] = L_shr( Q[2], (Word16) 1 ) ;
Q[2] = L_add( Q[2], (Word32) 0x20000000L ) ;

循环求之后的3个系数
for ( i = 2 ; i < LpcOrder/2 ; i ++ ) {

/* Compute coefficient (i+1) */
Acc0 = P[i] ;//lsc 这个值已经在之前的循环迭代中缩放到合适的值
Acc0 = L_mls( Acc0, Lsp[2*i+0] ) ;//lsc 2c*p[2],因为缩小1/2, L_mls能把结果值缩小2^-15
Acc0 = L_add( Acc0, P[i-1] ) ;//lsc p[1] + c*p[2]
P[i+1] = Acc0 ;//lsc 是当前循环迭代的 1/2(这么说有点费解,但笔者也找不到更适合的词了,不是真实值的1/2,在当前循环缩放因子的基础上再缩1/2)

//q的计算同上
Acc1 = Q[i] ;
Acc1 = L_mls( Acc1, Lsp[2*i+1] ) ;
Acc1 = L_add( Acc1, Q[i-1] ) ;
Q[i+1] = Acc1 ;

/* Compute coefficients i, i-1, ..., 2 */
for ( j = i ; j >= 2 ; j -- ) {
Acc0 = P[j-1] ;
Acc0 = L_mls( Acc0, Lsp[2*i+0] ) ;//lsc Lsp[2*i+0]是需要乘2的,但要缩小1/2,所以就不用乘了
Acc0 = L_add( Acc0, L_shr(P[j], (Word16) 1 ) ) ;//p[2]缩小1/2
Acc0 = L_add( Acc0, L_shr(P[j-2], (Word16) 1 ) ) ;//p[0]缩小1/2
P[j] = Acc0 ;//lsc 最终是缩小了1/2

//q的计算同上
Acc1 = Q[j-1] ;
Acc1 = L_mls( Acc1, Lsp[2*i+1] ) ;
Acc1 = L_add( Acc1, L_shr(Q[j], (Word16) 1 ) ) ;
Acc1 = L_add( Acc1, L_shr(Q[j-2], (Word16) 1 ) ) ;
Q[j] = Acc1 ;
}

/* Compute coefficients 1, 0 */
P[0] = L_shr( P[0], (Word16) 1 ) ;//这里更新p[0],每次缩小1/2,
Q[0] = L_shr( Q[0], (Word16) 1 ) ;

Acc0 = L_deposit_h( Lsp[2*i+0] ) ;//这里的lsp被扩大了2^30
Acc0 = L_shr( Acc0, (Word16) i ) ;//这里为了把Lsp[2*i + 0]调至与p[0]相同的数量级,因为p[0]本质是1, 也就是2^(30-i)
Acc0 = L_add( Acc0, P[1] ) ;//数据量相同的情况下才能相加,即 p[1] += (c*p[0])
Acc0 = L_shr( Acc0, (Word16) 1 ) ;//缩小至1/2
P[1] = Acc0 ;

Acc1 = L_deposit_h( Lsp[2*i+1] ) ;
Acc1 = L_shr( Acc1, (Word16) i ) ;
Acc1 = L_add( Acc1, Q[1] ) ;
Acc1 = L_shr( Acc1, (Word16) 1 ) ;
Q[1] = Acc1 ;
}

计算到这一步,笔者推断出p q的系数是扩大2^25的

根据A(z) = 1/2 {P(z) (1+ z^{-1}) + Q(z) (1 - z^{-1})}
求出lpc系数

//lsc 注意itu认为 p[0]是最高次(因为对称性,也可以认为是最低次,但itu假设成最高次也完全合理)
//lsc 清楚这点才能顺畅地理解以下的代码(笔者初看时把顺序弄反了,结果怎么看都不对)
for ( i = 0 ; i < LpcOrder/2 ; i ++ ) {
Acc0 = P[i] ;
Acc0 = L_add( Acc0, P[i+1] ) ;//lsc 到这里完成计算P(z) (1 + z^{-1})
Acc0 = L_sub( Acc0, Q[i] ) ;//lsc 这里是计算Q(z)
Acc0 = L_add( Acc0, Q[i+1] ) ;//lsc 这里是计算Q(z) * z^(-1)
Acc0 = L_shl( Acc0, (Word16) 3 ) ;//这里可以看出扩了2^4倍,因为还有一个1/2因子没有乘进去
Lsp[i] = negate( round( Acc0 ) ) ;//29-16=13

//对称性,同时计算别外一半,因为对称性,这里的加减顺序自然是相反的
Acc1 = P[i] ;
Acc1 = L_add( Acc1, P[i+1] ) ;
Acc1 = L_add( Acc1, Q[i] ) ;
Acc1 = L_sub( Acc1, Q[i+1] ) ;
Acc1 = L_shl( Acc1, (Word16) 3 ) ;
Lsp[LpcOrder-1-i] = negate( round( Acc1 ) ) ;
}

//lsc 到这一步,笔者推断出是扩大了2^13倍
//从之后的计算冲激响应也可以看出来,确实是扩大了2^13,因为我们看到冲激响应的脉冲输入是0x04000000L(它代表1,值为2^26),正好是13的两倍
//考虑到信号其中有冲激响应的乘法,这就正好对上了

到此为止,完成了lsp到lpc的转换,看上去很多,但其实都是根据公式平铺直叙,弄清楚放大倍数,代码就迎刃而解了

接下来的工作就是根据lpc构造冲激响应,进行自适应激励码本与固定码本搜索了
笔者将在后继的章节继续分析
林绍川
2011.6.23 于杭州