g729ev_tdac_mdct.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
*/

#include <stdlib.h>

#include "stl.h"
#include "G729EV_MAIN_defines.h"
#include "G729EV_TDAC_tfr.h"
#include "G729EV_TDAC_mdct.h"

/*--------------------------------------------------------------------------*
 *  Function  G729EV_TDAC_mdct()                                            *
 *  ~~~~~~~~~~~~~~~~~~~~~~~~~~                                              *
 *  Compute MDCT                                                            *
 * The MDCT is computed using the fast algorithm described in               *
 * - "A fast algorithm for the implementation of Filter Banks based on Time *
 * Domain Aliasing Cancellation", P. Duhamel et al., ICASSP91, pp 2209-2212 *
 * Petit, Conference ICASSP 1991, pp 2209 - 2212,                           *
 *                                                                          *
 * note about MDCT normalization factor                                     *
 *   -if norm_shift = -15, input samples are windowed to compute norm_shift *
 *    coefficient  but IFFT is not performed. norm_shift is then returned   *
 *   -if norm_shift != -15, mdct coefficients are computed using this norm  *
 *--------------------------------------------------------------------------*/
Word16 G729EV_TDAC_mdct(Word16 * mem,     /* (i)     old input samples */
                        Word16 * input,   /* (i)     input samples */
                        Word16 * ykr,     /* (o)     MDCT coefficients */
                        Word16 norm_shift /* (i)     normalization factor */
    )
{
  Word32    ACC0;               /* ACC */
  Word32    ACC1;               /* ACC */

  Word16    xin[G729EV_TDAC_L_WIN];
  Word16    ycr[G729EV_TDAC_L_WIN4];
  Word16    yci[G729EV_TDAC_L_WIN4];

  const Word16 *ptr_h1;         /* pointer on window */
  const Word16 *ptr_h2;         /* pointer on window */
  Word16   *ptr_x1;             /* pointer on input samples */
  Word16   *ptr_x2;             /* pointer on input samples */
  Word16   *ptr_ycr;            /* pointer on ycr */
  Word16   *ptr_yci;            /* pointer on yci */

  Word16    k;
  Word16    i;
  Word16    tmp16_ycr;
  Word16    sign;               /* sign for FFT/IFFT */

  /* variables for norm_shift evaluation */
  Word16    temp;
  Word16    temp1;
  Word16    temp2;
  Word16    temp3;

  /********************************************************************************/
  /* MDCT Computation                                                             */
  /********************************************************************************/

  FOR(i = 0; i < G729EV_TDAC_L_WIN2; i++)
  {
    xin[i] = mem[i];
#if(WMOPS)
    move16();
#endif
  }

  FOR(i = 0; i < G729EV_TDAC_L_WIN2; i++)
  {
    xin[i + G729EV_TDAC_L_WIN2] = input[i];
#if(WMOPS)
    move16();
#endif
  }

  IF(add(norm_shift, 15) == 0)  /* Evaluate norm_shift */
  {
    /* 1 --> computing Xc[n] (complex terms) */
    ptr_h1 = G729EV_TDAC_h + G729EV_TDAC_L_WIN2 - 1;  /* Middle of the window */
    ptr_x1 = xin + G729EV_TDAC_L_WIN2;

    ptr_h2 = G729EV_TDAC_h + G729EV_TDAC_L_WIN; /* End of the window */
    ptr_x2 = xin + G729EV_TDAC_L_WIN;
    k = (Word16) 0;

#if (WMOPS)
    move16();
#endif

    ptr_yci = yci;
    ptr_ycr = ycr;

    FOR(i = 0; i < G729EV_TDAC_L_WIN4; i++)
    {
      ACC0 = L_mult(G729EV_TDAC_h[k], xin[k]);
      ACC0 = L_msu(ACC0, ptr_h1[sub(-1, k)], ptr_x1[sub(-1, k)]);
      *ptr_ycr++ = round(ACC0);

      ACC0 = L_mult(ptr_h2[sub(-1, k)], ptr_x2[sub(-1, k)]);
      ACC0 = L_mac(ACC0, ptr_h1[k], ptr_x1[k]);
      *ptr_yci++ = round(ACC0);

      k = add(k, 2);
    }

    temp1 = (Word16) 0;
#if(WMOPS)
    move16();
#endif

    FOR(i = 0; i < G729EV_TDAC_L_WIN4; i++)
    {
      temp2 = abs_s(ycr[i]);    /* real part */
      temp3 = abs_s(yci[i]);    /* imaginary part */
      temp = sub(temp3, temp2); /* keep bigger value between R and I */

      if (temp > 0)             /* temp2 = max(temp2, temp3) */
      {
        temp2 = temp3;
#if (WMOPS)
        move16();
#endif
      }

      temp = sub(temp1, temp2); /* temp1 = max(temp1, temp2) */
      if (temp < 0)
      {
        temp1 = temp2;
#if (WMOPS)
        move16();
#endif
      }
    }

    temp = sub(temp1, 14000);
    IF(temp >= 0)
    {
      norm_shift = (Word16) 0;
#if (WMOPS)
      move16();
#endif
    }
    ELSE
    {
      ACC1 = L_mult(temp1, 9587);
      ACC1 = L_shr(ACC1, 20);
      temp2 = extract_l(ACC1);
      temp = norm_s(temp2);
      IF(temp == 0)
      {
        norm_shift = (Word16) 8;
#if (WMOPS)
        move16();
#endif
      }
      ELSE
      {
        norm_shift = sub(temp, 6);
      }

    }
    return norm_shift;          /* Quit function and return norm_shift value */
  }
  ELSE                          /* MDCT operations */
  {

    /* 1 --> computing Xc[n] (complex terms) (divided by two to avoid overflow) */

    ptr_h1 = G729EV_TDAC_h + G729EV_TDAC_L_WIN2;  /* Middle of the window */
    ptr_x1 = xin + G729EV_TDAC_L_WIN2;

    ptr_h2 = G729EV_TDAC_h + G729EV_TDAC_L_WIN; /* End of the window */
    ptr_x2 = xin + G729EV_TDAC_L_WIN;

    k = (Word16) 0;

#if (WMOPS)
    move16();
#endif

    ptr_yci = yci;
    ptr_ycr = ycr;

    FOR(i = 0; i < G729EV_TDAC_L_WIN4; i++)
    {
      ACC0 = L_mult(G729EV_TDAC_h[k], xin[k]);
      ACC0 = L_msu(ACC0, ptr_h1[sub(-1, k)], ptr_x1[sub(-1, k)]);
      *ptr_ycr++ = round(ACC0);

      ACC0 = L_mult(ptr_h2[sub(-1, k)], ptr_x2[sub(-1, k)]);
      ACC0 = L_mac(ACC0, ptr_h1[k], ptr_x1[k]);
      *ptr_yci++ = round(ACC0);

      k = add(k, 2);
    }

    /* 2 --> Product of Xc[n] by W**i + norm_shift */

    ptr_yci = yci;
    ptr_ycr = ycr;

    FOR(k = 0; k < G729EV_TDAC_L_WIN4; k++)
    {
      tmp16_ycr = *ptr_ycr;

#if (WMOPS)
      move16();
#endif

      ACC0 = L_mult(tmp16_ycr, G729EV_TDAC_wcos[k]);
      ACC0 = L_msu(ACC0, yci[k], G729EV_TDAC_wsin[k]);
      ACC0 = L_shl(ACC0, sub(norm_shift, 1));
      *ptr_ycr++ = round(ACC0);


      ACC0 = L_mult(tmp16_ycr, G729EV_TDAC_wsin[k]);
      ACC0 = L_mac(ACC0, yci[k], G729EV_TDAC_wcos[k]);
      ACC0 = L_shl(ACC0, sub(norm_shift, 1));
      *ptr_yci++ = round(ACC0);
    }

    /* 3 --> Inverse TFD */

    sign = (Word16) - 1;
#if (WMOPS)
    move16();
#endif

    G729EV_TDAC_tfr(ycr, yci, sign);

    /* 4 --> Final product and rearranging the results */

    ptr_x1 = ykr;
    ptr_x2 = ykr + G729EV_TDAC_L_WIN2 - 1;

    FOR(k = 0; k < G729EV_TDAC_L_WIN4; k++)
    {
      tmp16_ycr = ycr[k];

#if (WMOPS)
      move16();
#endif
      /* symetry of coeff k-1 and N-k */

      ACC0 = L_mult(yci[k], G729EV_TDAC_weti[k]);
      ACC0 = L_msu(ACC0, tmp16_ycr, G729EV_TDAC_wetr[k]);
      ACC0 = L_shr(ACC0, 1);
      *ptr_x2-- = round(ACC0);
      ptr_x2--;

      ACC0 = L_mult(tmp16_ycr, G729EV_TDAC_weti[k]);
      ACC0 = L_mac(ACC0, yci[k], G729EV_TDAC_wetr[k]);
      ACC0 = L_shr(ACC0, 1);
      *ptr_x1++ = round(ACC0);
      ptr_x1++;
    }

    FOR(i = 0; i < G729EV_TDAC_L_WIN2; i++)
    {
      mem[i] = input[i];
#if(WMOPS)
      move16();
#endif
    }
  }
  return 0;

}                               /* END MDCT */


/*--------------------------------------------------------------------------*
 *  Function  G729EV_TDAC_inv_mdct()                                        *
 *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                                          *
 *  Compute inverse MDCT and overlap-add                                    *
 * The underlying algorithm is described in                                 *
 * - "A fast algorithm for the implementation of Filter Banks based on Time *
 * Domain Aliasing Cancellation", P. Duhamel et al., ICASSP91, pp 2209-2212 *
 * Petit, Conference ICASSP 1991, pp 2209 - 2212,                           *
 *--------------------------------------------------------------------------*/
void G729EV_TDAC_inv_mdct(Word16 * xr,      /* (o)     output samples */
                          Word16 * ykq,     /* (i)     MDCT coefficients */
                          Word16 * ycim1,   /* (i)     previous MDCT memory */
                          Word16 norm_shift /* (i)     norm_shift value defined by coder */
    )
{
  Word32    ACC0;

  Word16    ycr[G729EV_TDAC_L_WIN4];
  Word16    yci[G729EV_TDAC_L_WIN4];
  Word16    sig_cur[G729EV_TDAC_L_WIN2];
  Word16    sig_next[G729EV_TDAC_L_WIN2];

  Word16   *ptr_yci;
  Word16   *ptr_ycr;
  Word16   *ptr1;
  Word16   *ptr2;
  Word16   *ptr1_next;
  Word16   *ptr2_next;
  const Word16 *ptr_h;          /* Pointer on window */

  Word16    k;
  Word16    sign;               /* sign for FFT / IFFT */
  Word16    tmp16;

    /*******************************************************************************/
  /* Inverse MDCT computation                                                    */
  /*******************************************************************************/

  /* 1 --> Input rotation = Product by wetrm1 */

  ptr1 = ykq;
  ptr2 = ykq + (G729EV_TDAC_L_WIN2 - 1);
  ptr_yci = yci;
  ptr_ycr = ycr;

  FOR(k = 0; k < G729EV_TDAC_L_WIN4; k++)
  {
    ACC0 = L_mult(*ptr2, G729EV_TDAC_wetrm1[k]);
    ACC0 = L_negate(ACC0);
    ACC0 = L_msu(ACC0, *ptr1, G729EV_TDAC_wetim1[k]);
    ACC0 = L_shl(ACC0, 1);
    *ptr_ycr++ = round(ACC0);

    ACC0 = L_mult(*ptr1++, G729EV_TDAC_wetrm1[k]);
    ACC0 = L_msu(ACC0, *ptr2--, G729EV_TDAC_wetim1[k]);
    ACC0 = L_shl(ACC0, 1);
    *ptr_yci++ = round(ACC0);

    ptr1++;
    ptr2--;
  }

  /* 2 --> Forward FFT : size = 160 */

#if (WMOPS)
  move16();
#endif

  sign = (Word16) 1;
  G729EV_TDAC_tfr(ycr, yci, sign);

  /* 3 --> Output rotation : product by a complex exponent */

  ptr_yci = yci;
  ptr_ycr = ycr;

  FOR(k = 0; k < G729EV_TDAC_L_WIN4; k++)
  {
    tmp16 = *ptr_ycr;

#if (WMOPS)
    move16();
#endif

    ACC0 = L_mult(tmp16, G729EV_TDAC_wcos[k]);
    ACC0 = L_mac(ACC0, yci[k], G729EV_TDAC_wsin[k]);
    ACC0 = L_shr(ACC0, sub(norm_shift, 1));
    *ptr_ycr++ = round(ACC0);

    ACC0 = L_mult(yci[k], G729EV_TDAC_wcos[k]);
    ACC0 = L_msu(ACC0, tmp16, G729EV_TDAC_wsin[k]);
    ACC0 = L_shr(ACC0, sub(norm_shift, 1));
    *ptr_yci++ = round(ACC0);
  }

  /* 4 --> Overlap and windowing (in one step) - equivalent to complex product */

  ptr1 = sig_cur;
  ptr2 = sig_cur + G729EV_TDAC_L_WIN2 - 1;
  ptr1_next = sig_next;
  ptr2_next = sig_next + G729EV_TDAC_L_WIN2 - 1;

  FOR(k = 0; k < G729EV_TDAC_L_WIN4; k++)
  {
    *ptr1++ = ycr[k];
    *ptr2-- = -ycr[k];
    *ptr1_next++ = yci[k];
    *ptr2_next-- = yci[k];
    *ptr1++;
    ptr1_next++;
    ptr2--;
    ptr2_next--;

#if (WMOPS)
    move16();
    move16();
    move16();
    move16();
#endif
  }

  ptr_h = G729EV_TDAC_h + G729EV_TDAC_L_WIN2;
  FOR(k = 0; k < G729EV_TDAC_L_WIN2; k++)
  {
    ACC0 = L_mult(sig_cur[k], G729EV_TDAC_h[k]);
    ACC0 = L_mac(ACC0, ycim1[k], *(--ptr_h));
    ACC0 = L_shl(ACC0, 1);
    xr[k] = round(ACC0);

    ycim1[k] = sig_next[k];
#if (WMOPS)
    move16();
    move16();
#endif

  }
  return;
}