g723源码详细分析-11-多脉码激励编码

上一节分析了g723低速率下的固定码本搜索,
在高速率下面,对去了自适应激励成份后的语音信号,
g723采用的是多脉冲编码方式.本文现做一个简要的闸述

在分析代码之前,仍然是要做一些公式推导,否则不易读懂代码
过程其实与固定码本搜索时的方式类似.

多脉冲激励的定义:
激励是由有限个经过最优估值的脉冲构成(每个脉冲的位置与增益都是
经过最优估值的)
那么很容易得出相应的误差公式,假设脉冲的个数是m
                    M
e[n]= t[n] -  Σ  g(k) * h[n-n(k)]
                  k=1

n(k) 代表第k个脉冲的位置 g(k)代表第k个脉冲的增益
误差均方:
       N                 N            M
E = Σ  e[n]^2 =  Σ  (t[n] -  Σ  g(k) * h[n-n(k)])^2 ---- <式1>
     n=1             n=1         k=1

对
 M
 Σ  g(k) * h[n-n(k)]
k=1
这个式子进行一些处理,将其变形为
 N
 Σ  g(l) * h[n-l]
l=1
当l=n(k) g(l)为对应位置的脉冲值g(k),
其它的l取值时g(l)为零

<式1>此时可以变形为
    N            N           N
E = Σ  e[n]^2 =  Σ  (t[n] -  Σ  g[l] * h[n-l])^2
   n=1          n=1         l=1

       N                   N         N                         N      N      N
E = Σ t[n]^2 - 2 * Σ t[n] * Σ (g[l] * h[n-l]) +  Σ     Σ       Σ  g[j]*h[n-j]*g[l]*h[n-l]
     n=1              n=1       l=1                      n=1  j=1    l=1

    N             N        N                    N   N                N
  = Σ t[n]^2 -2 * Σ g[l] * Σ (t[n] * h[n-l]) +  Σ   Σ  g[j] * g[l] * Σ (h[n-j]*h[n-l])   --- 交换求和次序
   n=1           l=1      n=1                  l=1 j=1              n=1

    N                    M          N                                  M          M           N
  = Σ t[n]^2 - 2 * Σ g[k] * Σ (t[n] * h[n-n(k)]) +  Σ g[k] *  Σ  g[j] *  Σ (h[n-n(j)]*h[n-n(k)])

   n=1                k=1      n=1                              k=1       j=1        n=1

--- 变量替换,去除g(x)等于0的项

E分别对每个g(k)分别进行求导,g(k)最优估值必定要导数为零处,将会得到M个等式
 N                                M          N
 Σ (t[n] * h[n-n(k)]) =  Σ  g[j] * Σ (h[n-n(j)]*h[n-n(k)])             k = 1,2, ..., M
n=1                            j=1       n=1
可以化简为
                     M
 Rth(n(k)) = Σ g[j] Rhh(n(j),n(k)) ---<式2>         Rth是t与h的互相关m, Rhh是h的自相关
                     j=1
            
而<式1>可以化简为:
       N            M
E = Σ t[n]^2 - Σ g[k]*Rt(n(k)) ---<式3>
      n=1        k=1
   
由于n(k)是未知的,要求解这个方程是极其复杂的,
在难以求得最优解的情况下,采用一种简化算法,求得次优解

即每次只求一个最优解的位置,<式2>就变成这样:
Rth(n(1)) = g[1] Rhh(n(1),n(1)) ---<式4>
        Rth(n(1))
g[1] = --------------
       Rhh(n(1),n(1))
从<式3>可以看出,n(1)是满足下式最大值的位置
   Rth(n(1))^2
  --------------   ---<式5>
  Rhh(n(1),n(1))

这样可以很容易地求出n(1),并根据<式4>求出g[1]

然后从 t中扣除g[1]与h卷积后的成份,更新t
t[n] = t[n] - g[1] * h[n-n(1)] ---<式6>

用更新后的t继续进行下一轮的搜索,直到找出M个脉冲位置以及它们的增益.

以上是公式推导过程,与itu的文档中的描述相当接近了,不过itu采取的是更为化简的估值算法.
itu认为每个脉冲的增益都是相同的,用
    max(|Rth(i)|)
G= ------------   ---- 等同于itu文档中的公式25
    Rhh(0,0)
来对增益进行估值,
在G-3.2db ~ G + 6.4db这个范围内,将G分成4个等级值,进行搜索脉冲的位置
脉冲位置搜索的依据是<式5>,并且忽略了分母的差异,只搜索分子最大的
每搜索一次按照<式6>更新输入的残差信号(代码里直接更新自相关,道理是相同的)
最后,得到估值增益,以及脉冲位置

推倒完毕(^_^),现在进入代码:

Find_Best
这个函数描述了搜索脉冲位置与增益的过程

首先是构造冲激响应,在后面会提到,搜索会根据基音周期是否小于子帧长,再调整冲激响应
再进行一次搜索,然后与未调整冲激响应的搜索结果进行比较,取更接近的那一组脉冲位置与增益

    /* Update Impulse response */
    if ( Olp < (Word16) (SubFrLen-2) ) {
        Temp.UseTrn = (Word16) 1 ;
        Gen_Trn( Imr, ImpResp, Olp ) ;//lsc 相当于与olp为周期间隔脉进行卷积
    }
    else {
        Temp.UseTrn = (Word16) 0 ;
        for ( i = 0 ; i < SubFrLen ; i ++ )
            Imr[i] = ImpResp[i] ;
    }

计算冲激响应的自相关,以相冲激响应与目标向量的互相关,并对其进行归一化处理,这里不贴代码了
计算的结果存放在ImrCorr(冲激响应自相关) ErrBlk(冲激响应与目标向量的互相关)
ImrCorr是扩大2^(9 + Exp)倍的
ErrBlk是扩大2^(10 + Exp)倍的(这样便于后继的计算)

由于脉冲要么全位于奇位置,要么全位于偶位置,所以搜索分成奇偶两次
在搜索循环中,首先进行增益估值
增益估值使用的是FcbkGainTable这个码本表里的值,该码本表是与固定码本增益共享的
采用乘法运算,绕开费时的除法运算,即将增益码本FcbkGainTable里的值与冲激响应的零时自相关(即Rhh(0))相乘
然后与ErrBlk数组里绝对值最大的相减(也就是max(|Rth(i)|)),取差值最小的那个增益,然后在它附近的4个增益值中
进行搜索最优增益估值
代码如下:

    /* Search for the best sequence */ //lsc 这里是指奇偶两次循环
    for ( k = 0 ; k < Sgrid ; k ++ ) {

        Temp.GridId = (Word16) k ;

        /* Find maximum amplitude */ //lsc 找了最大的那个 |d[j]|
        Acc1 = (Word32) 0 ;
        for ( i = k ; i < SubFrLen ; i += Sgrid ) {
            Acc0 = L_abs( ErrBlk[i] ) ;
            if ( Acc0 >= Acc1 ) {
                Acc1 = Acc0 ;
                Temp.Ploc[0] = (Word16) i ;//lsc 互相关绝对值最大的那个的下标,这也是搜索过程中被搜索出的第一个脉冲的位置
            }
        }

        /* Quantize the maximum amplitude */
        Acc2 = Acc1 ;//lsc Acc2就是最大的 |d[j]|
        Acc1 = (Word32) 0x40000000L ;
        MaxAmpId = (Word16) (NumOfGainLev - MlqSteps) ;

        //lsc 这里再次利用乘法绕开除法,将 G* Rhh[0]得到的结果与 Acc2相减,差值最小的那个就是G
        for ( i = MaxAmpId ; i >= MlqSteps ; i -- ) {
            Acc0 = L_mult( FcbkGainTable[i], ImrCorr[0] ) ;
            Acc0 = L_sub( Acc0, Acc2 ) ;//lsc 同一数量级,所以可减
            Acc0 = L_abs( Acc0 ) ;
            if ( Acc0 < Acc1 ) {
                Acc1 = Acc0 ;
                MaxAmpId = (Word16) i ;
            }
        }
        MaxAmpId -- ;//lsc 得到的值在降一级搜索,见itu文档 MaxAmpId最小值是1

        //lsc 位置的搜索过程，我们看到一些简化,即G都是一样的,每个位置的搜索并不对增益做任何更新(照理应计算出G，并用其来更新互相关)
        //按G的值循环搜索
        for ( i = 1 ; i <=2*MlqSteps ; i ++ ) {

            for ( j = k ; j < SubFrLen ; j += Sgrid ) {
                WrkBlk[j] = ErrBlk[j] ;//lsc 取出偶数位/奇数位的 互相关
                OccPos[j] = (Word16) 0 ;
            }
            Temp.MampId = MaxAmpId - (Word16) MlqSteps + (Word16) i ;//lsc 这样做等最从0位置开始搜索

            MaxAmp = FcbkGainTable[Temp.MampId] ;

            if ( WrkBlk[Temp.Ploc[0]] >= (Word32) 0 )
                Temp.Pamp[0] = MaxAmp ;
            else
                Temp.Pamp[0] = negate(MaxAmp) ;

            OccPos[Temp.Ploc[0]] = (Word16) 1 ;//lsc 选中第一个位置,对均方差求导很容易得到这个结论

            for ( j = 1 ; j < Np ; j ++ ) {

                Acc1 = (Word32) 0xc0000000L ;

                for ( l = k ; l < SubFrLen ; l += Sgrid ) {

                    if ( OccPos[l] != (Word16) 0 )//lsc 表示这个位置已经被选中了
                        continue ;

                    Acc0 = WrkBlk[l] ;//lsc 这个是互相关
                    Acc0 = L_msu( Acc0, Temp.Pamp[j-1],
                            ImrCorr[abs_s((Word16)(l-Temp.Ploc[j-1]))] ) ;//lsc 这里是在更新互相关,因为分母都一样,所以只要分子最大,就可以选定这个位置
                    WrkBlk[l] = Acc0 ;
                    Acc0 = L_abs( Acc0 ) ;
                    if ( Acc0 > Acc1 ) {
                        Acc1 = Acc0 ;   
                        Temp.Ploc[j] = (Word16) l ;//lsc Ploc的意思是 pluse location,即脉冲位置,
                    }
                }

                //lsc 这里其实只是记录符号
                if ( WrkBlk[Temp.Ploc[j]] >= (Word32) 0 )
                    Temp.Pamp[j] = MaxAmp ;
                else
                    Temp.Pamp[j] = negate(MaxAmp) ;

                //lsc 选定位置
                OccPos[Temp.Ploc[j]] = (Word16) 1 ;
            }

            /* Compute error vector */
            for ( j = 0 ; j < SubFrLen ; j ++ )
                OccPos[j] = (Word16) 0 ;

            for ( j = 0 ; j < Np ; j ++ )
                OccPos[Temp.Ploc[j]] = Temp.Pamp[j] ;

            //lsc 计算卷积
            for ( l = SubFrLen-1 ; l >= 0 ; l -- ) {
                Acc0 = (Word32) 0 ;
                for ( j = 0 ; j <= l ; j ++ )
                    Acc0 = L_mac( Acc0, OccPos[j], Imr[l-j] ) ;
                Acc0 = L_shl( Acc0, (Word16) 2 ) ;
                OccPos[l] = extract_h( Acc0 ) ;
            }

            /* Evaluate error */
            //lsc 计算均方差
            Acc1 = (Word32) 0 ;
            for ( j = 0 ; j < SubFrLen ; j ++ ) {
                Acc1 = L_mac( Acc1, Tv[j], OccPos[j] ) ;
                Acc0 = L_mult( OccPos[j], OccPos[j] ) ;
                Acc1 = L_sub( Acc1, L_shr( Acc0, (Word16) 1 ) ) ;
            }

            //lsc 如果大，说明均方差小 (a^2 - 2ab + b^2) 很容易得到这个结论
            //则选这组位置
            if ( Acc1 > (*Best).MaxErr ) {
                (*Best).MaxErr = Acc1 ;
                (*Best).GridId = Temp.GridId ;
                (*Best).MampId = Temp.MampId ;
                (*Best).UseTrn = Temp.UseTrn ;
                for ( j = 0 ; j < Np ; j ++ ) {
                    (*Best).Pamp[j] = Temp.Pamp[j] ;
                    (*Best).Ploc[j] = Temp.Ploc[j] ;
                }
            }
        }
    }

致此,g723编码部分的代码分析绝大部分已经完成了,剩余的一些代码是更新自适应码本,打包字节等
这将在下一个章节做一个简述

这里做一个g723编码模块的简要小结

step1:LPC分析
    这一步骤就是要得出一个全极点的声道系统函数,并且依此生成声道的冲激响应
    它包括lpc参数计算(莱文森德宾递推公式) 转lsp,并量化反量化等
    
step2:
    语音数据处理,包括共振峰感知加权,基音周期估算(采用自相关方法)

step3:
    自适应激励码本搜索

step4:
    去除自适应激励贡献后的语音数据的激励码本搜索
    它包含两种速率下面的不同处理,分别为高速率下多脉冲激励码搜索,以及低速率下的固定码本搜索
    
    
                                                               林绍川
                                                               2011.10.06 于杭州