lpc.c

g729b_pc.rar是G729语音压缩标准的附件B的C语言实现标准源代码。

    old_A[i] = A[i] = round(L_shl(t0, 1));
  }
  old_rc[0] = rc[0];
  old_rc[1] = rc[1];

  return;
}




/*-------------------------------------------------------------*
 *  procedure Az_lsp:                                          *
 *            ~~~~~~                                           *
 *   Compute the LSPs from  the LPC coefficients  (order=10)   *
 *-------------------------------------------------------------*/

/* local function */

static Word16 Chebps_11(Word16 x, Word16 f[], Word16 n);
static Word16 Chebps_10(Word16 x, Word16 f[], Word16 n);

void Az_lsp(
  Word16 a[],        /* (i) Q12 : predictor coefficients              */
  Word16 lsp[],      /* (o) Q15 : line spectral pairs                 */
  Word16 old_lsp[]   /* (i)     : old lsp[] (in case not found 10 roots) */
)
{
 Word16 i, j, nf, ip;
 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
 Word16 x, y, sign, exp;
 Word16 *coef;
 Word16 f1[M/2+1], f2[M/2+1];
 Word32 t0, L_temp;
 Flag   ovf_coef;
 Word16 (*pChebps)(Word16 x, Word16 f[], Word16 n);

/*-------------------------------------------------------------*
 *  find the sum and diff. pol. F1(z) and F2(z)                *
 *    F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1)  *
 *                                                             *
 * f1[0] = 1.0;                                                *
 * f2[0] = 1.0;                                                *
 *                                                             *
 * for (i = 0; i< NC; i++)                                     *
 * {                                                           *
 *   f1[i+1] = a[i+1] + a[M-i] - f1[i] ;                       *
 *   f2[i+1] = a[i+1] - a[M-i] + f2[i] ;                       *
 * }                                                           *
 *-------------------------------------------------------------*/

 ovf_coef = 0;
 pChebps = Chebps_11;

 f1[0] = 2048;          /* f1[0] = 1.0 is in Q11 */
 f2[0] = 2048;          /* f2[0] = 1.0 is in Q11 */

 for (i = 0; i< NC; i++)
 {
   Overflow = 0;
   t0 = L_mult(a[i+1], 16384);          /* x = (a[i+1] + a[M-i]) >> 1        */
   t0 = L_mac(t0, a[M-i], 16384);       /*    -> From Q12 to Q11             */
   x  = extract_h(t0);
   if ( Overflow ) {
     ovf_coef = 1;      }

   Overflow = 0;
   f1[i+1] = sub(x, f1[i]);    /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
   if ( Overflow ) {
     ovf_coef = 1;      }

   Overflow = 0;
   t0 = L_mult(a[i+1], 16384);          /* x = (a[i+1] - a[M-i]) >> 1        */
   t0 = L_msu(t0, a[M-i], 16384);       /*    -> From Q12 to Q11             */
   x  = extract_h(t0);
   if ( Overflow ) {
     ovf_coef = 1;      }

   Overflow = 0;
   f2[i+1] = add(x, f2[i]);    /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
   if ( Overflow ) {
     ovf_coef = 1;      }
 }

 if ( ovf_coef ) {
   /*printf("===== OVF ovf_coef =====\n");*/

   pChebps = Chebps_10;

   f1[0] = 1024;          /* f1[0] = 1.0 is in Q10 */
   f2[0] = 1024;          /* f2[0] = 1.0 is in Q10 */

   for (i = 0; i< NC; i++)
   {
     t0 = L_mult(a[i+1], 8192);          /* x = (a[i+1] + a[M-i]) >> 1        */
     t0 = L_mac(t0, a[M-i], 8192);       /*    -> From Q11 to Q10             */
     x  = extract_h(t0);
     f1[i+1] = sub(x, f1[i]);            /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */

     t0 = L_mult(a[i+1], 8192);          /* x = (a[i+1] - a[M-i]) >> 1        */
     t0 = L_msu(t0, a[M-i], 8192);       /*    -> From Q11 to Q10             */
     x  = extract_h(t0);
     f2[i+1] = add(x, f2[i]);            /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
   }
 }

/*-------------------------------------------------------------*
 * find the LSPs using the Chebichev pol. evaluation           *
 *-------------------------------------------------------------*/

 nf=0;          /* number of found frequencies */
 ip=0;          /* indicator for f1 or f2      */

 coef = f1;

 xlow = grid[0];
 ylow = (*pChebps)(xlow, coef, NC);

 j = 0;
 while ( (nf < M) && (j < GRID_POINTS) )
 {
   j =add(j,1);
   xhigh = xlow;
   yhigh = ylow;
   xlow  = grid[j];
   ylow  = (*pChebps)(xlow,coef,NC);

   L_temp = L_mult(ylow ,yhigh);
   if ( L_temp <= (Word32)0)
   {

     /* divide 2 times the interval */

     for (i = 0; i < 2; i++)
     {
       xmid = add( shr(xlow, 1) , shr(xhigh, 1)); /* xmid = (xlow + xhigh)/2 */

       ymid = (*pChebps)(xmid,coef,NC);

       L_temp = L_mult(ylow,ymid);
       if ( L_temp <= (Word32)0)
       {
         yhigh = ymid;
         xhigh = xmid;
       }
       else
       {
         ylow = ymid;
         xlow = xmid;
       }
     }

    /*-------------------------------------------------------------*
     * Linear interpolation                                        *
     *    xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);            *
     *-------------------------------------------------------------*/

     x   = sub(xhigh, xlow);
     y   = sub(yhigh, ylow);

     if(y == 0)
     {
       xint = xlow;
     }
     else
     {
       sign= y;
       y   = abs_s(y);
       exp = norm_s(y);
       y   = shl(y, exp);
       y   = div_s( (Word16)16383, y);
       t0  = L_mult(x, y);
       t0  = L_shr(t0, sub(20, exp) );
       y   = extract_l(t0);            /* y= (xhigh-xlow)/(yhigh-ylow) in Q11 */

       if(sign < 0) y = negate(y);

       t0   = L_mult(ylow, y);                  /* result in Q26 */
       t0   = L_shr(t0, 11);                    /* result in Q15 */
       xint = sub(xlow, extract_l(t0));         /* xint = xlow - ylow*y */
     }

     lsp[nf] = xint;
     xlow    = xint;
     nf =add(nf,1);

     if(ip == 0)
     {
       ip = 1;
       coef = f2;
     }
     else
     {
       ip = 0;
       coef = f1;
     }
     ylow = (*pChebps)(xlow,coef,NC);

   }
 }

 /* Check if M roots found */

 if( sub(nf, M) < 0)
 {
    for(i=0; i<M; i++)
    {
      lsp[i] = old_lsp[i];
    }

 /* printf("\n !!Not 10 roots found in Az_lsp()!!!\n"); */
 }

 return;
}

/*--------------------------------------------------------------*
 * function  Chebps_11, Chebps_10:                              *
 *           ~~~~~~~~~~~~~~~~~~~~                               *
 *    Evaluates the Chebichev polynomial series                 *
 *--------------------------------------------------------------*
 *                                                              *
 *  The polynomial order is                                     *
 *     n = M/2   (M is the prediction order)                    *
 *  The polynomial is given by                                  *
 *    C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 *
 * Arguments:                                                   *
 *  x:     input value of evaluation; x = cos(frequency) in Q15 *
 *  f[]:   coefficients of the pol.                             *
 *                         in Q11(Chebps_11), in Q10(Chebps_10) *
 *  n:     order of the pol.                                    *
 *                                                              *
 * The value of C(x) is returned. (Saturated to +-1.99 in Q14)  *
 *                                                              *
 *--------------------------------------------------------------*/
static Word16 Chebps_11(Word16 x, Word16 f[], Word16 n)
{
  Word16 i, cheb;
  Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
  Word32 t0;

 /* Note: All computation are done in Q24. */

  b2_h = 256;                           /* b2 = 1.0 in Q24 DPF */
  b2_l = 0;

  t0 = L_mult(x, 512);                  /* 2*x in Q24          */
  t0 = L_mac(t0, f[1], 4096);           /* + f[1] in Q24       */
  L_Extract(t0, &b1_h, &b1_l);          /* b1 = 2*x + f[1]     */

  for (i = 2; i<n; i++)
  {
    t0 = Mpy_32_16(b1_h, b1_l, x);      /* t0 = 2.0*x*b1              */
    t0 = L_shl(t0, 1);
    t0 = L_mac(t0,b2_h,(Word16)-32768L);/* t0 = 2.0*x*b1 - b2         */
    t0 = L_msu(t0, b2_l, 1);
    t0 = L_mac(t0, f[i], 4096);         /* t0 = 2.0*x*b1 - b2 + f[i]; */

    L_Extract(t0, &b0_h, &b0_l);        /* b0 = 2.0*x*b1 - b2 + f[i]; */

    b2_l = b1_l;                        /* b2 = b1; */
    b2_h = b1_h;
    b1_l = b0_l;                        /* b1 = b0; */
    b1_h = b0_h;
  }

  t0 = Mpy_32_16(b1_h, b1_l, x);        /* t0 = x*b1;              */
  t0 = L_mac(t0, b2_h,(Word16)-32768L); /* t0 = x*b1 - b2          */
  t0 = L_msu(t0, b2_l, 1);
  t0 = L_mac(t0, f[i], 2048);           /* t0 = x*b1 - b2 + f[i]/2 */

  t0 = L_shl(t0, 6);                    /* Q24 to Q30 with saturation */
  cheb = extract_h(t0);                 /* Result in Q14              */


  return(cheb);
}


static Word16 Chebps_10(Word16 x, Word16 f[], Word16 n)
{
  Word16 i, cheb;
  Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
  Word32 t0;

 /* Note: All computation are done in Q23. */

  b2_h = 128;                           /* b2 = 1.0 in Q23 DPF */
  b2_l = 0;

  t0 = L_mult(x, 256);                  /* 2*x in Q23          */
  t0 = L_mac(t0, f[1], 4096);           /* + f[1] in Q23       */
  L_Extract(t0, &b1_h, &b1_l);          /* b1 = 2*x + f[1]     */

  for (i = 2; i<n; i++)
  {
    t0 = Mpy_32_16(b1_h, b1_l, x);      /* t0 = 2.0*x*b1              */
    t0 = L_shl(t0, 1);
    t0 = L_mac(t0,b2_h,(Word16)-32768L);/* t0 = 2.0*x*b1 - b2         */
    t0 = L_msu(t0, b2_l, 1);
    t0 = L_mac(t0, f[i], 4096);         /* t0 = 2.0*x*b1 - b2 + f[i]; */

    L_Extract(t0, &b0_h, &b0_l);        /* b0 = 2.0*x*b1 - b2 + f[i]; */

    b2_l = b1_l;                        /* b2 = b1; */
    b2_h = b1_h;
    b1_l = b0_l;                        /* b1 = b0; */
    b1_h = b0_h;
  }

  t0 = Mpy_32_16(b1_h, b1_l, x);        /* t0 = x*b1;              */
  t0 = L_mac(t0, b2_h,(Word16)-32768L); /* t0 = x*b1 - b2          */
  t0 = L_msu(t0, b2_l, 1);
  t0 = L_mac(t0, f[i], 2048);           /* t0 = x*b1 - b2 + f[i]/2 */

  t0 = L_shl(t0, 7);                    /* Q23 to Q30 with saturation */
  cheb = extract_h(t0);                 /* Result in Q14              */


  return(cheb);
}