speex源码分析-3-自适应激励

本节来分析speex的自适应激励,
与g723 729一样,自适应激励最主要的课题就是基音周期的搜索,
这个自然是通过自相关算法来解决的.


而自适应码本则是由固定码本不断迭代计算出来的


首先来看一下,开环基音增益ol_gain.
直观地解释这个变量的含义,就是搜索出自适应激励码本后,要根据
这个增益做一个比例缩放.


如我们在之前章节中提到的ol_gain是由原始激励平均能量的开方,
它的计算见:http://blog.csdn.net/lsccsl/article/details/8126338


speex会先对ol_gain做一个量化,(便于网络报文传输能使用更少的bit)
   {
      int qe = (int)(floor(.5+3.5*log(ol_gain*1.0/SIG_SCALING)));
      if (qe<0)
         qe=0;
      if (qe>31)
         qe=31;
      ol_gain = exp(qe/3.5)*SIG_SCALING;
      speex_bits_pack(bits, qe, 5);
   }
从算法上看,就是取对数,减少量化的bit位,报文中实际传输的是qe的值


之后的代码是针对每一子帧进行的
      /* LSP interpolation (quantized and unquantized) */
      lsp_interpolate(st->old_lsp, lsp, interp_lsp, st->lpcSize, sub, st->nbSubframes);//lsc lsp插值
      lsp_interpolate(st->old_qlsp, qlsp, interp_qlsp, st->lpcSize, sub, st->nbSubframes);


      /* Make sure the filters are stable */
      lsp_enforce_margin(interp_lsp, st->lpcSize, LSP_MARGIN);//lsc lsp稳定性
      lsp_enforce_margin(interp_qlsp, st->lpcSize, LSP_MARGIN);
这个是lsp插值和稳定性调整,注意每个子帧的搜索系数会调整


lsp转换成lpc
      /* Compute interpolated LPCs (quantized and unquantized) */
      lsp_to_lpc(interp_lsp, interp_lpc, st->lpcSize,stack);//lsc lsp转lpc


      lsp_to_lpc(interp_qlsp, interp_qlpc, st->lpcSize, stack);
interp_lpc即为当前子帧插值lpc系数  
interp_qlpc为当前子帧量化后的插值lpc系数


计算感知加权滤波器
      /* Compute bandwidth-expanded (unquantized) LPCs for perceptual weighting */
      bw_lpc(st->gamma1, interp_lpc, bw_lpc1, st->lpcSize);
      if (st->gamma2>=0)
         bw_lpc(st->gamma2, interp_lpc, bw_lpc2, st->lpcSize);
      else
      {
         for (i=0;i<st->lpcSize;i++)
            bw_lpc2[i]=0;
      }
      
计算lpc量化后,对量的激励real_exc,


      /*FIXME: This will break if we change the window size */
      speex_assert(st->windowSize-st->frameSize == st->subframeSize);
      if (sub==0)
      {
         for (i=0;i<st->subframeSize;i++)
            real_exc[i] = sw[i] = st->winBuf[i];//lsc 这里又把sw的值还原了,还原成原始语音信号
      } else {
         for (i=0;i<st->subframeSize;i++)//lsc 这里又把sw的值还原了,还原成原始语音信号
            real_exc[i] = sw[i] = in[i+((sub-1)*st->subframeSize)];
      }
      fir_mem16(real_exc, interp_qlpc, real_exc, st->subframeSize, st->lpcSize, st->mem_exc2, stack);//lsc real_exc是原始语音语号通过反量化的lpc系数弄成的激励.
      
real_exc将被做为解码激励的一个参照用于计算固定码本激励时的增益,这个将在之后提到,现在只要记住它是怎么得来的即可


计算冲激响应
      if (st->complexity==0)
         response_bound >>= 1;//lsc 这话句的意思是乘2,因为如有感加权,冲激响应将会有20点,是没有感加权的2倍
      compute_impulse_response(interp_qlpc, bw_lpc1, bw_lpc2, syn_resp, response_bound, st->lpcSize, stack);//lsc 计算出冲激响应 syn_resp号


计算零输入响应
      iir_mem16(ringing, interp_qlpc, ringing, st->subframeSize, st->lpcSize, mem, stack);
      for (i=0;i<st->lpcSize;i++)
         mem[i]=SHL32(st->mem_sw[i],1);
      filter_mem16(ringing, bw_lpc1, bw_lpc2, ringing, st->subframeSize, st->lpcSize, mem, stack);//lsc 计算零输入响应,要用到之前存储的语音信号
      
这里简单地分析一下iir_mem16函数:
void iir_mem16(const spx_word16_t *x, const spx_coef_t *den, spx_word16_t *y, int N, int ord, spx_mem_t *mem, char *stack)
{//lsc den分母 x:输入 y:输出 N:帧长 ord:阶 mem:内存更新,从代码上看,主要用于下帧的计算(零输入响应等)
  int i,j;
    spx_word16_t yi,nyi;
//lsc mem[10] m[9]...m[0]  m[10]=-az(10)*y[n-10]  m[9]=m[10]-az[9]y[n-9] ... m[0]=m[1]-az[1]y[n-1]
//lsc 每轮计算时,mem[0]即完成了iir部分的运算,加上x[i],就完成了fir部分的运算,而每轮被mem被更新
//lsc 同时mem也会被函数之外的滤波器所使用,作为初始输入,如计算零输入响应
//lsc 这一点显得与g723不同,g723保留计算结果,而az系数用的是下一帧的,speex这种处理方法,计算下帧时,前几个输入仍用当前帧的,之后才完全使下一帧的az系数
  for (i=0;i<N;i++)
    {
      yi = EXTRACT16(SATURATE(ADD32(EXTEND32(x[i]),PSHR32(mem[0],LPC_SHIFT)),32767));
      nyi = NEG16(yi);
      for (j=0;j<ord-1;j++)
      {
        mem[j] = MAC16_16(mem[j+1],den[j],nyi);
      }
      mem[ord-1] = MULT16_16(den[ord-1],nyi);
      y[i] = yi;
    }
}
filter_mem16 函数与iir_mem16类似,只不过添加了分子部分的处理,即fir部分不大一样.


st->mem_sp st->mem_sw 将在每一子帧结束时更新,在后面会提到


计算目标信号,即原始语音信号感知加权后,扣减零输入响应,
      /* Compute weighted signal */
      for (i=0;i<st->lpcSize;i++)
         mem[i]=st->mem_sw[i];
      filter_mem16(sw, bw_lpc1, bw_lpc2, sw, st->subframeSize, st->lpcSize, mem, stack);//lsc 对原始语音信号计算感知加权,存入sw
      
      if (st->complexity==0)
         for (i=0;i<st->lpcSize;i++)
            st->mem_sw[i]=mem[i];
      
      /* Compute target signal (saturation prevents overflows on clipped input speech) */
      for (i=0;i<st->subframeSize;i++)
         target[i]=EXTRACT16(SATURATE(SUB32(sw[i],PSHR32(ringing[i],1)),32767));//lsc 扣减零输入响应,得到目标信号 targe(即扣减零输入响应的加权语音信号)




然后就是自适应码本搜索 pitch_search_3tap


从函数名字,大约可以猜出,这个是在基音延后附近搜索,每三个相邻点加权形成一个自适应激励(g723是五个点)
算法跟g723极其类似,除了加权采样点的个数不同之外
这个函数的参数较多,笔者就主要参数做一个简要说明
target:搜索的目标信号
sw:原始语音的感知加权信号
ak:声道系数滤波系数 awk1/awk2:感知加权滤波系数
exc:保存解码自适应激励
start:基音周期搜索的最小延后 end:最大基音搜索延后
p:阶
nsf:子帧长度
exc2:自适应激励码本(历史解码的自适应激励)
r:单位冲激响应(ak awk1 awk2)
par:自适应码本的搜索配置参数(自适应激励码本的加权系数码本表等)


先对当前子帧做一个开环基音搜索(g723的基音周期是针对当前帧,而speex则针对某个子帧单独计算基音周期):
open_loop_nbest_pitch采用的算法,就是自相关算法,这点与g723相同
有兴趣的读者可以参考:
http://blog.csdn.net/lsccsl/article/details/6425568
代码复杂,但表达的算法却是相对简单


pitch_gain_search_3tap
得到当前子帧的基音周期pitch后(会选出若干个备选的基音周期),
在每个备选基音延后 pitch-1 pitch pitch+1 计算这些延后的历史解码激励与当前综合滤波器(ak awk1 awk2)冲激响应的卷积


摘出部分代码片段:
      for (j=0;j<nsf;j++)
      {
         if (j-pp<0)
            e[j]=exc2[j-pp];//lsc exec2保存的是历史激励 正常情况
         else if (j-pp-pitch<0)
            e[j]=exc2[j-pp-pitch];//lsc 过头了,循环 j会取值到 0-39 而 pitch只有17,就需要循环取值
         else
            e[j]=0;//lsc 过了一个周期还不够,speex直接置0?
      }
以上为拷贝出 pitch-1 pitch pitch+1 延后的历史解码激励


计算它们的冲激响应:
      for (j=0;j<p;j++)
         mm[j] = 0;
      iir_mem16(e, ak, e, nsf, p, mm, stack);
      for (j=0;j<p;j++)
         mm[j] = 0;
      filter_mem16(e, awk1, awk2, e, nsf, p, mm, stack);
      for (j=0;j<nsf;j++)
         x[2][j] = e[j];//lsc 得到 -(pitch+1)处的历史激励冲激响应 a(z) * aw1(z)/aw2(z) 系统得到的输出


   for (i=1;i>=0;i--)
   {
      spx_word16_t e0=exc2[-pitch-1+i];
#ifdef FIXED_POINT
 ???
      /* Scale excitation down if needed (avoiding overflow) */
      if (scaledown)
         e0 = SHR16(e0,1);
#endif
      x[i][0]=MULT16_16_Q14(r[0], e0);
      for (j=0;j<nsf-1;j++)
         x[i][j+1]=ADD32(x[i+1][j],MULT16_16_P14(r[j+1], e0));//lsc 分别得到了 -pitch , -pitch + 1处的历史激励冲激响应输出
   }
这里用到了和g723相似的小技巧,减少了运算量


对加权增益码本的搜索是根据误差能量最小进行的,
假设三个基音延后输出矩阵分别为 L1 L2 L3
(a*L1 + b*L1 + c*L3 - t)(a*L1 + b*L1 + c*L3 - t) 展开
(a*L1)^2 + (b*L1)^2 + (c*L3)^2  + 2ab(L1 * L2) + 2ac(L1 * L2) + 2bc(L2 * L3) -2a(L1 *t) -2b(L2 * t) - 2c(L3 * t) + t^2
其中对了固定的子帧,t^2是常数项,只要关心前九项即可


增益码本的数据定义 gain_cdbk_nb
笔者摘选部码本中的部分码字,观察它们的规律
const signed char gain_cdbk_nb[512] = {//lsc 这里面是系数,类似于g723,每组4个系数之间是有联系的,简单的理解,就是将3个相邻的历史激励加权,与搜索向量比较,能量差最小的将被选中
-32, -32, -32, 0,//lsc 0 = 0/2
-28, -67, -5, 33,//lsc -28+32=4 -67+32=-35 -5+32=27  (4+35+27)/2=33
-42, -6, -32, 18,
-57, -10, -54, 35,
...//其它码字
}
这个码本每4个数字为一组,每组的前三个数字为增益加权系数,最后一个值是增益加权系数绝对值和除2
注意这个规律,每组最后的那个数字,比第一组的 0 第二组的 33,是用来在搜索中控制增益加权系数总和
g723也有对应的算法与speex这个处理相对应,但speex显然做了简化,看起来也更明了
(不理解的读者可以先不管最后的那个数字,之后的分析会让读者自然明白它的含义)


接下来的计算比较机械
   for (i=0;i<3;i++)
      corr[i]=inner_prod(x[i],new_target,nsf);//lsc 分别计算三个延后与搜索向量的内积
   for (i=0;i<3;i++)
      for (j=0;j<=i;j++)
         A[i][j]=A[j][i]=inner_prod(x[i],x[j],nsf);//lsc 三个延后之间的内积
以上计算各个矩阵相乘的结果


      C[0]=corr[2];//lsc 对应 Lx * t
      C[1]=corr[1];//lsc 对应 Lx * t
      C[2]=corr[0];//lsc 对应 Lx * t
      C[3]=A[1][2];//lsc Lx * Lx
      C[4]=A[0][1];//lsc Lx * Lx
      C[5]=A[0][2];//lsc Lx * Lx      
      C[6]=A[2][2];//lsc Lx * Lx
      C[7]=A[1][1];//lsc Lx * Lx
      C[8]=A[0][0];//lsc Lx * Lx
以上将它们保存到C16这个数组


      C[6]*=.5*(1+.02*plc_tuning);
      C[7]*=.5*(1+.02*plc_tuning);
      C[8]*=.5*(1+.02*plc_tuning);//lsc 除2,因为其它六个是有系数2的
以上处理系数为2的项


pitch_gain_search_3tap_vq:
这个函数开始在加权增益码本里搜索
  for (i=0;i<gain_cdbk_size;i++) {
         
    ptr = gain_cdbk+4*i;
    g[0]=ADD16((spx_word16_t)ptr[0],32);
    g[1]=ADD16((spx_word16_t)ptr[1],32);
    g[2]=ADD16((spx_word16_t)ptr[2],32);
    gain_sum = (spx_word16_t)ptr[3];
         
    sum = compute_pitch_error(C16, g, pitch_control);
         
    if (sum>best_sum && gain_sum<=max_gain) {//lsc 值越大,说明是更好的码本 gain_sum<=max_gain 三个采样点的放大倍数绝对值之和不能大于某个阀值
      best_sum=sum;
      best_cdbk=i;
    }


compute_pitch_error:计算 -[(a*L1)^2 + (b*L1)^2 + (c*L3)^2  + 2ab(L1 * L2) + 2ac(L1 * L2) + 2bc(L2 * L3) -2a(L1 *t) - 2b(L2 * t) - 2c(L3 * t)]
static inline spx_word32_t compute_pitch_error(spx_word16_t *C, spx_word16_t *g, spx_word16_t pitch_control)
{
   spx_word32_t sum = 0;
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[0],pitch_control),C[0]));
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[1],pitch_control),C[1]));
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[2],pitch_control),C[2]));
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[1]),C[3]));
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[1]),C[4]));
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[0]),C[5]));
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[0]),C[6]));
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[1],g[1]),C[7]));
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[2]),C[8]));
   return sum;
}


在这里可以看到gain_sum max_gain的作用, gain_sum就是码字里的第4个数
max_gain 是通过上一子帧计算出来的:
   *cumul_gain = 0.03125*MAX32(1024,*cumul_gain)*params->gain_cdbk[4*best_gain_index+3];
params->gain_cdbk[4*best_gain_index+3]:就是码字中的第4个数


对应每个备选的基音延后,都找出最佳的增益码字,以及它们的误差,选误差最小的那个基音延后,
根据加权码字与基音延后,以及自适应激励码本,计算出解码的自适应激励.
然后更搜索向量,扣除掉解码自适应激励的贡献成分,转给固定码本搜索.


解码当前延迟的激励:
   for (i=0;i<3;i++) 
   {
      int j;
      int tmp1, tmp3;
      int pp=pitch+1-i;
      tmp1=nsf;
      if (tmp1>pp)
         tmp1=pp;
      for (j=0;j<tmp1;j++)
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp]);
      tmp3=nsf;
      if (tmp3>pp+pitch)
         tmp3=pp+pitch;
      for (j=tmp1;j<tmp3;j++)
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp-pitch]);
   }


更新搜索向量
   for (i=0;i<nsf;i++)
   {
      spx_word32_t tmp = ADD32(ADD32(MULT16_16(gain[0],x[2][i]),MULT16_16(gain[1],x[1][i])),
                            MULT16_16(gain[2],x[0][i]));
      new_target[i] = SUB16(new_target[i], EXTRACT16(PSHR32(tmp,6)));
   }


计算误差信号的能量,即固定码本搜索对象的能量,作为函数外比较的依据
   err = inner_prod(new_target, new_target, nsf);
   
至此,完成了自适应激励的搜索,
总结一下speex窄带自适应搜索算法:
1 每个子帧计算基音周期(会算出若干个备选的)
2 每个备选基音周期延后,每相邻三个样点,进行加权增益(增益从一个码本表里搜索出),拟合成一个自适应激励.
3 解码自适应激励,更新搜索向量,准备进行固定码本搜索.






                                                     林绍川 2012-11-01 于杭州