/*
   ITU-T G.729 Speech Coder with Annex B    ANSI-C Source Code
   Version 1.3    Last modified: August 1997

   Copyright (c) 1996,
   AT&T, France Telecom, NTT, Universite de Sherbrooke, Lucent Technologies,
   Rockwell International
   All rights reserved.
*/

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

#include "typedef.h"
#include "basic_op.h"
#include "oper_32b.h"

#include "ld8k.h"
#include "tab_ld8k.h"

void Autocorr(
  Word16 x[],      /* (i)    : Input signal                      */
  Word16 m,        /* (i)    : LPC order                         */
  Word16 r_h[],    /* (o)    : Autocorrelations  (msb)           */
  Word16 r_l[],    /* (o)    : Autocorrelations  (lsb)           */
  Word16 *exp_R0
)
{
  Word16 i, j, norm;
  Word16 y[L_WINDOW];
  Word32 sum;

  extern Flag Overflow;

  /* Windowing of signal */

  for(i=0; i<L_WINDOW; i++)
  {
    y[i] = mult_r(x[i], hamwindow[i]);
  }

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

  /* should exp_R0 be initialized to 2 ??? */
  *exp_R0 = 1;    

  do {
    Overflow = 0;
    sum = 1;                   /* Avoid case of all zeros */
    for(i=0; i<L_WINDOW; i++)
      sum = L_mac(sum, y[i], y[i]);

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

    if(Overflow != 0)
    {
      for(i=0; i<L_WINDOW; i++)
      {
        y[i] = shr(y[i], 2);
      }
      *exp_R0 = add((*exp_R0), 4);
      Overflow = 1;
    }
  }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) */
  *exp_R0 = sub(*exp_R0, norm);

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

  for (i = 1; i <= m; i++)
  {
    sum = 0;
    for(j=0; j<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 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 Lag_window(
  Word16 m,         /* (i)     : LPC order                        */
  Word16 r_h[],     /* (i/o)   : Autocorrelations  (msb)          */
  Word16 r_l[]      /* (i/o)   : Autocorrelations  (lsb)          */
)
{
  Word16 i;
  Word32 x;

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

static Word16 old_A[M+1]={4096,0,0,0,0,0,0,0,0,0,0};
static Word16 old_rc[2]={0,0};

void 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 *Err       /* (o)     : Residual energy                          */
)
{
 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         */
 Word16 Ah[M+1], Al[M+1];      /* LPC coef. in double prec.                   */
 Word16 Anh[M+1], Anl[M+1];    /* LPC coef.for next iteration in double prec. */
 Word32 t0, t1, t2;            /* temporary variable                          */


/* 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;
  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 M                 *
 *--------------------------------------*/

  for(i= 2; i<=M; i++)
  {

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

    t0 = 0;
    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;

    /* Test for unstable filter. If unstable keep old A(z) */

    if (sub(abs_s(Kh), 32750) > 0)
    {
      for(j=0; j<=M; j++)
      {
        A[j] = old_A[j];
      }
      rc[0] = old_rc[0];        /* only two rc coefficients are needed */
      rc[1] = old_rc[1];
      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];
    }
  }

  *Err = shr(alp_h, alp_exp);

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

  A[0] = 4096;
  for(i=1; i<=M; i++)
  {