g723源码详细分析-16-舒适噪声生成

Cod_Cng
这个函数,先利用当前子帧的自相关系数和,即:
R[i]= Rsub0[i] + Rsub1[i] + Rsub2[i] + Rsub3[i]
i = 0 ~ 10
Rsub0[i]表示第一个子帧的自相关系数

利用这组自相关系数,通过莱文森德宾算法,估算当前帧的残差能量
即在德宾算法的最后一轮迭代产生的分母,就是残差能量的估值
代码片段
*/ //lsc CodCng.Acf 里保存的是4个子帧自相关系数的和
CodCng.Ener[0] = Durbin(curCoeff, &CodCng.Acf[1], CodCng.Acf[0], &temp);


量化sid帧的增益
如果前一帧是非静音帧,那么只用当前帧的估值残差能量来进行量化,
如果前一帧是静音帧或者sid帧,就取前三个残差能量加权进行估值了.
代码片段:
//lsc 根据最后三个帧的能量,估计增益
/*
* if first frame of silence => SID frame
*/
if(CodCng.PastFtyp == 1) {
*Ftyp = 2;
CodCng.NbEner = 1;
curQGain = Qua_SidGain(CodCng.Ener, CodCng.ShAcf, CodCng.NbEner);
}

else {
CodCng.NbEner++;
if(CodCng.NbEner > NbAvGain) CodCng.NbEner = NbAvGain;
curQGain = Qua_SidGain(CodCng.Ener, CodCng.ShAcf, CodCng.NbEner);

这里的增益其实也是一个大约的估值
估算匠公式见itu文档的A.4.3节的公式.

代码里并没有计算开方,而是将估值公式平方,然后针对平方后的值进行量化

增益估值量化函数,Qua_SidGain,
在Cod_cng中被调用时,nq!=0,走的时else路径 nq=0是在另一处调用时走的路径(Dec_cng,即解码舒适背景音)
我们先分析else路径的代码:
首先计算出这个表达式:
E[i]
2.70375^2 * ------------ i 取值为 1~3
i * 240
itu算法将
2.70375^2
---------- 预先算好,存在一个数组里
i * 240

TAB_LBC.C:约2978行
Word16 fact[4] = {273, 998, 499, 333};//lsc 这里面的值 从998 = 2.70375^2 * 32768 / (1 * 240) ,以此类推

所以循环时,直接用残差(函数参数:Ener)与fact数组对应元素相乘就可以了
代码片段如下:
sh1 = shEner[0];//lsc 将会存最小的能量归一化位移
for(i=1; i<nq; i++) {
if(shEner[i] < sh1) sh1 = shEner[i];
}
for(i=0, L_x=0L; i<nq; i++) {//lsc 根据能量的个数,每个能量乘上一个比例因子再相加
temp = sub(shEner[i], sh1);
temp = shr(Ener[i], temp);
temp = mult_r(fact[nq], temp);//lsc 这个乘法是没有引入缩放的,相应的fact是放大2^15
L_x = L_add(L_x, L_deposit_l(temp));
}
temp = sub(15, sh1);
L_x = L_shl(L_x, temp);//lsc 恢复成原值,因为残差能量本身是缩小2^-15的

这样,就得到了增益估值的原始值了 L_x
接下来量化增益估值,
量化是非均匀的,
L_x 在 (0 , 32) 时, 将有16个量化阶,也就是量化的精度为2 32=2*16
L_x 在 (32 , 96) 时, 将有16个量化阶,也就是量化的精度为4 96=2*16+4*16
L_x 在 (96 , 352) 时, 将有32个量化阶,也就是量化的精度为8 352=2*16+4*16+8*32

量化过程采用了二分查找的方式,最后会做一个类似四舍五入的处理;
代码片段如下:
/* Quantize L_x */
if(L_x >= L_bseg[2]) return(63);//lsc 超出了上限,返回63即可

/* Compute segment number iseg */
if(L_x >= L_bseg[1]) {
iseg = 2;
exp = 4;//lsc 2分查找 32个刻度,只要做4次循环
}
else {
exp = 3;//lsc 2分查找 16个刻度,只要做3次循环
if(L_x >= L_bseg[0]) iseg = 1;
else iseg = 0;
}

iseg_p1 = add(iseg,1);//lsc 这里是解析度 2 4 8
j = shl(1, exp);//lsc 让j处于区间的中间值 8 8 16
k = shr(j,1);//lsc 分别是 4 4 8

/* Binary search in segment iseg */
for(i=0; i<exp; i++) {
temp = add(base[iseg], shl(j, iseg_p1));//lsc 每个阶段,加上对应段的解析度
L_y = L_mult(temp, temp);//lsc 平方乘2
if(L_x >= L_y) j = add(j, k);//lsc 2分查找,跑到上半区间
else j = sub(j, k);//lsc 2分查找,跑到下半区间
k = shr(k, 1);//lsc 区间变小了
}//lsc 至此,基本上得到了量化后的值

temp = add(base[iseg], shl(j, iseg_p1));
L_y = L_mult(temp, temp);
L_y = L_sub(L_y, L_x);
if(L_y <= 0L) {//lsc 量化值偏小,要做一个类似四舍五入的操作
j2 = add(j, 1);
temp = add(base[iseg], shl(j2, iseg_p1));
L_acc = L_mult(temp, temp);
L_acc = L_sub(L_x, L_acc);
if(L_y > L_acc) temp = add(shl(iseg,4), j);
else temp = add(shl(iseg,4), j2);
}
else {//lsc 量化值偏大,要做一个类似四舍五入的操作
j2 = sub(j, 1);
temp = add(base[iseg], shl(j2, iseg_p1));
L_acc = L_mult(temp, temp);
L_acc = L_sub(L_x, L_acc);
if(L_y < L_acc) temp = add(shl(iseg,4), j);
else temp = add(shl(iseg,4), j2);
}
return(temp);//lsc 这里返回的是打包后量化值,2比bit的区间,4bit的解析度


在Cod_Cng函数中,增益估值被保存在curQGain这个局部变量当中,ok
它是做用来做滤波器相似度的一个判别依据

g723 帧类型有三种 0:静音(这种帧应该是不用传的), 1:非音静 2:sid帧
sid帧传递的信息实际为舒适背景音的滤波参数以及增益,
解码方根据sid帧来产生舒适背景音

这里就会涉及到,什么时候sid帧要被生成发送,因为滤波参数以及增益是会改变的
自然curQGain这个值就会成为判别的一个依据
代码片段如下:
temp = abs_s(sub(curQGain, CodCng.IRef));//lsc 如果增益误差过大,也认为当前的filterA不再适用,要发送sid帧
if(temp > ThreshGain) {
*Ftyp = 2;
}
else {
/* no transmission */
*Ftyp = 0;
}
CodCng.IRef保存的是前一帧的curQGain

另外还有一个依据,就是滤波器相似度(对应sid帧中的滤波器参数,其实就是量化后的lsp系数)
这里采用的算法是坂仓距离,笔者观察了代码后,认为,直接从坂仓距离来理解的话,是十分困难的
但从另外一个角度,这段代码就会变得非常容易理解,

笔者的理解如下:
Cod_Cng假设了一个平均波器filterA(由上一帧sid语音包产生的),filterA是否能表征当前语音帧的滤波器,是需要判断的,不相符,就要
重新传输filterA,即所谓的sid帧.
判断的依据,即是将语音信号通过filterA进行逆向滤波,这样可以得出残差信号,
根据残差信号的能量,与CodCng.Ener(实际残差信号能量)进行比较,误差在指范围内,即可认为当前的filterA是有效的.
否则,需要编解码双方更新filterA,即发送sid帧

我们设残差能量为
239
Σ x[n]^2 ---- <式1>
n=0
用lpc产生的滤器量容易得到
10
x[n]= Σ y[n-i] * a(i) 其中a(1)=1 其它的a(i)就是lpc系数
i=0
代入式1得
239 10
Σ ( Σ y[n-i] * a(i) ) ^ 2 (真希望能用纸和笔来直接写这个推导过程啊...@_@)
n=0 i=0

化为
239 10 10
Σ ( Σ y[n-i] * a(i) ) * ( Σ y[n-j] * a(j) )
n=0 i=0 j=0

接下来的做法,想必读者们都猜到了,变量替换,交换求和次序,就会得到代码中所用到的逆滤波算法过程,
239 10 10
Σ Σ Σ ( y[n-i] * a(i) * y[n-j] * a(j) )
n=0 i=0 j=0

10 10 239
Σ Σ a(i)*a(j) Σ ( y[n-i] * y[n-j] ) @_@ 手酸了...还是纸和笔好啊....
i=0 j=0 n=0

10 10
Σ Σ a(i)*a(j) * R(i-j)
i=0 j=0

再转化为

10 10
Σ Σ a(m+j)*a(j) * R(m)
m=0 j=0

10 10
Σ R(m) Σ a(m+j)*a(j)
m=0 j=0

10
Σ R(m) * Ra(m) ---- 终于结束了
m=0

这是大致的推倒过程
我们来观察LpcDiff这个函数,它就是来判断滤波器相似度的
代码片段如下:
Flag LpcDiff(Word16 *RC, Word16 ShRC, Word16 *ptrAcf, Word16 alpha)
{//lsc ptrAcf就是文档中的Rt[j] alpha是能量
Word32 L_temp0, L_temp1;
Word16 temp;
int i;
Flag diff;

L_temp0 = 0L;
for(i=0; i<=LpcOrder; i++) {
temp = shr(ptrAcf[i], 2); /* + 2 margin bits */
L_temp0 = L_mac(L_temp0, RC[i], temp);
}

temp = mult_r(alpha, FracThresh);
L_temp1 = L_add((Word32)temp, (Word32)alpha);//lsc 这里能量里面的的 Ener * 1.2136
temp = add(ShRC, 9); /* 9 = Lpc_justif. * 2 - 15 - 2 */
L_temp1 = L_shl(L_temp1, temp);

if(L_temp0 <= L_temp1) diff = 1; /* G723.1 maintenance April 2006*/
/* Before : if(L_temp0 < L_temp1) diff = 1; */
else diff = 0;
return(diff);
}
RC 其中RC就是 推倒过最后的Ra(m)数组 ptrAcf就是R(m)数组

Cod_Cng分析到这一步,接下去的就迎刃而解了
ComputePastAvFilter 计算出平均滤波器filterA
这个是根据当前保存的4帧自相关系数利用德宾算法计算出lpc系数(平均滤波器)
Ra(m)的生成,自然需要这些

不处于尾响阶段,需要更新相应的平均滤波器lpc系数,用于静音裁决逆滤波时使用,也就是上一节分析的静音检测时,要用于的滤波器
/* If adaptation enabled, fill noise filter */
if ( !VadStat.Aen ) {
for(i=0; i<LpcOrder; i++) VadStat.NLpc[i] = CodCng.SidLpc[i];
}

这里更新Ra(m)
/* Compute autocorr. of past average filter coefficients */ //lsc 根据平均滤波器,算RC,用于逆向滤波的,可以很容易推导,残差信号的能量,可以用这个系数与输入信号的自相关算出 Σx[n]^2 = Σ(Σy[n-i]a(i))^2 利用下标变换推导出
CalcRC(CodCng.SidLpc , CodCng.RC, &CodCng.ShRC);

更新后平均滤波器与当前滤波器进行比较,如果差异过大,则要修正更新
if(LpcDiff(CodCng.RC, CodCng.ShRC, CodCng.Acf, *CodCng.Ener) == 0){//lsc 如果相差太多,取当前帧的滤波器系数
for(i=0; i<LpcOrder; i++) {
CodCng.SidLpc[i] = curCoeff[i];
}
CalcRC(curCoeff, CodCng.RC, &CodCng.ShRC);//lsc 根据当前滤波器算RC
}

转成lsp,并量化,用于后继的振铃减法使用,做一些内存更新,如果保留当前增益估值,用于下次检验滤波器相似度使用
/*
* Compute SID frame codes
*/
/* Compute LspSid */
AtoLsp(CodCng.LspSid, CodCng.SidLpc, CodStat.PrevLsp);
Line->LspId = Lsp_Qnt(CodCng.LspSid, CodStat.PrevLsp);
Lsp_Inq(CodCng.LspSid, CodStat.PrevLsp, Line->LspId, 0);

Line->Sfs[0].Mamp = curQGain;
CodCng.IRef = curQGain;
CodCng.SidGain = Dec_SidGain(CodCng.IRef);

利用增益估值与一个随机数种子,来生成随机激励
/*
* Compute new excitation
*/
if(CodCng.PastFtyp == 1) {
CodCng.CurGain = CodCng.SidGain;
}
else {//lsc 增益为历史值与当前解码的增益按比例混合 7/8 + 1/8
CodCng.CurGain = extract_h(L_add( L_mult(CodCng.CurGain,0x7000),
L_mult(CodCng.SidGain,0x1000) ) ) ;
}
Calc_Exc_Rand(CodCng.CurGain, CodStat.PrevExc, DataExc,
&CodCng.RandSeed, Line);//lsc 采用随机的方式,生成舒适噪声的激励 参数为解码的CurGain,与随机值RandSeed,这两个均在编解码双方同步,这也将导致PrevExc在编解码双方同步

现在来看这个函数:Calc_Exc_Rand
随机激励由两部分构成,一部分是自适应,即来自历史解码激励
另一个部来自随机位置脉冲

首先取一段随机的自适应
/*
* generate LTP codes
*/
Line->Olp[0] = random_number(21, nRandom) + (Word16)123;
Line->Olp[1] = random_number(19, nRandom) + (Word16)123; /* G723.1 maintenance April 2006 */
/* Before : Line->Olp[1] = random_number(21, nRandom) + (Word16)123; */
for(i_subfr=0; i_subfr<SubFrames; i_subfr++) { /* in [1, NbFilt] */
Line->Sfs[i_subfr].AcGn = random_number(NbFilt, nRandom) + (Word16)1;//lsc 自适应码本增益,随机地取一组
}
Line->Sfs[0].AcLg = 1;
Line->Sfs[1].AcLg = 0;
Line->Sfs[2].AcLg = 1;
Line->Sfs[3].AcLg = 3;

可以看出,是从过去的 123 ~ 143这个区间随机抽取连续的240个历史激励
并随机地取自适应激励
由此可看出,编解码双方只要初始随机种子相同，就能保证构造出相同的随机激励

接下来构造随机脉冲,做法是随机地生成13位bit串
其中2bit做为两个子帧脉冲位置(奇或偶)
另外11bit做为两个子帧脉的符号(正或负)
代码片段如下:

/* Signs and Grids */
ptr_TabSign = TabSign;
ptr1 = offset;
for(iblk=0; iblk<SubFrames/2; iblk++) {
temp = random_number((1 << (NbPulsBlk+2)), nRandom);//lsc 取得13个随机的 0,1组合串
*ptr1++ = temp & (Word16)0x0001;//lsc 取一个随机位,对应与网格是否要移到奇数位置
temp = shr(temp, 1);
*ptr1++ = add( (Word16) SubFrLen, (Word16) (temp & 0x0001) );//lsc 再取一个随机位
for(i=0; i<NbPulsBlk; i++) {//lsc 剩下11个随机0,1串是用做符号的
*ptr_TabSign++= shl(sub((temp & (Word16)0x0002), 1), 14);//lsc 要么是1,要么是-1,然后左移14位
temp = shr(temp, 1);
}
}

接下来随机地指定脉冲所在的位置
/* Positions */
ptr_TabPos = TabPos;
for(i_subfr=0; i_subfr<SubFrames; i_subfr++) {

for(i=0; i<(SubFrLen/Sgrid); i++) tmp[i] = (Word16)i;
temp = (SubFrLen/Sgrid);
for(i=0; i<Nb_puls[i_subfr]; i++) {
j = random_number(temp, nRandom);//lsc 0-30之间的随机数,脉冲位置是随机的,分别会有5~6个随机脉冲
*ptr_TabPos++ = add(shl(tmp[(int)j],1), offset[i_subfr]);//lsc 位置乘2加上随机偏移1个位置
temp = sub(temp, 1);//lsc 随机数范围减1
tmp[(int)j] = tmp[(int)temp];//lsc 这里将被选中的位置排除,避免下一次被随中
}
}

接下来,就是确认随机脉冲的增益了,
依据就是itu A.4节的公式 A-19,对增益gf的要求是,构造出来的随机激励与CodCng.CurGain(舒适噪音的增益估值)
的平方差因尽量小(或者为0)
根据二次方程的求根公式，如果delta（就是二次方程的判别式）大于0,则求根
如果小于0，这时就直接计算该二次多项式的最小值

这些工作完成之后,Cod_Cng的工作也就基本结束了

这章分析完,g723的代码分析只剩后置滤波与丢包补偿需要整理了。

断断续续，整理分析9723代码居然长达半年多，
这半年事情实在太多了，为宝贝儿子的出生，我日夜地祈祷
笔者在完成这些后，希望能在2012.12.21之前写一个g723的编解码。

林绍川
2012.01.13于杭州