g723源码详细分析-9-自适应码本搜索

各种语音压缩算法的前半部分,即lpc部分是大同小异的,
均是假设声道是一个十阶的全极点滤波器,通过自适应算法，
求出10个lpc系数,然后就是编码残差信号,即激励信号,
各种算法不同之处，在于对激励信号编码方式以及码本结构。
有各种不同的激励信号的编码,
如自适应码本编码,多脉冲编码,等等
最后在解码时，实际上是做一个语音合成，解码出激励源与声道滤
波器.将激励源与滤波器的脉冲响应进行卷积得到语音数据

g723采用的是两级编码,第一级为自适应码本激励,它体现了语音
信号在时域上的相关性,第二级编码分为两种速率.
高速率编码时，采用的是多脉冲激励.
低速率编码时，采用的是伪随机固定码本激励.体现了语音信号的
随机性

Find_Acbk

查找自适应码本
自适应码本是由固定码本经过加权迭代而来的,
它体现的语音数据在时域上的相关性,
具体的作法是,将历史搜索出来的激励码保存在一个数组里,
用每次编码形成的激励源,更新这个数组.
这个数组作为一个码本,在对新帧进行激励源编码时,在基音周期附近
从这个数组搜索出60个最佳激励.

g723采用的五阶预测方式,即每个自适应激励是由历史的五个激励来加权计算出来的

整个函数的大致业务流程说完了,接下来就来看看代码.

首先对基音周期的范围做个限制,
//lsc 偶数帧 Olp被限为 PitchMin, PitchMax-5, 奇数帧时，因为已被限制，不用重复操作了(基音周期只计算了偶数帧)
if ( (Sfc & (Word16)1) == (Word16) 0 ) {
if ( Olp == (Word16) PitchMin )
Olp = add( Olp, (Word16) 1 ) ;
if ( Olp > (Word16) (PitchMax-5) )
Olp = (Word16)(PitchMax-5) ;
}

Get_Rez( Word16 *Tv, Word16 *PrevExc, Word16 Lag )
从历史码本中获取余量信号
最佳自适应码本的搜索是在基音周期附近搜索得到的
也就是往回退Olp个位置附件, CodStat.PrevExc这个全局变量
保存了历史解码来出的激励信号,
Get_Rez 的第一个参数Tv是输出参数,保存返回值数组,因为每个激励是由五个激励
加权得到的,所的Tv的大小是64(最后一个激励由 60-64这五个激励加权得到)
第二个参数 PrevExc 就是CodStat.PrevExc ，它保存了PitchMax个历史的激励源
Lag表示基音周期,也就是取残差信号时,如果基音周期大于64，则顺序取64个,如果小
于64，则在取大于基音周期的那个激励时,会采用取模来计算下标的方式,
如：基音周期如果是 32,则第33个激励又是从"0"开始，说起来比较绕，其实也简单,
就是让激励比基音周期循环取值

Get_Res的主要代码片段如下:
for ( i = 0 ; i < ClPitchOrd/2 ; i ++ )
Tv[i] = PrevExc[PitchMax - (int) Lag - ClPitchOrd/2 + i] ;


for ( i = 0 ; i < SubFrLen+ClPitchOrd/2 ; i ++ )
Tv[ClPitchOrd/2+i] = PrevExc[PitchMax - (int)Lag + i%(int)Lag] ;//lsc i%Lag,这个处理让激励呈周期性

接下来说明一下自适应激励搜索的过程
首先存在一个历史的解码激励,共143个,这个值恰为基音周期的最大值
在基音周期的 -1 0 +1 +2(即Lp-1, Lp Lp + 1, Lp +2,Lp为基音周期) 附近搜索最佳激励.
一个激励信号是由五个连续的历史激励加权形成的,比如
e(i) = A1 * e(i-L) + A2 * e(i-L + 1) + A3 * e(i-L + 2) + A4 * e(i-L + 3) + A4 * e(i-L + 4)
L是当前搜索的位置,取值为基音周期 -1 0 +1 +2
最佳激励搜索是按与目标向量的欧氏距离最短(最小)得出.即:
(tv[i] - e[i])^2,代入展开得到
tv[i]^2 - 2*tv[i]*e[i] + e[i]^2, 第一项在搜索过是不变的,所以只要求后面两项的最小值即可
将e[i]代入展开得
2*tv[i] * (A0 * e(i-L) + A1 * e(i-L + 1) + A2 * e(i-L + 2) + A3 * e(i-L + 3) + A4 * e(i-L + 4))
- (A0 * e(i-L) + A1 * e(i-L + 1) + A2 * e(i-L + 2) + A3 * e(i-L + 3) + A4 * e(i-L + 4))^2

因为有60个元素,为了简化,以下将60个tv[i]记作TV,以此类推5阶的预测分量的矩阵分别为
E[i] i = 0...4,
那么上式可以简化表示成
2*TV[i] * (A0 * E[0] + A1 * E[1] + ... + A4 * E(4))
- (A0 * E[0] + A1 * E[1] + ... + A4 * E[4])^2 ------- <公式1>
将后一项展开
会得到5个平方项,与10个互相关的2二次项,即:
E[0]^2 + E[1]^2 + ... + E[4]^2 + 2 Σ (E[i]*E[j]) --- i=0...4, j=0...4


即搜索出以上表达式的最大值 --- 注意这里的确是求最大值,符号位变了,在itu实现的723算法里,的确
也是这么做的

我们可以看到,这需要计算e[i]的能量(对量2次方项),零时自相关(对应e[i-L + ?]*e[i-L + ?])
tv[i]与e[i]零时自相关

A1 A2 A3 A4 A5 这些取值自一个码本表的,搜索过程中也要搜索码本表(即:AcbkGainTablePtr)


Find_Acbk 这个函数大部分就是在完成这些工作

现在来看一下代码片段

在取出64个历史自适应激励后(因为是五阶,所以多取4个)
计算卷积,
计算分成5轮进行,这里根据每个阶的特点,简化计算。
首先是计算第五阶与冲激响应的卷积,这是按正常的方式计算的
/* Filter the last one using the impulse response */ //lsc 自适应激励与冲激响应卷积
for ( i = 0 ; i < SubFrLen ; i ++ ) {
Acc0 = (Word32) 0 ;
for ( j = 0 ; j <= i ; j ++ )
Acc0 = L_mac( Acc0, RezBuf[ClPitchOrd-1+j], ImpResp[i-j] ) ;
FltBuf[ClPitchOrd-1][i] = round( Acc0 ) ;//lsc 值的缩放为 13 + 1 - 16 = -2
}


然后从第4阶开始,依次更新到第1阶
/* Update all the others */ //lsc 计算出另外4个分量与综合滤波器冲激响应的卷积
for ( i = ClPitchOrd-2 ; i >= 0 ; i -- ) {
FltBuf[i][0] = mult_r( RezBuf[i], (Word16) 0x2000 ) ;//lsc 先放大 13 然后mult_r 13 -15 = -2 对应 y4[0] 冲激响应的第一项肯定是1,看到一冲激响应的计算过程就知道了
for ( j = 1 ; j < SubFrLen ; j ++ ) {
Acc0 = L_deposit_h( FltBuf[i+1][j-1] ) ;//lsc 还原成 13 对应着y5[i]
Acc0 = L_mac( Acc0, RezBuf[i], ImpResp[j] ) ;// 13 + 1 = 14同一个数量级 对应着R4[0] * I[i]
FltBuf[i][j] = round( Acc0 ) ;//14-16 =-2
}
}
这里笔者举一个例子,就很容易理解这段代码了
但前4轮的卷积,为了节省计算量,是在最后一轮的基础上迭代出来的,以第4轮为例
我们可以列出这样的,
第4轮残差 R4[0] R4[1] R4[2] R4[3] ...
第5轮残差 R5[0] R5[1] R5[2] R5[3] ...
其中 R4[1] = R5[0] 以此类推,因为我们是这么取的
第5轮卷积 y5 会像这样
y5[0] = R5[0] * I[0]
y5[1] = R5[0] * I[1] + R5[1] * I[0]
y5[2] = R5[0] * I[2] + R5[1] * I[1] + R5[2] * I[0]

而第4轮卷积y4,会像这样
y4[0] = R4[0] * I[0]
y4[1] = R4[0] * I[1] + R4[1] * I[0]
y4[2] = R4[0] * I[2] + R4[1] * I[1] + R4[2] * I[0]
根据R4[i+1]=R5[i]我们对第4轮的卷积作个替换
y4[0] = R4[0] * I[0]
y4[1] = R4[0] * I[1] + R5[0] * I[0]
y4[2] = R4[0] * I[2] + R5[0] * I[1] + R5[1] * I[0]

对比y5,我们可以看出y4可以更新如下
y4[0] = R4[0] * I[0]
y4[1] = R4[0] * I[1] + y5[0]
y4[2] = R4[0] * I[2] + y5[1]
以此类推出y4[3] 等等,即节省了乘法计算量

这样,就完成了整个卷积的计算过程,即,已经完成了五个阶的E[0],E[1],E[2],E[3],E[4]的计算

计算 TV 与 E[i]点乘,这个对应的是 2 * TV[i] * E[i]项
/* Compute the cross with the signal */ //lsc 计算与自适应激励响应与原始信号的零时自相关
for ( i = 0 ; i < ClPitchOrd ; i ++ ) {
Acc1 = (Word32) 0 ;
for ( j = 0 ; j < SubFrLen ; j ++ ) {
Acc0 = L_mult( Tv[j], FltBuf[i][j] ) ;
Acc1 = L_add( Acc1, L_shr( Acc0, (Word16) 1 ) ) ;//lsc -1 -2 = -3 Tv是缩小1/2的
}
*lPnt ++ = L_shl( Acc1, (Word16) 1 ) ;//lsc 这里在乘2 //lsc -2
}

计算 E[i]^2,即能量,对应搜索公式中的平方项,比搜索公式中看出,这一项不需要乘2的
/* Compute the energies */
for ( i = 0 ; i < ClPitchOrd ; i ++ ) {
Acc1 = (Word32) 0 ;
for ( j = 0 ; j < SubFrLen ; j ++ )
Acc1 = L_mac( Acc1, FltBuf[i][j], FltBuf[i][j] ) ; //lsc -2 -2 +1 = -3
*lPnt ++ = Acc1 ;// lsc 不用乘2，对应式中的 e1^2之类的项
}

最后计算 E[i] * E[j] 项,这项是要乘2的
/* Compute the between crosses */ //lsc这里计算 e1*e2之类的项,排列组合,有10这样的项
for ( i = 1 ; i < ClPitchOrd ; i ++ ) {
for ( j = 0 ; j < i ; j ++ ) {
Acc1 = (Word32) 0 ;
for ( l = 0 ; l < SubFrLen ; l ++ ) {
Acc0 = L_mult( FltBuf[i][l], FltBuf[j][l] ) ;
Acc1 = L_add( Acc1, L_shr( Acc0, (Word16) 1 ) ) ; //lsc -2 -2 = -4
}
*lPnt ++ = L_shl( Acc1, (Word16) 2 ) ;//lsc 这里在乘2 //lsc 此处的含义仍然是乘2 -4+2 =-2
}
}

计算结果被存在 CorBuf 这个数组里,每20个为一组(TV与E[i]贡献5个,平方项5个,交叉乘10(排列组合),共20个)


接下来归一化:
/* Find Max and normalize */
Acc1 = (Word32) 0 ;
for ( i = 0 ; i < Hb*20 ; i ++ ) {
Acc0 = L_abs(CorBuf[i]) ;
if ( Acc0 > Acc1 )
Acc1 = Acc0 ;
}

//lsc 归一化
Exp = norm_l( Acc1 ) ;
/* Convert to shorts */
for ( i = 0 ; i < Hb*20 ; i ++ ) {
Acc0 = L_shl( CorBuf[i], Exp ) ;
CorVct[i] = round( Acc0 ) ;
}
归一化的结果存在CorVct数组里

误差测试,这一块是确定搜索系数码本的边界,前面根据前面误差,来调整增益加权
这一块需要跟后面的结合起来分析,先知道函数功能


这一段相对简单,算出欧式距离最小的(即<公式1>的最大值,保留相应的系数加权码本表的索引)
AcbkGainTable085 , AcbkGainTable170 , 这两个数组,就是系数加权码本表,
从注释很容易看出,每20个为一组,相应的值都事先计算好,直接套到以下的循环里即可
for ( k = 0 ; k < (int) Hb ; k ++ ) {

/* Select Quantization tables */
l = 0 ;
if ( WrkRate == Rate63 ){
if ( (Sfc & (Word16) 1) == (Word16) 0 ) {
if ( (int)Olp-Pstep+k >= SubFrLen-2 )l ++ ;
}
else {
if ( (int)Olp >= SubFrLen-2 ) l ++ ;
}
}
else {
l = 1;
}


sPnt = AcbkGainTablePtr[l] ;

for ( i = 0 ; i < (int) Bound[l] ; i ++ ) {

Acc0 = (Word32) 0 ;
for ( j = 0 ; j < 20 ; j ++ )
Acc0 = L_add( Acc0, L_shr( L_mult(CorVct[k*20+j], *sPnt ++),
(Word16) 1 ) ) ;//lsc 这里将能量与零时相关等20个因子与加权系数码本表里相乘(系数只有5个,经过处理,以20个为一组)

if ( Acc0 > Acc1 ) {//lsc 取最大的,做为自适应激励
Acc1 = Acc0 ;
Gid = (Word16) i ;//lsc 增益加权码本表的索引
Lid = (Word16) k ;//lsc 基音延迟
}
}
}

接下来解码自适应激励,并将其从原始激励中扣除,将得到类似随机信号的残差,用于固定码本搜索
/* Decode the Acbk contribution and subtract it */
Decod_Acbk( RezBuf, PrevExc, Olp, Lid, Gid ) ;

//lsc 从目标向量里扣除自适量激励部分,得到了近似于随机信号的残差,进入固定随机码本搜索
for ( i = 0 ; i < SubFrLen ; i ++ ) {
Acc0 = L_deposit_h( Tv[i] ) ;
Acc0 = L_shr( Acc0, (Word16) 1 ) ;

for ( j = 0 ; j <= i ; j ++ )
Acc0 = L_msu( Acc0, RezBuf[j], ImpResp[i-j] ) ;
Acc0 = L_shl( Acc0, (Word16) 1 ) ;
Tv[i] = round( Acc0 ) ;
}

至此,自适应码本搜索计算完成
总结一下:
首先,自适应码本初始都会0,它是由于固定码本经过反复迭代而来的,
每次对激励编码后,都进行解码,保存成为自适应码本

自适应码本的搜索,是取连续64个历史值,由相邻的5个历史解码激励,加权形成
一个当前的自适应激励,并通过计算欧式距离最小的原则,确定取哪60个历史激励源,
并确定加权系数

最近太忙,都没时间,花了一个多月,才整理完这个函数的分析
笔者将在下一章节分析固定码本搜索,723编码的分析将接近尾声

林绍川
2011.8.23 于杭州