g729ev_g729_lpc.c

最新的ITU-T的宽带语音编解码标准G.729.1,是对原先的G.729的最好的调整.码流输出速率可以进行自适应调整.满足未来通信要求.希望对大家有所帮助.

/* ITU-T G.729EV Optimization/Characterization Candidate                         */
/* Version:       1.0.a                                                          */
/* Revision Date: June 28, 2006                                                  */

/*
   ITU-T G.729EV Optimization/Characterization Candidate    ANSI-C Source Code
   Copyright (c) 2006 
   France Telecom, Matsushita Electric, Mindspeed, Siemens AG, ETRI, VoiceAge Corp.
   All rights reserved
*/

/*-----------------------------------------------------*
 * Function G729EV_G729_AutocorrLPC()                  *
 *                                                     *
 *   Compute autocorrelations of signal with windowing *
 *                                                     *
 *-----------------------------------------------------*/

#include "stl.h"
#include "G729EV_MAIN_OPER_32B.h"

#include "G729EV_G729_ld8k.h"
#include "G729EV_G729_TAB_ld8k.h"

#ifdef WMOPS
#include "count.h"
#endif

void G729EV_G729_AutocorrLPC(Word16 x[],  /* (i)    : Input signal                      */
                             Word16 m,    /* (i)    : LPC order                         */
                             Word16 r_h[],/* (o)    : Autocorrelations  (msb)           */
                             Word16 r_l[] /* (o)    : Autocorrelations  (lsb)           */
    )
{
  Word32    sum;

  Word16    y[G729EV_G729_L_WINDOW];

  Word16    i, j, norm;

  extern Flag Overflow;

  /* Windowing of signal */

  FOR(i = 0; i < G729EV_G729_L_WINDOW; i++)
  {
#ifdef WMOPS
    move16();
#endif
    y[i] = mult_r(x[i], hamwindow[i]);
  }

  /* Compute r[0] and test for overflow */

  DO
  {
    Overflow = 0;
    sum = 1;                    /* Avoid case of all zeros */
#ifdef WMOPS
    move16();
    move32();
#endif
    FOR(i = 0; i < G729EV_G729_L_WINDOW; i++)
      sum = L_mac(sum, y[i], y[i]);

    /* If overflow divide y[] by 4 */

    IF(Overflow != 0)
    {
      FOR(i = 0; i < G729EV_G729_L_WINDOW; i++)
      {
#ifdef WMOPS
        move16();
#endif
        y[i] = shr(y[i], 2);
      }
    }
  }
  WHILE(Overflow != 0);

  /* Normalization of r[0] */

  norm = norm_l(sum);
  sum = L_shl(sum, norm);
  L_Extract(sum, &r_h[0], &r_l[0]); /* Put in DPF format (see oper_32b) */

  /* r[1] to r[m] */

  FOR(i = 1; i <= m; i++)
  {
    sum = 0;
#ifdef WMOPS
    move32();
#endif
    FOR(j = 0; j < G729EV_G729_L_WINDOW - i; j++)
      sum = L_mac(sum, y[j], y[j + i]);

    sum = L_shl(sum, norm);
    L_Extract(sum, &r_h[i], &r_l[i]);
  }
  return;
}


/*-------------------------------------------------------*
 * Function G729EV_G729_Lag_window()                                 *
 *                                                       *
 * Lag_window on autocorrelations.                       *
 *                                                       *
 * r[i] *= lag_wind[i]                                   *
 *                                                       *
 *  r[i] and lag_wind[i] are in special double precision.*
 *  See "oper_32b.c" for the format                      *
 *                                                       *
 *-------------------------------------------------------*/

void G729EV_G729_Lag_window(Word16 m,     /* (i)     : LPC order                        */
                            Word16 r_h[], /* (i/o)   : Autocorrelations  (msb)          */
                            Word16 r_l[]  /* (i/o)   : Autocorrelations  (lsb)          */
    )
{
  Word32    x;

  Word16    i;


  FOR(i = 1; i <= m; i++)
  {
    x = Mpy_32(r_h[i], r_l[i], lag_h[i - 1], lag_l[i - 1]);
    L_Extract(x, &r_h[i], &r_l[i]);
  }
  return;
}


/*___________________________________________________________________________
 |                                                                           |
 |      LEVINSON-DURBIN algorithm in double precision                        |
 |      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                        |
 |---------------------------------------------------------------------------|
 |                                                                           |
 | Algorithm                                                                 |
 |                                                                           |
 |       R[i]    autocorrelations.                                           |
 |       A[i]    filter coefficients.                                        |
 |       K       reflection coefficients.                                    |
 |       Alpha   prediction gain.                                            |
 |                                                                           |
 |       Initialization:                                                     |
 |               A[0] = 1                                                    |
 |               K    = -R[1]/R[0]                                           |
 |               A[1] = K                                                    |
 |               Alpha = R[0] * (1-K**2]                                     |
 |                                                                           |
 |       Do for  i = 2 to M                                                  |
 |                                                                           |
 |            S =  SUM ( R[j]*A[i-j] ,j=1,i-1 ) +  R[i]                      |
 |                                                                           |
 |            K = -S / Alpha                                                 |
 |                                                                           |
 |            An[j] = A[j] + K*A[i-j]   for j=1 to i-1                       |
 |                                      where   An[i] = new A[i]             |
 |            An[i]=K                                                        |
 |                                                                           |
 |            Alpha=Alpha * (1-K**2)                                         |
 |                                                                           |
 |       END                                                                 |
 |                                                                           |
 | Remarks on the dynamics of the calculations.                              |
 |                                                                           |
 |       The numbers used are in double precision in the following format :  |
 |       A = AH <<16 + AL<<1.  AH and AL are 16 bit signed integers.         |
 |       Since the LSB's also contain a sign bit, this format does not       |
 |       correspond to standard 32 bit integers.  We use this format since   |
 |       it allows fast execution of multiplications and divisions.          |
 |                                                                           |
 |       "DPF" will refer to this special format in the following text.      |
 |       See oper_32b.c                                                      |
 |                                                                           |
 |       The R[i] were normalized in routine AUTO (hence, R[i] < 1.0).       |
 |       The K[i] and Alpha are theoretically < 1.0.                         |
 |       The A[i], for a sampling frequency of 8 kHz, are in practice        |
 |       always inferior to 16.0.                                            |
 |                                                                           |
 |       These characteristics allow straigthforward fixed-point             |
 |       implementation.  We choose to represent the parameters as           |
 |       follows :                                                           |
 |                                                                           |
 |               R[i]    Q31   +- .99..                                      |
 |               K[i]    Q31   +- .99..                                      |
 |               Alpha   Normalized -> mantissa in Q31 plus exponent         |
 |               A[i]    Q27   +- 15.999..                                   |
 |                                                                           |
 |       The additions are performed in 32 bit.  For the summation used      |
 |       to calculate the K[i], we multiply numbers in Q31 by numbers        |
 |       in Q27, with the result of the multiplications in Q27,              |
 |       resulting in a dynamic of +- 16.  This is sufficient to avoid       |
 |       overflow, since the final result of the summation is                |
 |       necessarily < 1.0 as both the K[i] and Alpha are                    |
 |       theoretically < 1.0.                                                |
 |___________________________________________________________________________|
*/


/* Last A(z) for case of unstable filter */

void G729EV_G729_Levinson(Word16 Rh[],  /* (i)     : Rh[M+1] Vector of autocorrelations (msb) */
                          Word16 Rl[],  /* (i)     : Rl[M+1] Vector of autocorrelations (lsb) */
                          Word16 A[],   /* (o) Q12 : A[M]    LPC coefficients  (m = 10)       */
                          Word16 rc[],  /* (o) Q15 : rc[M]   Reflection coefficients.         */
                          Word16 old_A[], Word16 old_rc[])
{
  Word32    t0, t1, t2;             /* temporary variable                          */

  Word16    Ah[G729EV_G729_M + 1], Al[G729EV_G729_M + 1]; /* LPC coef. in double prec.                   */
  Word16    Anh[G729EV_G729_M + 1], Anl[G729EV_G729_M + 1]; /* LPC coef.for next iteration in double prec. */

  Word16    i, j;
  Word16    hi, lo;
  Word16    Kh, Kl;                 /* reflection coefficient; hi and lo           */
  Word16    alp_h, alp_l, alp_exp;  /* Prediction gain; hi lo and exponent         */



/* K = A[1] = -R[1] / R[0] */

  t1 = L_Comp(Rh[1], Rl[1]);    /* R[1] in Q31      */
  t2 = L_abs(t1);               /* abs R[1]         */
  t0 = Div_32(t2, Rh[0], Rl[0]);/* R[1]/R[0] in Q31 */
  if (t1 > 0)
    t0 = L_negate(t0);          /* -R[1]/R[0]       */
  L_Extract(t0, &Kh, &Kl);      /* K in DPF         */
  rc[0] = Kh;
#ifdef WMOPS
  move16();
#endif
  t0 = L_shr(t0, 4);            /* A[1] in Q27      */
  L_Extract(t0, &Ah[1], &Al[1]);/* A[1] in DPF      */

/*  Alpha = R[0] * (1-K**2) */

  t0 = Mpy_32(Kh, Kl, Kh, Kl);  /* K*K      in Q31 */
  t0 = L_abs(t0);               /* Some case <0 !! */
  t0 = L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K  in Q31 */
  L_Extract(t0, &hi, &lo);      /* DPF format      */
  t0 = Mpy_32(Rh[0], Rl[0], hi, lo);  /* Alpha in Q31    */

/* Normalize Alpha */

  alp_exp = norm_l(t0);
  t0 = L_shl(t0, alp_exp);
  L_Extract(t0, &alp_h, &alp_l);  /* DPF format    */

/*--------------------------------------*
 * ITERATIONS  I=2 to G729EV_G729_M                 *
 *--------------------------------------*/

  FOR(i = 2; i <= G729EV_G729_M; i++)
  {

    /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) +  R[i] */

    t0 = 0;
#ifdef WMOPS
    move32();
#endif
    FOR(j = 1; j < i; j++)
      t0 = L_add(t0, Mpy_32(Rh[j], Rl[j], Ah[i - j], Al[i - j]));

    t0 = L_shl(t0, 4);          /* result in Q27 -> convert to Q31 */
    /* No overflow possible            */
    t1 = L_Comp(Rh[i], Rl[i]);
    t0 = L_add(t0, t1);         /* add R[i] in Q31                 */

    /* K = -t0 / Alpha */

    t1 = L_abs(t0);
    t2 = Div_32(t1, alp_h, alp_l);  /* abs(t0)/Alpha                   */
    if (t0 > 0)
      t2 = L_negate(t2);            /* K =-t0/Alpha                    */
    t2 = L_shl(t2, alp_exp);        /* denormalize; compare to Alpha   */
    L_Extract(t2, &Kh, &Kl);        /* K in DPF                        */
    rc[i - 1] = Kh;
#ifdef WMOPS
    move16();
#endif

    /* Test for unstable filter. If unstable keep old A(z) */
    IF(sub(abs_s(Kh), 32750) > 0)
    {
      FOR(j = 0; j <= G729EV_G729_M; j++)
      {
        A[j] = old_A[j];
#ifdef WMOPS
        move16();
#endif
      }
      rc[0] = old_rc[0];        /* only two rc coefficients are needed */
      rc[1] = old_rc[1];
#ifdef WMOPS
      move16();
      move16();
#endif
      return;
    }

    /*------------------------------------------*
     *  Compute new LPC coeff. -> An[i]         *
     *  An[j]= A[j] + K*A[i-j]     , j=1 to i-1 *
     *  An[i]= K                                *
     *------------------------------------------*/


    FOR(j = 1; j < i; j++)
    {
      t0 = Mpy_32(Kh, Kl, Ah[i - j], Al[i - j]);
      t0 = L_add(t0, L_Comp(Ah[j], Al[j]));
      L_Extract(t0, &Anh[j], &Anl[j]);
    }
    t2 = L_shr(t2, 4);                /* t2 = K in Q31 ->convert to Q27  */
    L_Extract(t2, &Anh[i], &Anl[i]);  /* An[i] in Q27                    */

    /*  Alpha = Alpha * (1-K**2) */

    t0 = Mpy_32(Kh, Kl, Kh, Kl);  /* K*K      in Q31 */
    t0 = L_abs(t0);               /* Some case <0 !! */
    t0 = L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K  in Q31 */
    L_Extract(t0, &hi, &lo);      /* DPF format      */
    t0 = Mpy_32(alp_h, alp_l, hi, lo);  /* Alpha in Q31    */

    /* Normalize Alpha */

    j = norm_l(t0);
    t0 = L_shl(t0, j);
    L_Extract(t0, &alp_h, &alp_l);  /* DPF format    */
    alp_exp = add(alp_exp, j);      /* Add normalization to alp_exp */

    /* A[j] = An[j] */

    FOR(j = 1; j <= i; j++)
    {
      Ah[j] = Anh[j];
      Al[j] = Anl[j];
#ifdef WMOPS
      move16();
      move16();
#endif
    }
  }

  /* Truncate A[i] in Q27 to Q12 with rounding */

  A[0] = 4096;
#ifdef WMOPS
  move16();
#endif
  FOR(i = 1; i <= G729EV_G729_M; i++)
  {
    t0 = L_Comp(Ah[i], Al[i]);
#ifdef WMOPS
    move16();
    move16();
#endif
    old_A[i] = A[i] = round(L_shl(t0, 1));
  }
  old_rc[0] = rc[0];
  old_rc[1] = rc[1];
#ifdef WMOPS
  move16();
  move16();
#endif

  return;
}



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

/* local function */


void G729EV_G729_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) */
    )
{
  Word32    t0, L_temp;

  Word16    f1[G729EV_G729_M / 2 + 1], f2[G729EV_G729_M / 2 + 1];

  Word16(*pChebps) (Word16 x, Word16 f[], Word16 n);

  Word16   *coef;

  Word16    i, j, nf, ip;
  Word16    xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
  Word16    x, y, sign, exp;
  Flag      ovf_coef;