/***************************************************************************
 *
 *   File Name:  dec_func.c
 *
 *   Derivation:
 *	This library module was derived from sp_dec.c in the original
 *	GSM 06.06 code drop.  Compared to the original, this version
 *	excludes the main speechDecoder() function and two static
 *	functions that make use of speech decoder state, leaving only
 *	those functions that can be used by both decoder and encoder
 *	engines and are not tied to either state structure.
 *
 *   Original comments follow:
 *
 *   Purpose:
 *      Contains all functions for decoding speech.  It does not
 *      include those routines needed to decode channel information.
 *
 *      Since the GSM half-rate speech coder is an analysis-by-synthesis
 *      coder, many of the routines in this file are also called by the
 *      encoder.  Functions are included for coded-parameter lookup,
 *      LPC filter coefficient interpolation, excitation vector lookup
 *      and construction, vector quantized gain lookup, and LPC synthesis
 *      filtering.  In addition, some post-processing functions are
 *      included.
 *
 *     Below is a listing of all the functions appearing in the file.
 *     The functions are arranged according to their purpose.  Under
 *     each heading, the ordering is hierarchical.
 *
 *     The entire speech decoder, under which all these routines fall,
 *     except were noted:
 *     speechDecoder()
 *
 *     Spectral Smoothing of LPC:
 *       a_sst()
 *         aFlatRcDp()
 *         rcToCorrDpL()
 *         aToRc()
 *         rcToADp()
 *     VSELP codevector construction:
 *       b_con()
 *       v_con()
 *     LTP vector contruction:
 *       fp_ex()
 *         get_ipjj()
 *       lagDecode()
 *     LPC contruction
 *       getSfrmLpc()
 *         interpolateCheck()
 *         res_eng()
 *       lookupVq()
 *     Excitation scaling:
 *       rs_rr()
 *         g_corr1() (no scaling)
 *       rs_rrNs()
 *         g_corr1s() (g_corr1 with scaling)
 *       scaleExcite()
 *     Post filtering:
 *       pitchPreFilt()
 *         agcGain()
 *         lpcIir()
 *       r0BasedEnergyShft()
 *       spectralPostFilter()
 *         lpcFir()
 *
 *
 *     Routines not referenced by speechDecoder()
 *     Filtering routines:
 *       lpcIrZsIir()
 *       lpcZiIir()
 *       lpcZsFir()
 *       lpcZsIir()
 *       lpcZsIirP()
 *     Square root:
 *       sqroot()
 *
 **************************************************************************/

/*_________________________________________________________________________
 |                                                                         |
 |                            Include Files                                |
 |_________________________________________________________________________|
*/

#include "namespace.h"
#include "tw_gsmhr.h"
#include "typedefs.h"
#include "mathhalf.h"
#include "sp_rom.h"
#include "dec_func.h"


/*_________________________________________________________________________
 |                                                                         |
 |                              Local Defines                              |
 |_________________________________________________________________________|
*/

#define  P_INT_MACS   10
#define  ASCALE       0x0800
#define  ASHIFT       4
#define  DELTA_LEVELS 16
#define  GSP0_SCALE   1
#define  C_BITS_V     9                /* number of bits in any voiced VSELP
                                        * codeword */
#define  C_BITS_UV    7                /* number of bits in a unvoiced VSELP
                                        * codeword */
#define  MAXBITS      C_BITS_V         /* max number of bits in any VSELP
                                        * codeword */
#define  LTP_LEN      147              /* 147==0x93 length of LTP history */
#define  SQRT_ONEHALF 0x5a82           /* the 0.5 ** 0.5 */
#define  LPC_ROUND    0x00000800L      /* 0x8000 >> ASHIFT */
#define  AFSHIFT      2                /* number of right shifts to be
                                        * applied to the autocorrelation
                                        * sequence in aFlatRcDp     */

/***************************************************************************
 *
 *   FUNCTION NAME: aFlatRcDp
 *
 *   PURPOSE:
 *
 *     Given a Longword autocorrelation sequence, representing LPC
 *     information, aFlatRcDp converts the vector to one of NP
 *     Shortword reflection coefficients.
 *
 *   INPUT:
 *
 *
 *     pL_R[0:NP]    - An input Longword autocorrelation sequence, (pL_R[0] =
 *                     not necessarily 0x7fffffffL).  pL_R is altered in the
 *                     call, by being right shifted by global constant
 *                     AFSHIFT bits.
 *
 *                     The input array pL_R[] should be shifted left as much
 *                     as possible to improve precision.
 *
 *     AFSHIFT       - The number of right shifts to be applied to the
 *                     normalized autocorrelation sequence pL_R.
 *
 *   OUTPUT:
 *
 *     pswRc[0:NP-1] - A Shortword output vector of NP reflection
 *                     coefficients.
 *
 *   RETURN VALUE:
 *
 *     None
 *
 *   DESCRIPTION:
 *
 *     This routine transforms LPC information from one set of
 *     parameters to another.  It is better suited for fixed point
 *     implementations than the Levinson-Dubin recursion.
 *
 *     The function is called by a_sst(), and getNWCoefs().  In a_sst()
 *     direct form coefficients are converted to autocorrelations,
 *     and smoothed in that domain.  Conversion back to direct form
 *     coefficients is done by calling aFlatRc(), followed by rcToADp().
 *
 *     In getNwCoefs() again the conversion back to direct form
 *     coefficients is done by calling aFlatRc(), followed by rcToADp().
 *     In getNwCoefs() an autocorrelation sequence is generated from the
 *     impulse response of the weighting filters.
 *
 *     The fundamental recursion is derived from AFLAT, which is
 *     described in section 4.1.4.1.
 *
 *     Unlike in AFLAT where the reflection coefficients are known, here
 *     they are the unknowns.  PBar and VBar for j==0 are initially
 *     known, as is rSub1.  From this information the next set of P's
 *     and V's are generated.  At the end of the recursion the next,
 *     reflection coefficient rSubj (pswRc[j]) can be calcluated by
 *     dividing Vsubj by Psubj.
 *
 *     Precision is crucial in this routine.  At each stage, a
 *     normalization is performed prior to the reflection coefficient
 *     calculation.  In addition, to prevent overflow, the
 *     autocorrelation sequence is scaled down by ASHIFT (4) right
 *     shifts.
 *
 *
 *   REFERENCES: Sub_Clause 4.1.9 and 4.2.1  of GSM Recomendation 06.20
 *
 *   KEYWORDS: reflection coefficients, AFLAT, aflat, recursion, LPC
 *   KEYWORDS: autocorrelation
 *
 *************************************************************************/

  void   aFlatRcDp(Longword *pL_R, Shortword *pswRc)
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword pL_pjNewSpace[NP];
  Longword pL_pjOldSpace[NP];
  Longword pL_vjNewSpace[2 * NP - 1];
  Longword pL_vjOldSpace[2 * NP - 1];

  Longword *pL_pjOld;
  Longword *pL_pjNew;
  Longword *pL_vjOld;
  Longword *pL_vjNew;
  Longword *pL_swap;

  Longword L_temp;
  Longword L_sum;
  Shortword swRc,
         swRcSq,
         swTemp,
         swTemp1,
         swAbsTemp1,
         swTemp2;
  int    i,
         j;


/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  pL_pjOld = pL_pjOldSpace;
  pL_pjNew = pL_pjNewSpace;
  pL_vjOld = pL_vjOldSpace + NP - 1;
  pL_vjNew = pL_vjNewSpace + NP - 1;


  /* Extract the 0-th reflection coefficient */
  /*-----------------------------------------*/

  swTemp1 = round(pL_R[1]);
  swTemp2 = round(pL_R[0]);
  swAbsTemp1 = abs_s(swTemp1);
  if (swTemp2 <= 0 || sub(swAbsTemp1, swTemp2) >= 0)
  {
    j = 0;
    for (i = j; i < NP; i++)
    {
      pswRc[i] = 0;
    }
    return;
  }

  swRc = divide_s(swAbsTemp1, swTemp2);/* return division result */

  if (sub(swTemp1, swAbsTemp1) == 0)
    swRc = negate(swRc);               /* negate reflection Rc[j] */

  pswRc[0] = swRc;                     /* copy into the output Rc array */

  for (i = 0; i <= NP; i++)
  {
    pL_R[i] = L_shr(pL_R[i], AFSHIFT);
  }

  /* Initialize the pjOld and vjOld recursion arrays */
  /*-------------------------------------------------*/

  for (i = 0; i < NP; i++)
  {
    pL_pjOld[i] = pL_R[i];
    pL_vjOld[i] = pL_R[i + 1];
  }
  for (i = -1; i > -NP; i--)
    pL_vjOld[i] = pL_R[-(i + 1)];


  /* Compute the square of the j=0 reflection coefficient */
  /*------------------------------------------------------*/

  swRcSq = mult_r(swRc, swRc);

  /* Update pjNew and vjNew arrays for lattice stage j=1 */
  /*-----------------------------------------------------*/

  /* Updating pjNew: */
  /*-------------------*/

  for (i = 0; i <= NP - 2; i++)
  {
    L_temp = L_mpy_ls(pL_vjOld[i], swRc);
    L_sum = L_add(L_temp, pL_pjOld[i]);
    L_temp = L_mpy_ls(pL_pjOld[i], swRcSq);
    L_sum = L_add(L_temp, L_sum);
    L_temp = L_mpy_ls(pL_vjOld[-i], swRc);
    pL_pjNew[i] = L_add(L_sum, L_temp);
  }

  /* Updating vjNew: */
  /*-------------------*/

  for (i = -NP + 2; i <= NP - 2; i++)
  {
    L_temp = L_mpy_ls(pL_vjOld[-i - 1], swRcSq);
    L_sum = L_add(L_temp, pL_vjOld[i + 1]);
    L_temp = L_mpy_ls(pL_pjOld[(((i + 1) >= 0) ? i + 1 : -(i + 1))], swRc);
    L_temp = L_shl(L_temp, 1);
    pL_vjNew[i] = L_add(L_temp, L_sum);
  }



  j = 0;

  /* Compute reflection coefficients Rc[1],...,Rc[9] */
  /*-------------------------------------------------*/

  for (j = 1; j < NP; j++)
  {

    /* Swap pjNew and pjOld buffers */
    /*------------------------------*/

    pL_swap = pL_pjNew;
    pL_pjNew = pL_pjOld;
    pL_pjOld = pL_swap;

    /* Swap vjNew and vjOld buffers */
    /*------------------------------*/

    pL_swap = pL_vjNew;
    pL_vjNew = pL_vjOld;
    pL_vjOld = pL_swap;

    /* Compute the j-th reflection coefficient */
    /*-----------------------------------------*/

    swTemp = norm_l(pL_pjOld[0]);      /* get shift count */
    swTemp1 = round(L_shl(pL_vjOld[0], swTemp));        /* normalize num.  */
    swTemp2 = round(L_shl(pL_pjOld[0], swTemp));        /* normalize den.  */

    /* Test for invalid divide conditions: a) devisor < 0 b) abs(divident) >
     * abs(devisor) If either of these conditions is true, zero out
     * reflection coefficients for i=j,...,NP-1 and return. */

    swAbsTemp1 = abs_s(swTemp1);
    if (swTemp2 <= 0 || sub(swAbsTemp1, swTemp2) >= 0)
    {
      i = j;
      for (i = j; i < NP; i++)
      {
        pswRc[i] = 0;
      }
      return;
    }

    swRc = divide_s(swAbsTemp1, swTemp2);       /* return division result */
    if (sub(swTemp1, swAbsTemp1) == 0)
      swRc = negate(swRc);             /* negate reflection Rc[j] */
    swRcSq = mult_r(swRc, swRc);       /* compute Rc^2 */
    pswRc[j] = swRc;                   /* copy Rc[j] to output array */

    /* Update pjNew and vjNew arrays for the next lattice stage if j < NP-1 */
    /*---------------------------------------------------------------------*/

    /* Updating pjNew: */
    /*-----------------*/

    for (i = 0; i <= NP - j - 2; i++)
    {
      L_temp = L_mpy_ls(pL_vjOld[i], swRc);
      L_sum = L_add(L_temp, pL_pjOld[i]);
      L_temp = L_mpy_ls(pL_pjOld[i], swRcSq);
      L_sum = L_add(L_temp, L_sum);
      L_temp = L_mpy_ls(pL_vjOld[-i], swRc);
      pL_pjNew[i] = L_add(L_sum, L_temp);
    }

    /* Updating vjNew: */
    /*-----------------*/

    for (i = -NP + j + 2; i <= NP - j - 2; i++)
    {
      L_temp = L_mpy_ls(pL_vjOld[-i - 1], swRcSq);
      L_sum = L_add(L_temp, pL_vjOld[i + 1]);
      L_temp = L_mpy_ls(pL_pjOld[(((i + 1) >= 0) ? i + 1 : -(i + 1))], swRc);
      L_temp = L_shl(L_temp, 1);
      pL_vjNew[i] = L_add(L_temp, L_sum);
    }
  }
  return;
}

/***************************************************************************
 *
 *   FUNCTION NAME: aToRc
 *
 *   PURPOSE:
 *
 *     This subroutine computes a vector of reflection coefficients, given
 *     an input vector of direct form LPC filter coefficients.
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the LPC filter (global constant)
 *
 *     swAshift
 *                     The number of right shifts applied externally
 *                     to the direct form filter coefficients.
 *
 *     pswAin[0:NP-1]
 *                     The input vector of direct form LPC filter
 *                     coefficients.
 *
 *   OUTPUTS:
 *
 *     pswRc[0:NP-1]
 *                     Array containing the reflection coefficients.
 *
 *   RETURN VALUE:
 *
 *     siUnstableFlt
 *                     If stable reflection coefficients 0, 1 if unstable.
 *
 *
 *   DESCRIPTION:
 *
 *     This function performs the conversion from direct form
 *     coefficients to reflection coefficients. It is used in a_sst()
 *     and interpolateCheck().  In a_sst() reflection coefficients used
 *     as a transitional data format.  aToRc() is used for this
 *     conversion.
 *
 *     When performing interpolation, a stability check must be
 *     performed. interpolateCheck() does this by calling aToRc().
 *
 *     First coefficients are shifted down by iAshift. NP, the filter
 *     order is 10. The a's and rc's each have NP elements in them. An
 *     elaborate algorithm description can be found on page 443, of
 *     "Digital Processing of Speech Signals" by L.R. Rabiner and R.W.
 *     Schafer; Prentice-Hall; Englewood Cliffs, NJ (USA).  1978.
 *
 *   REFERENCES: Sub_Clause 4.1.6, and 4.2.3 of GSM Recomendation 06.20
 *
 *   KEYWORDS: reflectioncoefficients, parcors, conversion, atorc, ks, as
 *   KEYWORDS: parcorcoefficients, lpc, flat, vectorquantization
 *
 *************************************************************************/

short aToRc(Shortword swAshift, Shortword pswAin[], Shortword pswRc[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Constants                                    |
 |_________________________________________________________________________|
*/

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Shortword pswTmpSpace[NP],
         pswASpace[NP],
         swNormShift,
         swActShift,
         swNormProd,
         swRcOverE,
         swDiv,
        *pswSwap,
        *pswTmp,
        *pswA;

  Longword L_temp;

  short int siUnstableFlt,
         i,
         j;                            /* Loop control variables */

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Initialize starting addresses for temporary buffers */
  /*-----------------------------------------------------*/

  pswA = pswASpace;
  pswTmp = pswTmpSpace;

  /* Copy the direct form filter coefficients to a temporary array */
  /*---------------------------------------------------------------*/

  for (i = 0; i < NP; i++)
  {
    pswA[i] = pswAin[i];
  }

  /* Initialize the flag for filter stability check */
  /*------------------------------------------------*/

  siUnstableFlt = 0;

  /* Start computation of the reflection coefficients, Rc[9],...,Rc[1] */
  /*-------------------------------------------------------------------*/

  for (i = NP - 1; i >= 1; i--)
  {

    pswRc[i] = shl(pswA[i], swAshift); /* write Rc[i] to output array */

    /* Check the stability of i-th reflection coefficient */
    /*----------------------------------------------------*/

    siUnstableFlt = siUnstableFlt | isSwLimit(pswRc[i]);

    /* Precompute intermediate variables for needed for the computation */
    /* of direct form filter of order i-1                               */
    /*------------------------------------------------------------------*/

    if (sub(pswRc[i], SW_MIN) == 0)
    {
      siUnstableFlt = 1;
      swRcOverE = 0;
      swDiv = 0;
      swActShift = 2;
    }
    else
    {
      L_temp = LW_MAX;                 /* Load ~1.0 into accum */
      L_temp = L_msu(L_temp, pswRc[i], pswRc[i]);       /* 1.-Rc[i]*Rc[i]  */
      swNormShift = norm_l(L_temp);
      L_temp = L_shl(L_temp, swNormShift);
      swNormProd = extract_h(L_temp);
      swActShift = add(2, swNormShift);
      swDiv = divide_s(0x2000, swNormProd);
      swRcOverE = mult_r(pswRc[i], swDiv);
    }
    /* Check stability   */
    /*---------------------*/
    siUnstableFlt = siUnstableFlt | isSwLimit(swRcOverE);

    /* Compute direct form filter coefficients corresponding to */
    /* a direct form filter of order i-1                        */
    /*----------------------------------------------------------*/

    for (j = 0; j <= i - 1; j++)
    {
      L_temp = L_mult(pswA[j], swDiv);
      L_temp = L_msu(L_temp, pswA[i - j - 1], swRcOverE);
      L_temp = L_shl(L_temp, swActShift);
      pswTmp[j] = round(L_temp);
      siUnstableFlt = siUnstableFlt | isSwLimit(pswTmp[j]);
    }

    /* Swap swA and swTmp buffers */
    /*----------------------------*/

    pswSwap = pswA;
    pswA = pswTmp;
    pswTmp = pswSwap;
  }

  /* Compute reflection coefficient Rc[0] */
  /*--------------------------------------*/

  pswRc[0] = shl(pswA[0], swAshift);   /* write Rc[0] to output array */

  /* Check the stability of 0-th reflection coefficient */
  /*----------------------------------------------------*/

  siUnstableFlt = siUnstableFlt | isSwLimit(pswRc[0]);

  return (siUnstableFlt);
}

/***************************************************************************
 *
 *   FUNCTION NAME: b_con
 *
 *   PURPOSE:
 *     Expands codeword into an one dimensional array. The 0/1 input is
 *     changed to an element with magnitude +/- 0.5.
 *
 *     input  output
 *
 *       0    -0.5
 *       1    +0.5
 *
 *   INPUT:
 *
 *      swCodeWord
 *                     Input codeword, information in the LSB's
 *
 *      siNumBits
 *                     number of bits in the input codeword and number
 *                     of elements in output vector
 *
 *      pswVectOut[0:siNumBits]
 *
 *                     pointer to bit array
 *
 *   OUTPUT:
 *
 *      pswVectOut[0:siNumBits]
 *
 *                     signed bit array
 *
 *   RETURN:
 *
 *      none
 *
 *   REFERENCES: Sub_Clause 4.1.10 and 4.2.1 of GSM Recomendation 06.20
 *
 *   KEYWORDS: b_con, codeword, expansion
 *
 *************************************************************************/

void   b_con(Shortword swCodeWord, short siNumBits,
                    Shortword pswVectOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  short int siLoopCnt;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  for (siLoopCnt = 0; siLoopCnt < siNumBits; siLoopCnt++)
  {

    if (swCodeWord & 1)                /* temp accumulator get 0.5 */
      pswVectOut[siLoopCnt] = (Shortword) 0x4000;
    else                               /* temp accumulator gets -0.5 */
      pswVectOut[siLoopCnt] = (Shortword) 0xc000;

    swCodeWord = shr(swCodeWord, 1);
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: fp_ex
 *
 *   PURPOSE:
 *
 *     Looks up a vector in the adaptive excitation codebook (long-term
 *     predictor).
 *
 *   INPUTS:
 *
 *     swOrigLagIn
 *
 *                     Extended resolution lag (lag * oversampling factor)
 *
 *     pswLTPState[-147:39]
 *
 *                     Adaptive codebook (with space at end for looked up
 *                     vector).  Indicies [-147:-1] are the history, [0:39]
 *                     are for the looked up vector.
 *
 *     psrPitchIntrpFirBase[0:59]
 *     ppsrPVecIntFilt[0:9][0:5] ([tap][phase])
 *
 *                     Interpolating FIR filter coefficients.
 *
 *   OUTPUTS:
 *
 *     pswLTPState[0:39]
 *
 *                     Array containing the contructed output vector
 *
 *   RETURN VALUE:
 *     none
 *
 *   DESCRIPTION:
 *
 *     The adaptive codebook consists of the history of the excitation.
 *     The vector is looked up by going back into this history
 *     by the amount of the input lag.  If the input lag is fractional,
 *     then the samples to be looked up are interpolated from the existing
 *     samples in the history.
 *
 *     If the lag is less than the length of the vector to be generated
 *     (i.e. less than the subframe length), then the lag is doubled
 *     after the first n samples have been looked up (n = input lag).
 *     In this way, the samples being generated are not part of the
 *     codebook.  This is described in section 4.1.8.
 *
 *   REFERENCES: Sub_Clause 4.1.8.5 and 4.2.1  of GSM Recomendation 06.20
 *
 *   Keywords: pitch, excitation vector, long term filter, history,
 *   Keywords: fractional lag, get_ipjj
 *
 *************************************************************************/



void   fp_ex(Shortword swOrigLagIn,
                    Shortword pswLTPState[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Temp;
  Shortword swIntLag,
         swRemain,
         swRunningLag;
  short int siSampsSoFar,
         siSampsThisPass,
         i,
         j;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Loop: execute until all samples in the vector have been looked up */
  /*-------------------------------------------------------------------*/

  swRunningLag = swOrigLagIn;
  siSampsSoFar = 0;
  while (siSampsSoFar < S_LEN)
  {

    /* Get integer lag and remainder.  These are used in addressing */
    /* the LTP state and the interpolating filter, respectively     */
    /*--------------------------------------------------------------*/

    get_ipjj(swRunningLag, &swIntLag, &swRemain);


    /* Get the number of samples to look up in this pass */
    /*---------------------------------------------------*/

    if (sub(swIntLag, S_LEN) < 0)
      siSampsThisPass = swIntLag - siSampsSoFar;
    else
      siSampsThisPass = S_LEN - siSampsSoFar;

    /* Look up samples by interpolating (fractional lag), or copying */
    /* (integer lag).                                                */
    /*---------------------------------------------------------------*/

    if (swRemain == 0)
    {

      /* Integer lag: copy samples from history */
      /*----------------------------------------*/

      for (i = siSampsSoFar; i < siSampsSoFar + siSampsThisPass; i++)
        pswLTPState[i] = pswLTPState[i - swIntLag];
    }
    else
    {

      /* Fractional lag: interpolate to get samples */
      /*--------------------------------------------*/

      for (i = siSampsSoFar; i < siSampsSoFar + siSampsThisPass; i++)
      {

        /* first tap with rounding offset */
        /*--------------------------------*/
        L_Temp = L_mac((long) 32768,
                       pswLTPState[i - swIntLag - P_INT_MACS / 2],
                       ppsrPVecIntFilt[0][swRemain]);

        for (j = 1; j < P_INT_MACS - 1; j++)
        {

          L_Temp = L_mac(L_Temp,
                         pswLTPState[i - swIntLag - P_INT_MACS / 2 + j],
                         ppsrPVecIntFilt[j][swRemain]);

        }

        pswLTPState[i] = extract_h(L_mac(L_Temp,
                             pswLTPState[i - swIntLag + P_INT_MACS / 2 - 1],
                                ppsrPVecIntFilt[P_INT_MACS - 1][swRemain]));
      }
    }

    /* Done with this pass: update loop controls */
    /*-------------------------------------------*/

    siSampsSoFar += siSampsThisPass;
    swRunningLag = add(swRunningLag, swOrigLagIn);
  }
}

/***************************************************************************
 *
 *    FUNCTION NAME: g_corr1 (no scaling)
 *
 *    PURPOSE:
 *
 *     Calculates energy in subframe vector.  Differs from g_corr1s,
 *     in that there the estimate of the maximum possible
 *     energy is < 1.0
 *
 *
 *    INPUT:
 *
 *       pswIn[0:39]
 *                     A subframe vector.
 *
 *
 *    OUTPUT:
 *
 *       *pL_out
 *                     A Longword containing the normalized energy
 *                     in the input vector.
 *
 *    RETURN:
 *
 *       swOut
 *                     Number of right shifts which the accumulator was
 *                     shifted to normalize it.  Negative number implies
 *                     a left shift, and therefore an energy larger than
 *                     1.0.
 *
 *    REFERENCES: Sub_Clause 4.1.10.2 and 4.2.1 of GSM Recomendation 06.20
 *
 *    KEYWORDS: energy, autocorrelation, correlation, g_corr1
 *
 *
 *************************************************************************/

Shortword g_corr1(Shortword *pswIn, Longword *pL_out)
{


/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_sum;
  Shortword swEngyLShft;
  int    i;


/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/


  /* Calculate energy in subframe vector (40 samples) */
  /*--------------------------------------------------*/

  L_sum = L_mult(pswIn[0], pswIn[0]);
  for (i = 1; i < S_LEN; i++)
  {
    L_sum = L_mac(L_sum, pswIn[i], pswIn[i]);
  }



  if (L_sum != 0)
  {

    /* Normalize the energy in the output Longword */
    /*---------------------------------------------*/

    swEngyLShft = norm_l(L_sum);
    *pL_out = L_shl(L_sum, swEngyLShft);        /* normalize output
                                                 * Longword */
  }
  else
  {

    /* Special case: energy is zero */
    /*------------------------------*/

    *pL_out = L_sum;
    swEngyLShft = 0;
  }

  return (swEngyLShft);
}

/***************************************************************************
 *
 *    FUNCTION NAME: g_corr1s (g_corr1 with scaling)
 *
 *    PURPOSE:
 *
 *     Calculates energy in subframe vector.  Differs from g_corr1,
 *     in that there is an estimate of the maximum possible energy in the
 *     vector.
 *
 *    INPUT:
 *
 *       pswIn[0:39]
 *                     A subframe vector.
 *
 *       swEngyRShft
 *
 *                     Number of right shifts to apply to the vectors energy
 *                     to ensure that it remains less than 1.0
 *                     (swEngyRShft is always positive or zero)
 *
 *    OUTPUT:
 *
 *       *pL_out
 *                     A Longword containing the normalized energy
 *                     in the input vector.
 *
 *    RETURN:
 *
 *       swOut
 *                     Number of right shifts which the accumulator was
 *                     shifted to normalize it.  Negative number implies
 *                     a left shift, and therefore an energy larger than
 *                     1.0.
 *
 *    REFERENCES: Sub-Clause 4.1.8, 4.2.1, 4.2.2, and 4.2.4
 *                of GSM Recomendation 06.20
 *
 *    keywords: energy, autocorrelation, correlation, g_corr1
 *
 *
 *************************************************************************/

Shortword g_corr1s(Shortword pswIn[], Shortword swEngyRShft,
                          Longword *pL_out)
{


/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_sum;
  Shortword swTemp,
         swEngyLShft;
  Shortword swInputRShft;

  int    i;


/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/


  /* Calculate energy in subframe vector (40 samples) */
  /*--------------------------------------------------*/

  if (sub(swEngyRShft, 1) <= 0)
  {

    /* use the energy shift factor, although it is an odd shift count */
    /*----------------------------------------------------------------*/

    swTemp = shr(pswIn[0], swEngyRShft);
    L_sum = L_mult(pswIn[0], swTemp);
    for (i = 1; i < S_LEN; i++)
    {
      swTemp = shr(pswIn[i], swEngyRShft);
      L_sum = L_mac(L_sum, pswIn[i], swTemp);
    }

  }
  else
  {

    /* convert energy shift factor to an input shift factor */
    /*------------------------------------------------------*/

    swInputRShft = shift_r(swEngyRShft, -1);
    swEngyRShft = shl(swInputRShft, 1);

    swTemp = shr(pswIn[0], swInputRShft);
    L_sum = L_mult(swTemp, swTemp);
    for (i = 1; i < S_LEN; i++)
    {
      swTemp = shr(pswIn[i], swInputRShft);
      L_sum = L_mac(L_sum, swTemp, swTemp);
    }
  }

  if (L_sum != 0)
  {

    /* Normalize the energy in the output Longword */
    /*---------------------------------------------*/

    swTemp = norm_l(L_sum);
    *pL_out = L_shl(L_sum, swTemp);    /* normalize output Longword */
    swEngyLShft = sub(swTemp, swEngyRShft);
  }
  else
  {

    /* Special case: energy is zero */
    /*------------------------------*/

    *pL_out = L_sum;
    swEngyLShft = 0;
  }

  return (swEngyLShft);
}

/***************************************************************************
 *
 *   FUNCTION NAME: getSfrmLpc
 *
 *   PURPOSE:
 *
 *     Given frame information from past and present frame, interpolate
 *     (or copy) the frame based LPC coefficients into subframe
 *     lpc coeffs, i.e. the ones which will be used by the subframe
 *     as opposed to those coded and transmitted.
 *
 *   INPUTS:
 *
 *     siSoftInterpolation
 *
 *                     interpolate 1/0, a coded parameter.
 *
 *     swPrevR0,swNewR0
 *
 *                     Rq0 for the last frame and for this frame.
 *                     These are the decoded values, not the codewords.
 *
 *     Previous lpc coefficients from the previous frame:
 *       in all filters below array[0] is the t=-1 element array[9]
 *       t=-10 element.
 *
 *     pswPrevFrmKs[0:9]
 *
 *                     decoded version of the rc's tx'd last frame
 *
 *     pswPrevFrmAs[0:9]
 *
 *                     the above K's converted to A's.  i.e. direct
 *                     form coefficients.
 *
 *     pswPrevFrmPFNum[0:9], pswPrevFrmPFDenom[0:9]
 *
 *                     numerator and denominator coefficients used in the
 *                     postfilter
 *
 *     Current lpc coefficients from the current frame:
 *
 *     pswNewFrmKs[0:9], pswNewFrmAs[0:9],
 *     pswNewFrmPFNum[0:9], pswNewFrmPFDenom[0:9] same as above.
 *
 *   OUTPUTS:
 *
 *     psnsSqrtRs[0:3]
 *
 *                      a normalized number (struct NormSw)
 *                      containing an estimate of RS for each subframe.
 *                      (number and a shift)
 *
 *     ppswSynthAs[0:3][0:9]
 *
 *                      filter coefficients used by the synthesis filter.
 *
 *     ppswPFNumAs[0:3][0:9]
 *
 *                      filter coefficients used by the postfilters
 *                      numerator.
 *
 *     ppswPFDenomAs[0:3][0:9]
 *
 *                      filter coefficients used by postfilters denominator.
 *
 *   RETURN VALUE:
 *
 *     None
 *
 *   DESCRIPTION:
 *
 *     For interpolated subframes, the direct form coefficients
 *     are converted to reflection coeffiecients to check for
 *     filter stability. If unstable, the uninterpolated coef.
 *     are used for that subframe.
 *
 *     Interpolation is described in section 4.1.6, "Soft Interpolation
 *     of the Spectral Parameters"
 *
 *    REFERENCES: Sub_clause 4.2.1 of GSM Recomendation 06.20
 *
 *   KEYWORDS: soft interpolation, int_lpc, interpolate, atorc,res_eng,i_mov
 *
 *************************************************************************/

void   getSfrmLpc(short int siSoftInterpolation,
                         Shortword swPrevR0, Shortword swNewR0,
          /* last frm */ Shortword pswPrevFrmKs[], Shortword pswPrevFrmAs[],
                         Shortword pswPrevFrmPFNum[],
                         Shortword pswPrevFrmPFDenom[],

            /* this frm */ Shortword pswNewFrmKs[], Shortword pswNewFrmAs[],
                         Shortword pswNewFrmPFNum[],
                         Shortword pswNewFrmPFDenom[],

                   /* output */ struct NormSw *psnsSqrtRs,
                         Shortword *ppswSynthAs[], Shortword *ppswPFNumAs[],
                         Shortword *ppswPFDenomAs[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                           Local Static Variables                        |
 |_________________________________________________________________________|
*/


/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  short int siSfrm,
         siStable,
         i;

  Longword L_Temp1,
         L_Temp2;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  if (siSoftInterpolation)
  {
    /* yes, interpolating */
    /* ------------------ */

    siSfrm = 0;

    siStable = interpolateCheck(pswPrevFrmKs, pswPrevFrmAs,
                                pswPrevFrmAs, pswNewFrmAs,
                                psrOldCont[siSfrm], psrNewCont[siSfrm],
                                swPrevR0,
                                &psnsSqrtRs[siSfrm],
                                ppswSynthAs[siSfrm]);
    if (siStable)
    {

      /* interpolate between direct form coefficient sets */
      /* for both numerator and denominator coefficients  */
      /* assume output will be stable                     */
      /* ------------------------------------------------ */

      for (i = 0; i < NP; i++)
      {
        L_Temp1 = L_mult(pswNewFrmPFNum[i], psrNewCont[siSfrm]);
        ppswPFNumAs[siSfrm][i] = mac_r(L_Temp1, pswPrevFrmPFNum[i],
                                       psrOldCont[siSfrm]);
        L_Temp2 = L_mult(pswNewFrmPFDenom[i], psrNewCont[siSfrm]);
        ppswPFDenomAs[siSfrm][i] = mac_r(L_Temp2, pswPrevFrmPFDenom[i],
                                         psrOldCont[siSfrm]);
      }
    }
    else
    {
      /* this subframe is unstable */
      /* ------------------------- */
      for (i = 0; i < NP; i++)
      {
        ppswPFNumAs[siSfrm][i] = pswPrevFrmPFNum[i];
        ppswPFDenomAs[siSfrm][i] = pswPrevFrmPFDenom[i];
      }
    }
    for (siSfrm = 1; siSfrm < N_SUB - 1; siSfrm++)
    {

      siStable = interpolateCheck(pswNewFrmKs, pswNewFrmAs,
                                  pswPrevFrmAs, pswNewFrmAs,
                                  psrOldCont[siSfrm], psrNewCont[siSfrm],
                                  swNewR0,
                                  &psnsSqrtRs[siSfrm],
                                  ppswSynthAs[siSfrm]);
      if (siStable)
      {

        /* interpolate between direct form coefficient sets */
        /* for both numerator and denominator coefficients  */
        /* assume output will be stable                     */
        /* ------------------------------------------------ */

        for (i = 0; i < NP; i++)
        {
          L_Temp1 = L_mult(pswNewFrmPFNum[i], psrNewCont[siSfrm]);
          ppswPFNumAs[siSfrm][i] = mac_r(L_Temp1, pswPrevFrmPFNum[i],
                                         psrOldCont[siSfrm]);
          L_Temp2 = L_mult(pswNewFrmPFDenom[i], psrNewCont[siSfrm]);
          ppswPFDenomAs[siSfrm][i] = mac_r(L_Temp2, pswPrevFrmPFDenom[i],
                                           psrOldCont[siSfrm]);
        }
      }
      else
      {
        /* this subframe has unstable filter coeffs, would like to
         * interpolate but can not  */
        /* -------------------------------------- */
        for (i = 0; i < NP; i++)
        {
          ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i];
          ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i];
        }
      }
    }
    /* the last subframe never interpolate */
    /* ----------------------------------- */
    siSfrm = 3;
    for (i = 0; i < NP; i++)
    {
      ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i];
      ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i];
      ppswSynthAs[siSfrm][i] = pswNewFrmAs[i];
    }

    res_eng(pswNewFrmKs, swNewR0, &psnsSqrtRs[siSfrm]);

  }
  /* SoftInterpolation == 0  - no interpolation */
  /* ------------------------------------------ */
  else
  {
    siSfrm = 0;
    for (i = 0; i < NP; i++)
    {
      ppswPFNumAs[siSfrm][i] = pswPrevFrmPFNum[i];
      ppswPFDenomAs[siSfrm][i] = pswPrevFrmPFDenom[i];
      ppswSynthAs[siSfrm][i] = pswPrevFrmAs[i];
    }

    res_eng(pswPrevFrmKs, swPrevR0, &psnsSqrtRs[siSfrm]);

    /* for subframe 1 and all subsequent sfrms, use result from new frm */
    /* ---------------------------------------------------------------- */


    res_eng(pswNewFrmKs, swNewR0, &psnsSqrtRs[1]);

    for (siSfrm = 1; siSfrm < N_SUB; siSfrm++)
    {


      psnsSqrtRs[siSfrm].man = psnsSqrtRs[1].man;
      psnsSqrtRs[siSfrm].sh = psnsSqrtRs[1].sh;

      for (i = 0; i < NP; i++)
      {
        ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i];
        ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i];
        ppswSynthAs[siSfrm][i] = pswNewFrmAs[i];
      }
    }
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: get_ipjj
 *
 *   PURPOSE:
 *
 *     This subroutine calculates IP, the single-resolution lag rounded
 *     down to the nearest integer, and JJ, the remainder when the
 *     extended resolution lag is divided by the oversampling factor
 *
 *   INPUTS:
 *
 *     swLagIn
 *                     extended resolution lag as an integer, i.e.
 *                     fractional lag x oversampling factor
 *
 *   OUTPUTS:
 *
 *     *pswIp
 *                     fractional lag rounded down to nearest integer, IP
 *
 *     *pswJj
 *                     the remainder JJ
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     ip = integer[lag/OS_FCTR]
 *     jj = integer_round[((lag/OS_FCTR)-ip)*(OS_FCTR)]
 *          if the rounding caused an 'overflow'
 *            set remainder jj to 0 and add 'carry' to ip
 *
 *     This routine is involved in the mechanics of fractional and
 *     integer LTP searchs.  The LTP is described in section 5.
 *
 *   REFERENCES: Sub-clause 4.1.8 and 4.2.2 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lag, fractional, remainder, ip, jj, get_ipjj
 *
 *************************************************************************/

void   get_ipjj(Shortword swLagIn,
                       Shortword *pswIp, Shortword *pswJj)
{

/*_________________________________________________________________________
 |                                                                         |
 |                              Local Constants                            |
 |_________________________________________________________________________|
*/

#define  OS_FCTR_INV  (Shortword)0x1555/* SW_MAX/OS_FCTR */

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Temp;

  Shortword swTemp,
         swTempIp,
         swTempJj;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* calculate ip */
  /* ------------ */

  L_Temp = L_mult(OS_FCTR_INV, swLagIn);        /* lag/OS_FCTR */
  swTempIp = extract_h(L_Temp);

  /* calculate jj */
  /* ------------ */

  swTemp = extract_l(L_Temp);          /* loose ip */
  swTemp = shr(swTemp, 1);             /* isolate jj fraction */
  swTemp = swTemp & SW_MAX;
  L_Temp = L_mult(swTemp, OS_FCTR);    /* ((lag/OS_FCTR)-ip))*(OS_FCTR) */
  swTemp = round(L_Temp);              /* round and pick-off jj */
  if (sub(swTemp, OS_FCTR) == 0)
  {                                    /* if 'overflow ' */
    swTempJj = 0;                      /* set remainder,jj to 0 */
    swTempIp = add(swTempIp, 1);       /* 'carry' overflow into ip */
  }
  else
  {
    swTempJj = swTemp;                 /* read-off remainder,jj */
  }

  /* return ip and jj */
  /* ---------------- */

  *pswIp = swTempIp;
  *pswJj = swTempJj;
}

/***************************************************************************
 *
 *   FUNCTION NAME: interpolateCheck
 *
 *   PURPOSE:
 *
 *     Interpolates between direct form coefficient sets.
 *     Before releasing the interpolated coefficients, they are checked.
 *     If unstable, the "old" parameters are used.
 *
 *   INPUTS:
 *
 *     pswRefKs[0:9]
 *                     decoded version of the rc's tx'd last frame
 *
 *     pswRefCoefsA[0:9]
 *                     above K's converted to direct form coefficients
 *
 *     pswOldCoefsA[0:9]
 *                     array of old Coefseters
 *
 *     pswNewCoefsA[0:9]
 *                     array of new Coefseters
 *
 *     swOldPer
 *                     amount old coefs supply to the output
 *
 *     swNewPer
 *                     amount new coefs supply to the output
 *
 *     ASHIFT
 *                     shift for reflection coef. conversion
 *
 *     swRq
 *                     quantized energy to use for subframe
 * *
 *   OUTPUTS:
 *
 *     psnsSqrtRsOut
 *                     output pointer to sqrt(RS) normalized
 *
 *     pswCoefOutA[0:9]
 *                     output coefficients
 *
 *   RETURN VALUE:
 *
 *     siInterp_flg
 *                     temporary subframe interpolation flag
 *                     0 - coef. interpolated, 1 -coef. not interpolated
 *
 *   DESCRIPTION:
 *
 *     For interpolated subframes, the direct form coefficients
 *     are converted to reflection coefficients to check for
 *     filter stability. If unstable, the uninterpolated coef.
 *     are used for that subframe.  Section 4.1.6 describes
 *     interpolation.
 *
 *   REFERENCES: Sub-clause 4.1.6 and 4.2.3 of GSM Recomendation 06.20
 *
 *   KEYWORDS: soft interpolation, int_lpc, interpolate, atorc,res_eng,i_mov
 *
 *************************************************************************/

short int interpolateCheck(Shortword pswRefKs[],
                                  Shortword pswRefCoefsA[],
                         Shortword pswOldCoefsA[], Shortword pswNewCoefsA[],
                                  Shortword swOldPer, Shortword swNewPer,
                                  Shortword swRq,
                                  struct NormSw *psnsSqrtRsOut,
                                  Shortword pswCoefOutA[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Shortword pswRcTemp[NP];

  Longword L_Temp;

  short int siInterp_flg,
         i;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Interpolation loop, NP is order of LPC filter */
  /* --------------------------------------------- */

  for (i = 0; i < NP; i++)
  {
    L_Temp = L_mult(pswNewCoefsA[i], swNewPer);
    pswCoefOutA[i] = mac_r(L_Temp, pswOldCoefsA[i], swOldPer);
  }

  /* Convert to reflection coefficients and check stability */
  /* ------------------------------------------------------ */

  if (aToRc(ASHIFT, pswCoefOutA, pswRcTemp) != 0)
  {

    /* Unstable, use uninterpolated parameters and compute RS update the
     * state with the frame data closest to this subfrm */
    /* --------------------------------------------------------- */

    res_eng(pswRefKs, swRq, psnsSqrtRsOut);

    for (i = 0; i < NP; i++)
    {
      pswCoefOutA[i] = pswRefCoefsA[i];
    }
    siInterp_flg = 0;
  }
  else
  {

    /* Stable, compute RS */
    /* ------------------ */
    res_eng(pswRcTemp, swRq, psnsSqrtRsOut);

    /* Set temporary subframe interpolation flag */
    /* ----------------------------------------- */
    siInterp_flg = 1;
  }

  /* Return subframe interpolation flag */
  /* ---------------------------------- */
  return (siInterp_flg);
}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcFir
 *
 *   PURPOSE:
 *
 *     The purpose of this function is to perform direct form fir filtering
 *     assuming a NP order filter and given state, coefficients, and input.
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter
 *
 *     pswInput[0:S_LEN-1]
 *
 *                     input array of points to be filtered.
 *                     pswInput[0] is the oldest point (first to be filtered)
 *                     pswInput[siLen-1] is the last point filtered (newest)
 *
 *     pswCoef[0:NP-1]
 *
 *                     array of direct form coefficients
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *     pswState[0:NP-1]
 *
 *                     array of the filter state following form of pswCoef
 *                     pswState[0] = state of filter for delay n = -1
 *                     pswState[NP-1] = state of filter for delay n = -NP
 *
 *   OUTPUTS:
 *
 *     pswState[0:NP-1]
 *
 *                     updated filter state, ready to filter
 *                     pswInput[siLen], i.e. the next point
 *
 *     pswFiltOut[0:S_LEN-1]
 *
 *                     the filtered output
 *                     same format as pswInput, pswFiltOut[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] += state[j]*coef[j] (state is taken from either input
 *                                       state[] or input in[] arrays)
 *        rescale(out[i])
 *        out[i] += in[i]
 *     update final state array using in[]
 *
 *   REFERENCES: Sub-clause 4.1.7 and 4.2.4 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, fir, lpcFir, inversefilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set, i_dir_mod
 *
 *************************************************************************/

void   lpcFir(Shortword pswInput[], Shortword pswCoef[],
                     Shortword pswState[], Shortword pswFiltOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* filter 1st sample */
  /* ----------------- */

  /* sum past state outputs */
  /* ---------------------- */
  /* 0th coef, with rounding */
  L_Sum = L_mac(LPC_ROUND, pswState[0], pswCoef[0]);

  for (siStage = 1; siStage < NP; siStage++)
  {                                    /* remaining coefs */
    L_Sum = L_mac(L_Sum, pswState[siStage], pswCoef[siStage]);
  }

  /* add input to partial output */
  /* --------------------------- */

  L_Sum = L_shl(L_Sum, ASHIFT);
  L_Sum = L_msu(L_Sum, pswInput[0], 0x8000);

  /* save 1st output sample */
  /* ---------------------- */

  pswFiltOut[0] = extract_h(L_Sum);

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_mac(LPC_ROUND, pswInput[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_mac(L_Sum, pswInput[siSmp - siStage - 1], pswCoef[siStage]);
    }

    /* sum past states, if any */
    /* ----------------------- */

    for (siStage = siSmp; siStage < NP; siStage++)
    {
      L_Sum = L_mac(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]);
    }

    /* add input to partial output */
    /* --------------------------- */

    L_Sum = L_shl(L_Sum, ASHIFT);
    L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);

    /* save current output sample */
    /* -------------------------- */

    pswFiltOut[siSmp] = extract_h(L_Sum);
  }

  /* save final state */
  /* ---------------- */

  for (siStage = 0; siStage < NP; siStage++)
  {
    pswState[siStage] = pswInput[S_LEN - siStage - 1];
  }

}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcIir
 *
 *   PURPOSE:
 *
 *     The purpose of this function is to perform direct form IIR filtering
 *     assuming a NP order filter and given state, coefficients, and input
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter
 *
 *     pswInput[0:S_LEN-1]
 *
 *                     input array of points to be filtered
 *                     pswInput[0] is the oldest point (first to be filtered)
 *                     pswInput[siLen-1] is the last point filtered (newest)
 *
 *     pswCoef[0:NP-1]
 *                     array of direct form coefficients
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *     pswState[0:NP-1]
 *
 *                     array of the filter state following form of pswCoef
 *                     pswState[0] = state of filter for delay n = -1
 *                     pswState[NP-1] = state of filter for delay n = -NP
 *
 *   OUTPUTS:
 *
 *     pswState[0:NP-1]
 *
 *                     updated filter state, ready to filter
 *                     pswInput[siLen], i.e. the next point
 *
 *     pswFiltOut[0:S_LEN-1]
 *
 *                     the filtered output
 *                     same format as pswInput, pswFiltOut[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] -= state[j]*coef[j] (state is taken from either input
 *                                       state[] or prior out[] arrays)
 *        rescale(out[i])
 *        out[i] += in[i]
 *     update final state array using out[]
 *
 *   REFERENCES: Sub-clause 4.1.7 and 4.2.4 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
 *
 *************************************************************************/

void   lpcIir(Shortword pswInput[], Shortword pswCoef[],
                     Shortword pswState[], Shortword pswFiltOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* filter 1st sample */
  /* ----------------- */

  /* sum past state outputs */
  /* ---------------------- */
  /* 0th coef, with rounding */
  L_Sum = L_msu(LPC_ROUND, pswState[0], pswCoef[0]);

  for (siStage = 1; siStage < NP; siStage++)
  {                                    /* remaining coefs */
    L_Sum = L_msu(L_Sum, pswState[siStage], pswCoef[siStage]);
  }

  /* add input to partial output */
  /* --------------------------- */

  L_Sum = L_shl(L_Sum, ASHIFT);
  L_Sum = L_msu(L_Sum, pswInput[0], 0x8000);

  /* save 1st output sample */
  /* ---------------------- */

  pswFiltOut[0] = extract_h(L_Sum);

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
                    pswCoef[siStage]);
    }

    /* sum past states, if any */
    /* ----------------------- */

    for (siStage = siSmp; siStage < NP; siStage++)
    {
      L_Sum = L_msu(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]);
    }

    /* add input to partial output */
    /* --------------------------- */

    L_Sum = L_shl(L_Sum, ASHIFT);
    L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);

    /* save current output sample */
    /* -------------------------- */

    pswFiltOut[siSmp] = extract_h(L_Sum);
  }

  /* save final state */
  /* ---------------- */

  for (siStage = 0; siStage < NP; siStage++)
  {
    pswState[siStage] = pswFiltOut[S_LEN - siStage - 1];
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcIrZsIir
 *
 *   PURPOSE:
 *
 *     The purpose of this function is to calculate the impulse response
 *     via direct form IIR filtering with zero state assuming a NP order
 *     filter and given coefficients
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter
 *
 *     pswCoef[0:NP-1]
 *                     array of direct form coefficients
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *   OUTPUTS:
 *
 *     pswFiltOut[0:S_LEN-1]
 *
 *                     the filtered output
 *                     same format as pswInput, pswFiltOut[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     This routine is called by getNWCoefs().
 *
 *     Because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] -= state[j]*coef[j] (state taken from prior output[])
 *        rescale(out[i])
 *
 *   REFERENCES: Sub-clause 4.1.8 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
 *
 *************************************************************************/

void   lpcIrZsIir(Shortword pswCoef[], Shortword pswFiltOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* output 1st sample */
  /* ----------------- */

  pswFiltOut[0] = 0x0400;

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
                    pswCoef[siStage]);
    }

    /* scale output */
    /* ------------ */

    L_Sum = L_shl(L_Sum, ASHIFT);

    /* save current output sample */
    /* -------------------------- */

    pswFiltOut[siSmp] = extract_h(L_Sum);
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcZiIir
 *
 *   PURPOSE:
 *     The purpose of this function is to perform direct form iir filtering
 *     with zero input assuming a NP order filter, and given state and
 *     coefficients
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter MUST be <= MAX_ZIS
 *
 *     pswCoef[0:NP-1]
 *
 *                     array of direct form coefficients.
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *     pswState[0:NP-1]
 *
 *                     array of the filter state following form of pswCoef
 *                     pswState[0] = state of filter for delay n = -1
 *                     pswState[NP-1] = state of filter for delay n = -NP
 *
 *   OUTPUTS:
 *
 *     pswFiltOut[0:S_LEN-1]
 *
 *                     the filtered output
 *                     same format as pswIn, pswFiltOut[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     The routine is called from sfrmAnalysis, and is used to let the
 *     LPC filters ring out.
 *
 *     because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] -= state[j]*coef[j] (state is taken from either input
 *                                       state[] or prior output[] arrays)
 *        rescale(out[i])
 *
 *   REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
 *
 *************************************************************************/

void   lpcZiIir(Shortword pswCoef[], Shortword pswState[],
                       Shortword pswFiltOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* filter 1st sample */
  /* ----------------- */

  /* sum past state outputs */
  /* ---------------------- */
  /* 0th coef, with rounding */
  L_Sum = L_msu(LPC_ROUND, pswState[0], pswCoef[0]);

  for (siStage = 1; siStage < NP; siStage++)
  {                                    /* remaining coefs */
    L_Sum = L_msu(L_Sum, pswState[siStage], pswCoef[siStage]);
  }

  /* scale output */
  /* ------------ */

  L_Sum = L_shl(L_Sum, ASHIFT);

  /* save 1st output sample */
  /* ---------------------- */

  pswFiltOut[0] = extract_h(L_Sum);

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
                    pswCoef[siStage]);
    }

    /* sum past states, if any */
    /* ----------------------- */

    for (siStage = siSmp; siStage < NP; siStage++)
    {
      L_Sum = L_msu(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]);
    }

    /* scale output */
    /* ------------ */

    L_Sum = L_shl(L_Sum, ASHIFT);

    /* save current output sample */
    /* -------------------------- */

    pswFiltOut[siSmp] = extract_h(L_Sum);
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcZsFir
 *
 *   PURPOSE:
 *     The purpose of this function is to perform direct form fir filtering
 *     with zero state, assuming a NP order filter and given coefficients
 *     and non-zero input.
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter
 *
 *     pswInput[0:S_LEN-1]
 *
 *                     input array of points to be filtered.
 *                     pswInput[0] is the oldest point (first to be filtered)
 *                     pswInput[siLen-1] is the last point filtered (newest)
 *
 *     pswCoef[0:NP-1]
 *
 *                     array of direct form coefficients
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *   OUTPUTS:
 *
 *     pswFiltOut[0:S_LEN-1]
 *
 *                     the filtered output
 *                     same format as pswInput, pswFiltOut[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     This routine is used in getNWCoefs().  See section 4.1.7.
 *
 *     because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] += state[j]*coef[j] (state taken from in[])
 *        rescale(out[i])
 *        out[i] += in[i]
 *
 *   REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, fir, lpcFir, inversefilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set, i_dir_mod
 *
 *************************************************************************/

void   lpcZsFir(Shortword pswInput[], Shortword pswCoef[],
                       Shortword pswFiltOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* output 1st sample */
  /* ----------------- */

  pswFiltOut[0] = pswInput[0];

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_mac(LPC_ROUND, pswInput[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_mac(L_Sum, pswInput[siSmp - siStage - 1],
                    pswCoef[siStage]);
    }

    /* add input to partial output */
    /* --------------------------- */

    L_Sum = L_shl(L_Sum, ASHIFT);
    L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);

    /* save current output sample */
    /* -------------------------- */

    pswFiltOut[siSmp] = extract_h(L_Sum);
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcZsIir
 *
 *   PURPOSE:
 *
 *     The purpose of this function is to perform direct form IIR filtering
 *     with zero state, assuming a NP order filter and given coefficients
 *     and non-zero input.
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter
 *
 *     pswInput[0:S_LEN-1]
 *
 *                     input array of points to be filtered
 *                     pswInput[0] is the oldest point (first to be filtered)
 *                     pswInput[siLen-1] is the last point filtered (newest)
 *
 *     pswCoef[0:NP-1]
 *                     array of direct form coefficients
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *   OUTPUTS:
 *
 *     pswFiltOut[0:S_LEN-1]
 *
 *                     the filtered output
 *                     same format as pswInput, pswFiltOut[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     This routine is used in the subframe analysis process.  It is
 *     called by sfrmAnalysis() and fnClosedLoop().  It is this function
 *     which performs the weighting of the excitation vectors.
 *
 *     because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] -= state[j]*coef[j] (state taken from prior out[])
 *        rescale(out[i])
 *        out[i] += in[i]
 *
 *   REFERENCES: Sub-clause 4.1.8.5 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
 *
 *************************************************************************/

void   lpcZsIir(Shortword pswInput[], Shortword pswCoef[],
                       Shortword pswFiltOut[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* output 1st sample */
  /* ----------------- */

  pswFiltOut[0] = pswInput[0];

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
                    pswCoef[siStage]);
    }

    /* add input to partial output */
    /* --------------------------- */

    L_Sum = L_shl(L_Sum, ASHIFT);
    L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);

    /* save current output sample */
    /* -------------------------- */

    pswFiltOut[siSmp] = extract_h(L_Sum);
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: lpcZsIirP
 *
 *   PURPOSE:
 *
 *     The purpose of this function is to perform direct form iir filtering
 *     with zero state, assuming a NP order filter and given coefficients
 *     and input
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the lpc filter
 *
 *     S_LEN
 *                     number of samples to filter
 *
 *     pswCommonIO[0:S_LEN-1]
 *
 *                     input array of points to be filtered
 *                     pswCommonIO[0] is oldest point (first to be filtered)
 *                     pswCommonIO[siLen-1] is last point filtered (newest)
 *
 *     pswCoef[0:NP-1]
 *                     array of direct form coefficients
 *                     pswCoef[0] = coeff for delay n = -1
 *                     pswCoef[NP-1] = coeff for delay n = -NP
 *
 *     ASHIFT
 *                     number of shifts input A's have been shifted down by
 *
 *     LPC_ROUND
 *                     rounding constant
 *
 *   OUTPUTS:
 *
 *     pswCommonIO[0:S_LEN-1]
 *
 *                     the filtered output
 *                     pswCommonIO[0] is oldest point
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     This function is called by geNWCoefs().  See section 4.1.7.
 *
 *     because of the default sign of the coefficients the
 *     formula for the filter is :
 *     i=0, i < S_LEN
 *        out[i] = rounded(state[i]*coef[0])
 *        j=1, j < NP
 *           out[i] += state[j]*coef[j] (state taken from prior out[])
 *        rescale(out[i])
 *        out[i] += in[i]
 *
 *   REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20
 *
 *   KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
 *   KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
 *
 *************************************************************************/

void   lpcZsIirP(Shortword pswCommonIO[], Shortword pswCoef[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Sum;
  short int siStage,
         siSmp;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* filter remaining samples */
  /* ------------------------ */

  for (siSmp = 1; siSmp < S_LEN; siSmp++)
  {

    /* sum past outputs */
    /* ---------------- */
    /* 0th coef, with rounding */
    L_Sum = L_mac(LPC_ROUND, pswCommonIO[siSmp - 1], pswCoef[0]);
    /* remaining coefs */
    for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
    {
      L_Sum = L_mac(L_Sum, pswCommonIO[siSmp - siStage - 1],
                    pswCoef[siStage]);
    }

    /* add input to partial output */
    /* --------------------------- */

    L_Sum = L_shl(L_Sum, ASHIFT);
    L_Sum = L_msu(L_Sum, pswCommonIO[siSmp], 0x8000);

    /* save current output sample */
    /* -------------------------- */

    pswCommonIO[siSmp] = extract_h(L_Sum);
  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: r0BasedEnergyShft
 *
 *   PURPOSE:
 *
 *     Given an R0 voicing level, find the number of shifts to be
 *     performed on the energy to ensure that the subframe energy does
 *     not overflow.  example if energy can maximally take the value
 *     4.0, then 2 shifts are required.
 *
 *   INPUTS:
 *
 *     swR0Index
 *                     R0 codeword (0-0x1f)
 *
 *   OUTPUTS:
 *
 *     none
 *
 *   RETURN VALUE:
 *
 *     swShiftDownSignal
 *
 *                    number of right shifts to apply to energy (0..6)
 *
 *   DESCRIPTION:
 *
 *     Based on the R0, the average frame energy, we can get an
 *     upper bound on the energy any one subframe can take on.
 *     Using this upper bound we can calculate what right shift is
 *     needed to ensure an unsaturated output out of a subframe
 *     energy calculation (g_corr).
 *
 *   REFERENCES: Sub-clause 4.1.9 and 4.2.1 of GSM Recomendation 06.20
 *
 *   KEYWORDS: spectral postfilter
 *
 *************************************************************************/

Shortword r0BasedEnergyShft(Shortword swR0Index)
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Shortword swShiftDownSignal;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  if (sub(swR0Index, 26) <= 0)
  {
    if (sub(swR0Index, 23) <= 0)
    {
      if (sub(swR0Index, 21) <= 0)
        swShiftDownSignal = 0;         /* r0  [0,  21] */
      else
        swShiftDownSignal = 1;         /* r0  [22, 23] */
    }
    else
    {
      if (sub(swR0Index, 24) <= 0)
        swShiftDownSignal = 2;         /* r0  [23, 24] */
      else
        swShiftDownSignal = 3;         /* r0  [24, 26] */
    }
  }
  else
  {                                    /* r0 index > 26 */
    if (sub(swR0Index, 28) <= 0)
    {
      swShiftDownSignal = 4;           /* r0  [26, 28] */
    }
    else
    {
      if (sub(swR0Index, 29) <= 0)
        swShiftDownSignal = 5;         /* r0  [28, 29] */
      else
        swShiftDownSignal = 6;         /* r0  [29, 31] */
    }
  }
  if (sub(swR0Index, 18) > 0)
    swShiftDownSignal = add(swShiftDownSignal, 2);

  return (swShiftDownSignal);
}

/***************************************************************************
 *
 *   FUNCTION NAME: rcToADp
 *
 *   PURPOSE:
 *
 *     This subroutine computes a vector of direct form LPC filter
 *     coefficients, given an input vector of reflection coefficients.
 *     Double precision is used internally, but 16 bit direct form
 *     filter coefficients are returned.
 *
 *   INPUTS:
 *
 *     NP
 *                     order of the LPC filter (global constant)
 *
 *     swAscale
 *                     The multiplier which scales down the direct form
 *                     filter coefficients.
 *
 *     pswRc[0:NP-1]
 *                     The input vector of reflection coefficients.
 *
 *   OUTPUTS:
 *
 *     pswA[0:NP-1]
 *                     Array containing the scaled down direct form LPC
 *                     filter coefficients.
 *
 *   RETURN VALUE:
 *
 *     siLimit
 *                     1 if limiting occured in computation, 0 otherwise.
 *
 *   DESCRIPTION:
 *
 *     This function performs the conversion from reflection coefficients
 *     to direct form LPC filter coefficients. The direct form coefficients
 *     are scaled by multiplication by swAscale. NP, the filter order is 10.
 *     The a's and rc's each have NP elements in them. Double precision
 *     calculations are used internally.
 *
 *        The equations are:
 *        for i = 0 to NP-1{
 *
 *          a(i)(i) = rc(i)              (scaling by swAscale occurs here)
 *
 *          for j = 0 to i-1
 *             a(i)(j) = a(i-1)(j) + rc(i)*a(i-1)(i-j-1)
 *       }
 *
 *     See page 443, of
 *     "Digital Processing of Speech Signals" by L.R. Rabiner and R.W.
 *     Schafer; Prentice-Hall; Englewood Cliffs, NJ (USA).  1978.
 *
 *   REFERENCES: Sub-clause 4.1.7 and 4.2.3 of GSM Recomendation 06.20
 *
 *  KEYWORDS: reflectioncoefficients, parcors, conversion, rctoadp, ks, as
 *  KEYWORDS: parcorcoefficients, lpc, flat, vectorquantization
 *
 *************************************************************************/

short  rcToADp(Shortword swAscale, Shortword pswRc[],
                      Shortword pswA[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword pL_ASpace[NP],
         pL_tmpSpace[NP],
         L_temp,
        *pL_A,
        *pL_tmp,
        *pL_swap;

  short int i,
         j,                            /* loop counters */
         siLimit;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Initialize starting addresses for temporary buffers */
  /*-----------------------------------------------------*/

  pL_A = pL_ASpace;
  pL_tmp = pL_tmpSpace;

  /* Initialize the flag for checking if limiting has occured */
  /*----------------------------------------------------------*/

  siLimit = 0;

  /* Compute direct form filter coefficients, pswA[0],...,pswA[9] */
  /*-------------------------------------------------------------------*/

  for (i = 0; i < NP; i++)
  {

    pL_tmp[i] = L_mult(swAscale, pswRc[i]);
    for (j = 0; j <= i - 1; j++)
    {
      L_temp = L_mpy_ls(pL_A[i - j - 1], pswRc[i]);
      pL_tmp[j] = L_add(L_temp, pL_A[j]);
      siLimit |= isLwLimit(pL_tmp[j]);
    }
    if (i != NP - 1)
    {
      /* Swap swA and swTmp buffers */

      pL_swap = pL_tmp;
      pL_tmp = pL_A;
      pL_A = pL_swap;
    }
  }

  for (i = 0; i < NP; i++)
  {
    pswA[i] = round(pL_tmp[i]);
    siLimit |= isSwLimit(pswA[i]);
  }
  return (siLimit);
}

/***************************************************************************
 *
 *   FUNCTION NAME: rcToCorrDpL
 *
 *   PURPOSE:
 *
 *     This subroutine computes an autocorrelation vector, given a vector
 *     of reflection coefficients as an input. Double precision calculations
 *     are used internally, and a double precision (Longword)
 *     autocorrelation sequence is returned.
 *
 *   INPUTS:
 *
 *     NP
 *                     LPC filter order passed in as a global constant.
 *
 *     swAshift
 *                     Number of right shifts to be applied to the
 *                     direct form filter coefficients being computed
 *                     as an intermediate step to generating the
 *                     autocorrelation sequence.
 *
 *     swAscale
 *                     A multiplicative scale factor corresponding to
 *                     swAshift; i.e. swAscale = 2 ^(-swAshift).
 *
 *     pswRc[0:NP-1]
 *                     An input vector of reflection coefficients.
 *
 *   OUTPUTS:
 *
 *     pL_R[0:NP]
 *                     An output Longword array containing the
 *                     autocorrelation vector where
 *                     pL_R[0] = 0x7fffffff; (i.e., ~1.0).
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     The algorithm used for computing the correlation sequence is
 *     described on page 232 of the book "Linear Prediction of Speech",
 *     by J.D.  Markel and A.H. Gray, Jr.; Springer-Verlag, Berlin,
 *     Heidelberg, New York, 1976.
 *
 *   REFERENCES: Sub_Clause 4.1.4 and 4.2.1  of GSM Recomendation 06.20
 *
 *   KEYWORDS: normalized autocorrelation, reflection coefficients
 *   KEYWORDS: conversion
 *
 **************************************************************************/

void   rcToCorrDpL(Shortword swAshift, Shortword swAscale,
                          const Shortword pswRc[], Longword pL_R[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword pL_ASpace[NP],
         pL_tmpSpace[NP],
         L_temp,
         L_sum,
        *pL_A,
        *pL_tmp,
        *pL_swap;

  short int i,
         j;                            /* loop control variables */

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Set R[0] = 0x7fffffff, (i.e., R[0] = 1.0) */
  /*-------------------------------------------*/

  pL_R[0] = LW_MAX;

  /* Assign an address onto each of the two temporary buffers */
  /*----------------------------------------------------------*/

  pL_A = pL_ASpace;
  pL_tmp = pL_tmpSpace;

  /* Compute correlations R[1],...,R[10] */
  /*------------------------------------*/

  for (i = 0; i < NP; i++)
  {

    /* Compute, as an intermediate step, the filter coefficients for */
    /* for an i-th order direct form filter (pL_tmp[j],j=0,i)        */
    /*---------------------------------------------------------------*/

    pL_tmp[i] = L_mult(swAscale, pswRc[i]);
    for (j = 0; j <= i - 1; j++)
    {
      L_temp = L_mpy_ls(pL_A[i - j - 1], pswRc[i]);
      pL_tmp[j] = L_add(L_temp, pL_A[j]);
    }

    /* Swap pL_A and pL_tmp buffers */
    /*------------------------------*/

    pL_swap = pL_A;
    pL_A = pL_tmp;
    pL_tmp = pL_swap;

    /* Given the direct form filter coefficients for an i-th order filter  */
    /* and the autocorrelation vector computed up to and including stage i */
    /* compute the autocorrelation coefficient R[i+1]                      */
    /*---------------------------------------------------------------------*/

    L_temp = L_mpy_ll(pL_A[0], pL_R[i]);
    L_sum = L_negate(L_temp);

    for (j = 1; j <= i; j++)
    {
      L_temp = L_mpy_ll(pL_A[j], pL_R[i - j]);
      L_sum = L_sub(L_sum, L_temp);
    }
    pL_R[i + 1] = L_shl(L_sum, swAshift);

  }
}

/***************************************************************************
 *
 *   FUNCTION NAME: res_eng
 *
 *   PURPOSE:
 *
 *     Calculates square root of subframe residual energy estimate:
 *
 *                     sqrt( R(0)(1-k1**2)...(1-k10**2) )
 *
 *   INPUTS:
 *
 *     pswReflecCoefIn[0:9]
 *
 *                     Array of reflection coeffcients.
 *
 *     swRq
 *
 *                     Subframe energy = sqrt(frame_energy * S_LEN/2**S_SH)
 *                     (quantized).
 *
 *   OUTPUTS:
 *
 *     psnsSqrtRsOut
 *
 *                     (Pointer to) the output residual energy estimate.
 *
 *   RETURN VALUE:
 *
 *     The shift count of the normalized residual energy estimate, as int.
 *
 *   DESCRIPTION:
 *
 *     First, the canonic product of the (1-ki**2) terms is calculated
 *     (normalizations are done to maintain precision).  Also, a factor of
 *     2**S_SH is applied to the product to offset this same factor in the
 *     quantized square root of the subframe energy.
 *
 *     Then the product is square-rooted, and multiplied by the quantized
 *     square root of the subframe energy.  This combined product is put
 *     out as a normalized fraction and shift count (mantissa and exponent).
 *
 *   REFERENCES: Sub-clause 4.1.7 and 4.2.1 of GSM Recomendation 06.20
 *
 *   KEYWORDS: residualenergy, res_eng, rs
 *
 *************************************************************************/

void   res_eng(Shortword pswReflecCoefIn[], Shortword swRq,
                      struct NormSw *psnsSqrtRsOut)
{
/*_________________________________________________________________________
 |                                                                         |
 |                              Local Constants                            |
 |_________________________________________________________________________|
*/

#define  S_SH          6               /* ceiling(log2(S_LEN)) */
#define  MINUS_S_SH    -S_SH


/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Product,
         L_Shift,
         L_SqrtResEng;

  Shortword swPartialProduct,
         swPartialProductShift,
         swTerm,
         swShift,
         swSqrtPP,
         swSqrtPPShift,
         swSqrtResEng,
         swSqrtResEngShift;

  short int i;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Form canonic product, maintain precision and shift count */
  /*----------------------------------------------------------*/

  /* (Start off with unity product (actually -1), and shift offset) */
  /*----------------------------------------------------------------*/
  swPartialProduct = SW_MIN;
  swPartialProductShift = MINUS_S_SH;

  for (i = 0; i < NP; i++)
  {

    /* Get next (-1 + k**2) term, form partial canonic product */
    /*---------------------------------------------------------*/


    swTerm = mac_r(LW_MIN, pswReflecCoefIn[i], pswReflecCoefIn[i]);

    L_Product = L_mult(swTerm, swPartialProduct);

    /* Normalize partial product, round */
    /*----------------------------------*/

    swShift = norm_s(extract_h(L_Product));
    swPartialProduct = round(L_shl(L_Product, swShift));
    swPartialProductShift = add(swPartialProductShift, swShift);
  }

  /* Correct sign of product, take square root */
  /*-------------------------------------------*/

  swPartialProduct = abs_s(swPartialProduct);

  swSqrtPP = sqroot(L_deposit_h(swPartialProduct));

  L_Shift = L_shr(L_deposit_h(swPartialProductShift), 1);

  swSqrtPPShift = extract_h(L_Shift);

  if (extract_l(L_Shift) != 0)
  {

    /* Odd exponent: shr above needs to be compensated by multiplying */
    /* mantissa by sqrt(0.5)                                          */
    /*----------------------------------------------------------------*/

    swSqrtPP = mult_r(swSqrtPP, SQRT_ONEHALF);
  }

  /* Form final product, the residual energy estimate, and do final */
  /* normalization                                                  */
  /*----------------------------------------------------------------*/

  L_SqrtResEng = L_mult(swRq, swSqrtPP);

  swShift = norm_l(L_SqrtResEng);
  swSqrtResEng = round(L_shl(L_SqrtResEng, swShift));
  swSqrtResEngShift = add(swSqrtPPShift, swShift);

  /* Return */
  /*--------*/
  psnsSqrtRsOut->man = swSqrtResEng;
  psnsSqrtRsOut->sh = swSqrtResEngShift;

  return;
}

/***************************************************************************
 *
 *   FUNCTION NAME: rs_rr
 *
 *   PURPOSE:
 *
 *     Calculates sqrt(RS/R(x,x)) using floating point format,
 *     where RS is the approximate residual energy in a given
 *     subframe and R(x,x) is the power in each long term
 *     predictor vector or in each codevector.
 *     Used in the joint optimization of the gain and the long
 *     term predictor coefficient.
 *
 *   INPUTS:
 *
 *     pswExcitation[0:39] - excitation signal array
 *     snsSqrtRs - structure sqrt(RS) normalized with mantissa and shift
 *
 *   OUTPUTS:
 *
 *     snsSqrtRsRr - structure sqrt(RS/R(x,x)) with mantissa and shift
 *
 *   RETURN VALUE:
 *
 *     None
 *
 *   DESCRIPTION:
 *
 *     Implemented as sqrt(RS)/sqrt(R(x,x)) where both sqrts
 *     are stored normalized (0.5<=x<1.0) and the associated
 *     shift.  See section 4.1.11.1 for details
 *
 *   REFERENCES: Sub-clause 4.1.11.1 and 4.2.1 of GSM
 *               Recomendation 06.20
 *
 *   KEYWORDS: rs_rr, sqroot
 *
 *************************************************************************/

void   rs_rr(Shortword pswExcitation[], struct NormSw snsSqrtRs,
                    struct NormSw *snsSqrtRsRr)
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/
  Longword L_Temp;
  Shortword swTemp,
         swTemp2,
         swEnergy,
         swNormShift,
         swShift;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  swEnergy = sub(shl(snsSqrtRs.sh, 1), 3);      /* shift*2 + margin ==
                                                 * energy. */


  if (swEnergy < 0)
  {

    /* High-energy residual: scale input vector during energy */
    /* calculation. The shift count of the energy of the      */
    /* residual estimate is used as an estimate of the shift  */
    /* count needed for the excitation energy                 */
    /*--------------------------------------------------------*/

    swNormShift = g_corr1s(pswExcitation, negate(swEnergy), &L_Temp);

  }
  else
  {

    /* Lower-energy residual: no overflow protection needed */
    /*------------------------------------------------------*/

    swNormShift = g_corr1(pswExcitation, &L_Temp);
  }

  /* Compute single precision square root of energy sqrt(R(x,x)) */
  /* ----------------------------------------------------------- */
  swTemp = sqroot(L_Temp);

  /* If odd no. of shifts compensate by sqrt(0.5) */
  /* -------------------------------------------- */
  if (swNormShift & 1)
  {

    /* Decrement no. of shifts in accordance with sqrt(0.5) */
    /* ---------------------------------------------------- */
    swNormShift = sub(swNormShift, 1);

    /* sqrt(R(x,x) = sqrt(R(x,x)) * sqrt(0.5) */
    /* -------------------------------------- */
    L_Temp = L_mult(0x5a82, swTemp);
  }
  else
  {
    L_Temp = L_deposit_h(swTemp);
  }

  /* Normalize again and update shifts */
  /* --------------------------------- */
  swShift = norm_l(L_Temp);
  swNormShift = add(shr(swNormShift, 1), swShift);
  L_Temp = L_shl(L_Temp, swShift);

  /* Shift sqrt(RS) to make sure less than divisor */
  /* --------------------------------------------- */
  swTemp = shr(snsSqrtRs.man, 1);

  /* Divide sqrt(RS)/sqrt(R(x,x)) */
  /* ---------------------------- */
  swTemp2 = divide_s(swTemp, round(L_Temp));

  /* Calculate shift for division, compensate for shift before division */
  /* ------------------------------------------------------------------ */
  swNormShift = sub(snsSqrtRs.sh, swNormShift);
  swNormShift = sub(swNormShift, 1);

  /* Normalize and get no. of shifts */
  /* ------------------------------- */
  swShift = norm_s(swTemp2);
  snsSqrtRsRr->sh = add(swNormShift, swShift);
  snsSqrtRsRr->man = shl(swTemp2, swShift);

}

/***************************************************************************
 *
 *   FUNCTION NAME: rs_rrNs
 *
 *   PURPOSE:
 *
 *     Calculates sqrt(RS/R(x,x)) using floating point format,
 *     where RS is the approximate residual energy in a given
 *     subframe and R(x,x) is the power in each long term
 *     predictor vector or in each codevector.
 *
 *     Used in the joint optimization of the gain and the long
 *     term predictor coefficient.
 *
 *   INPUTS:
 *
 *     pswExcitation[0:39] - excitation signal array
 *     snsSqrtRs - structure sqrt(RS) normalized with mantissa and shift
 *
 *   OUTPUTS:
 *
 *     snsSqrtRsRr - structure sqrt(RS/R(x,x)) with mantissa and shift
 *
 *   RETURN VALUE:
 *
 *     None
 *
 *   DESCRIPTION:
 *
 *     Implemented as sqrt(RS)/sqrt(R(x,x)) where both sqrts
 *     are stored normalized (0.5<=x<1.0) and the associated
 *     shift.
 *
 *   REFERENCES: Sub-clause 4.1.11.1 and 4.2.1 of GSM
 *               Recomendation 06.20
 *
 *   KEYWORDS: rs_rr, sqroot
 *
 *************************************************************************/

void   rs_rrNs(Shortword pswExcitation[], struct NormSw snsSqrtRs,
                      struct NormSw *snsSqrtRsRr)
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/
  Longword L_Temp;
  Shortword swTemp,
         swTemp2,
         swNormShift,
         swShift;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Lower-energy residual: no overflow protection needed */
  /*------------------------------------------------------*/

  swNormShift = g_corr1(pswExcitation, &L_Temp);


  /* Compute single precision square root of energy sqrt(R(x,x)) */
  /* ----------------------------------------------------------- */
  swTemp = sqroot(L_Temp);

  /* If odd no. of shifts compensate by sqrt(0.5) */
  /* -------------------------------------------- */
  if (swNormShift & 1)
  {

    /* Decrement no. of shifts in accordance with sqrt(0.5) */
    /* ---------------------------------------------------- */
    swNormShift = sub(swNormShift, 1);

    /* sqrt(R(x,x) = sqrt(R(x,x)) * sqrt(0.5) */
    /* -------------------------------------- */
    L_Temp = L_mult(0x5a82, swTemp);
  }
  else
  {
    L_Temp = L_deposit_h(swTemp);
  }

  /* Normalize again and update shifts */
  /* --------------------------------- */

  swShift = norm_l(L_Temp);
  swNormShift = add(shr(swNormShift, 1), swShift);
  L_Temp = L_shl(L_Temp, swShift);

  /* Shift sqrt(RS) to make sure less than divisor */
  /* --------------------------------------------- */
  swTemp = shr(snsSqrtRs.man, 1);

  /* Divide sqrt(RS)/sqrt(R(x,x)) */
  /* ---------------------------- */
  swTemp2 = divide_s(swTemp, round(L_Temp));

  /* Calculate shift for division, compensate for shift before division */
  /* ------------------------------------------------------------------ */
  swNormShift = sub(snsSqrtRs.sh, swNormShift);
  swNormShift = sub(swNormShift, 1);

  /* Normalize and get no. of shifts */
  /* ------------------------------- */
  swShift = norm_s(swTemp2);
  snsSqrtRsRr->sh = add(swNormShift, swShift);
  snsSqrtRsRr->man = shl(swTemp2, swShift);

}


/***************************************************************************
 *
 *   FUNCTION NAME: scaleExcite
 *
 *   PURPOSE:
 *
 *     Scale an arbitrary excitation vector (codevector or
 *     pitch vector)
 *
 *   INPUTS:
 *
 *     pswVect[0:39] - the unscaled vector.
 *     iGsp0Scale - an positive offset to compensate for the fact
 *                  that GSP0 table is scaled down.
 *     swErrTerm - rather than a gain being passed in, (beta, gamma)
 *                 it is calculated from this error term - either
 *                 Gsp0[][][0] error term A or Gsp0[][][1] error
 *                 term B. Beta is calculated from error term A,
 *                 gamma from error term B.
 *     snsRS - the RS_xx appropriate to pswVect.
 *
 *   OUTPUTS:
 *
 *      pswScldVect[0:39] - the output, scaled excitation vector.
 *
 *   RETURN VALUE:
 *
 *     swGain - One of two things.  Either a clamped value of 0x7fff if the
 *              gain's shift was > 0 or the rounded vector gain otherwise.
 *
 *   DESCRIPTION:
 *
 *     If gain > 1.0 then
 *         (do not shift gain up yet)
 *         partially scale vector element THEN shift and round save
 *      else
 *         shift gain correctly
 *         scale vector element and round save
 *         update state array
 *
 *   REFERENCES: Sub-clause 4.1.10.2 and 4.2.1 of GSM
 *               Recomendation 06.20
 *
 *   KEYWORDS: excite_vl, sc_ex, excitevl, scaleexcite, codevector, p_vec,
 *   KEYWORDS: x_vec, pitchvector, gain, gsp0
 *
 *************************************************************************/

Shortword scaleExcite(Shortword pswVect[],
                             Shortword swErrTerm, struct NormSw snsRS,
                             Shortword pswScldVect[])
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/
  Longword L_GainUs,
         L_scaled,
         L_Round_off;
  Shortword swGain,
         swGainUs,
         swGainShift,
         i,
         swGainUsShft;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/


  L_GainUs = L_mult(swErrTerm, snsRS.man);
  swGainUsShft = norm_s(extract_h(L_GainUs));
  L_GainUs = L_shl(L_GainUs, swGainUsShft);

  swGainShift = add(swGainUsShft, snsRS.sh);
  swGainShift = sub(swGainShift, GSP0_SCALE);


  /* gain > 1.0 */
  /* ---------- */

  if (swGainShift < 0)
  {
    swGainUs = round(L_GainUs);

    L_Round_off = L_shl((long) 32768, swGainShift);

    for (i = 0; i < S_LEN; i++)
    {
      L_scaled = L_mac(L_Round_off, swGainUs, pswVect[i]);
      L_scaled = L_shr(L_scaled, swGainShift);
      pswScldVect[i] = extract_h(L_scaled);
    }

    if (swGainShift == 0)
      swGain = swGainUs;
    else
      swGain = 0x7fff;
  }

  /* gain < 1.0 */
  /* ---------- */

  else
  {

    /* shift down or not at all */
    /* ------------------------ */
    if (swGainShift > 0)
      L_GainUs = L_shr(L_GainUs, swGainShift);

    /* the rounded actual vector gain */
    /* ------------------------------ */
    swGain = round(L_GainUs);

    /* now scale the vector (with rounding) */
    /* ------------------------------------ */

    for (i = 0; i < S_LEN; i++)
    {
      L_scaled = L_mac((long) 32768, swGain, pswVect[i]);
      pswScldVect[i] = extract_h(L_scaled);
    }
  }
  return (swGain);
}

/***************************************************************************
 *
 *   FUNCTION NAME: sqroot
 *
 *   PURPOSE:
 *
 *     The purpose of this function is to perform a single precision square
 *     root function on a Longword
 *
 *   INPUTS:
 *
 *     L_SqrtIn
 *                     input to square root function
 *
 *   OUTPUTS:
 *
 *     none
 *
 *   RETURN VALUE:
 *
 *     swSqrtOut
 *                     output to square root function
 *
 *   DESCRIPTION:
 *
 *      Input assumed to be normalized
 *
 *      The algorithm is based around a six term Taylor expansion :
 *
 *        y^0.5 = (1+x)^0.5
 *             ~= 1 + (x/2) - 0.5*((x/2)^2) + 0.5*((x/2)^3)
 *                - 0.625*((x/2)^4) + 0.875*((x/2)^5)
 *
 *      Max error less than 0.08 % for normalized input ( 0.5 <= x < 1 )
 *
 *   REFERENCES: Sub-clause 4.1.4.1, 4.1.7, 4.1.11.1, 4.2.1,
 *               4.2.2, 4.2.3 and 4.2.4 of GSM Recomendation 06.20
 *
 *   KEYWORDS: sqrt, squareroot, sqrt016
 *
 *************************************************************************/

Shortword sqroot(Longword L_SqrtIn)
{

/*_________________________________________________________________________
 |                                                                         |
 |                              Local Constants                            |
 |_________________________________________________________________________|
*/

#define    PLUS_HALF          0x40000000L       /* 0.5 */
#define    MINUS_ONE          0x80000000L       /* -1 */
#define    TERM5_MULTIPLER    0x5000   /* 0.625 */
#define    TERM6_MULTIPLER    0x7000   /* 0.875 */

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Temp0,
         L_Temp1;

  Shortword swTemp,
         swTemp2,
         swTemp3,
         swTemp4,
         swSqrtOut;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* determine 2nd term x/2 = (y-1)/2 */
  /* -------------------------------- */

  L_Temp1 = L_shr(L_SqrtIn, 1);        /* L_Temp1 = y/2 */
  L_Temp1 = L_sub(L_Temp1, PLUS_HALF); /* L_Temp1 = (y-1)/2 */
  swTemp = extract_h(L_Temp1);         /* swTemp = x/2 */

  /* add contribution of 2nd term */
  /* ---------------------------- */

  L_Temp1 = L_sub(L_Temp1, MINUS_ONE); /* L_Temp1 = 1 + x/2 */

  /* determine 3rd term */
  /* ------------------ */

  L_Temp0 = L_msu(0L, swTemp, swTemp); /* L_Temp0 = -(x/2)^2 */
  swTemp2 = extract_h(L_Temp0);        /* swTemp2 = -(x/2)^2 */
  L_Temp0 = L_shr(L_Temp0, 1);         /* L_Temp0 = -0.5*(x/2)^2 */

  /* add contribution of 3rd term */
  /* ---------------------------- */

  L_Temp0 = L_add(L_Temp1, L_Temp0);   /* L_Temp0 = 1 + x/2 - 0.5*(x/2)^2 */

  /* determine 4rd term */
  /* ------------------ */

  L_Temp1 = L_msu(0L, swTemp, swTemp2);/* L_Temp1 = (x/2)^3 */
  swTemp3 = extract_h(L_Temp1);        /* swTemp3 = (x/2)^3 */
  L_Temp1 = L_shr(L_Temp1, 1);         /* L_Temp1 = 0.5*(x/2)^3 */

  /* add contribution of 4rd term */
  /* ---------------------------- */

  /* L_Temp1 = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 */

  L_Temp1 = L_add(L_Temp0, L_Temp1);

  /* determine partial 5th term */
  /* -------------------------- */

  L_Temp0 = L_mult(swTemp, swTemp3);   /* L_Temp0 = (x/2)^4 */
  swTemp4 = round(L_Temp0);            /* swTemp4 = (x/2)^4 */

  /* determine partial 6th term */
  /* -------------------------- */

  L_Temp0 = L_msu(0L, swTemp2, swTemp3);        /* L_Temp0 = (x/2)^5 */
  swTemp2 = round(L_Temp0);            /* swTemp2 = (x/2)^5 */

  /* determine 5th term and add its contribution */
  /* ------------------------------------------- */

  /* L_Temp0 = -0.625*(x/2)^4 */

  L_Temp0 = L_msu(0L, TERM5_MULTIPLER, swTemp4);

  /* L_Temp1 = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 - 0.625*(x/2)^4 */

  L_Temp1 = L_add(L_Temp0, L_Temp1);

  /* determine 6th term and add its contribution */
  /* ------------------------------------------- */

  /* swSqrtOut = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 */
  /* - 0.625*(x/2)^4 + 0.875*(x/2)^5     */

  swSqrtOut = mac_r(L_Temp1, TERM6_MULTIPLER, swTemp2);

  /* return output */
  /* ------------- */

  return (swSqrtOut);
}

/***************************************************************************
 *
 *   FUNCTION NAME: v_con
 *
 *   PURPOSE:
 *
 *     This subroutine constructs a codebook excitation
 *     vector from basis vectors
 *
 *   INPUTS:
 *
 *     pswBVects[0:siNumBVctrs*S_LEN-1]
 *
 *                     Array containing a set of basis vectors.
 *
 *     pswBitArray[0:siNumBVctrs-1]
 *
 *                     Bit array dictating the polarity of the
 *                     basis vectors in the output vector.
 *                     Each element of the bit array is either -0.5 or +0.5
 *
 *     siNumBVctrs
 *                     Number of bits in codeword
 *
 *   OUTPUTS:
 *
 *     pswOutVect[0:39]
 *
 *                     Array containing the contructed output vector
 *
 *   RETURN VALUE:
 *
 *     none
 *
 *   DESCRIPTION:
 *
 *     The array pswBitArray is used to multiply each of the siNumBVctrs
 *     basis vectors.  The input pswBitArray[] is an array whose
 *     elements are +/-0.5.  These multiply the VSELP basis vectors and
 *     when summed produce a VSELP codevector.  b_con() is the function
 *     used to translate a VSELP codeword into pswBitArray[].
 *
 *
 *   REFERENCES: Sub-clause 4.1.10 and 4.2.1 of GSM Recomendation 06.20
 *
 *   KEYWORDS: v_con, codeword, reconstruct, basis vector, excitation
 *
 *************************************************************************/

void   v_con(const Shortword pswBVects[], Shortword pswOutVect[],
                    Shortword pswBitArray[], short int siNumBVctrs)
{

/*_________________________________________________________________________
 |                                                                         |
 |                            Automatic Variables                          |
 |_________________________________________________________________________|
*/

  Longword L_Temp;

  short int siSampleCnt,
         siCVectCnt;

/*_________________________________________________________________________
 |                                                                         |
 |                              Executable Code                            |
 |_________________________________________________________________________|
*/

  /* Sample loop  */
  /*--------------*/
  for (siSampleCnt = 0; siSampleCnt < S_LEN; siSampleCnt++)
  {

    /* First element of output vector */
    L_Temp = L_mult(pswBitArray[0], pswBVects[0 * S_LEN + siSampleCnt]);

    /* Construct output vector */
    /*-------------------------*/
    for (siCVectCnt = 1; siCVectCnt < siNumBVctrs; siCVectCnt++)
    {
      L_Temp = L_mac(L_Temp, pswBitArray[siCVectCnt],
                     pswBVects[siCVectCnt * S_LEN + siSampleCnt]);
    }

    /* store the output vector sample */
    /*--------------------------------*/
    L_Temp = L_shl(L_Temp, 1);
    pswOutVect[siSampleCnt] = extract_h(L_Temp);
  }
}