语音压缩

开发平台：
C/C++

LSP.C：源码内容
							/*
**
** File:    lsp.c
**
** Description: Functions that implement line spectral pair 
**      (LSP) operations.  
**
** Functions:
**
**  Converting between linear predictive coding (LPC) coefficients
**  and LSP frequencies:
**
**      AtoLsp()
**      LsptoA()
**
**  Vector quantization (VQ) of LSP frequencies:
**
**      Lsp_Qnt()
**      Lsp_Svq()
**      Lsp_Inq()
**
**  Interpolation of LSP frequencies:
**
**      Lsp_Int()
*/
/*
    ITU-T G.723 Speech Coder   ANSI-C Source Code     Version 5.00
    copyright (c) 1995, AudioCodes, DSP Group, France Telecom,
    Universite de Sherbrooke.  All rights reserved.
*/
#include <stdio.h>
#include "typedef.h"
#include "basop.h"
#include "cst_lbc.h"
#include "tab_lbc.h"
#include "lsp.h"
/*
**
** Function:            AtoLsp()
**
** Description:     Transforms 10 LPC coefficients to the 10
**          corresponding LSP frequencies for a subframe.
**          This transformation is done once per frame,
**          for subframe 3 only.  The transform algorithm
**          generates sum and difference polynomials from
**          the LPC coefficients.  It then evaluates the
**          sum and difference polynomials at uniform
**          intervals of pi/256 along the unit circle.
**          Intervals where a sign change occurs are
**          interpolated to find the zeros of the
**          polynomials, which are the LSP frequencies.
**
** Links to text:   Section 2.5
**
** Arguments:       
**
**  Word16 *LspVect     Empty Buffer
**  Word16 Lpc[]        Unquantized LPC coefficients (10 words)
**  Word16 PrevLsp[]    LSP frequencies from the previous frame (10 words)
**
** Outputs:
**
**  Word16 LspVect[]    LSP frequencies for the current frame (10 words)
**
** Return value:        None
**
**/
void AtoLsp( Word16 *LspVect, Word16 *Lpc, Word16 *PrevLsp )
{
    int   i,j,k ;
    Word32   Lpq[LpcOrder+2] ;
    Word16   Spq[LpcOrder+2] ;
    Word16   Exp   ;
    Word16   LspCnt ;
    Word32   PrevVal,CurrVal   ;
    Word32   Acc0,Acc1   ;
 /*
  * Perform a bandwidth expansion on the LPC coefficients.  This
  * scales the poles of the LPC synthesis filter by a factor of
  * 0.994.
  */
    for ( i = 0 ; i < LpcOrder ; i ++ )
        LspVect[i] = mult_r( Lpc[i], BandExpTable[i] ) ;
 /*
  * Compute the sum and difference polynomials with the roots at z =
  * -1 (sum) or z = +1 (difference) removed.  Let these polynomials
  * be P(z) and Q(z) respectively, and let their coefficients be
  * {p_i} amd {q_i}.  The coefficients are stored in the array Lpq[]
  * as follows: p_0, q_0, p_1, q_1, ..., p_5, q_5.  There is no need
  * to store the other coefficients because of symmetry.
  */
 /*
  * Set p_0 = q_0 = 1.  The LPC coefficients are already scaled by
  *  1/4.  P(z) and Q(z) are scaled by an additional scaling factor of
  *  1/16, for an overall factor of 1/64 = 0x02000000L.
  */
    Lpq[0] = Lpq[1] = (Word32) 0x02000000L ;
 /*
  * This loop computes the coefficients of P(z) and Q(z).  The long
  * division (to remove the real zeros) is done recursively.
  */
    for ( i = 0 ; i < LpcOrder/2 ; i ++ ) {
        /* P(z) */
        Acc0 = L_negate( Lpq[2*i+0] ) ;
        Acc1 = L_deposit_h( LspVect[i] ) ;
        Acc1 = L_shr( Acc1, (Word16) 4 ) ;
        Acc0 = L_sub( Acc0, Acc1 ) ;
        Acc1 = L_deposit_h( LspVect[LpcOrder-1-i] ) ;
        Acc1 = L_shr( Acc1, (Word16) 4 ) ;
        Acc0 = L_sub( Acc0, Acc1 ) ;
        Lpq[2*i+2] = Acc0 ;
        /* Q(z) */
        Acc0 = Lpq[2*i+1] ;
        Acc1 = L_deposit_h( LspVect[i] ) ;
        Acc1 = L_shr( Acc1, (Word16) 4 ) ;
        Acc0 = L_sub( Acc0, Acc1 ) ;
        Acc1 = L_deposit_h( LspVect[LpcOrder-1-i] ) ;
        Acc1 = L_shr( Acc1, (Word16) 4 ) ;
        Acc0 = L_add( Acc0, Acc1 ) ;
        Lpq[2*i+3] = Acc0 ;
    }
 /*
  * Divide p_5 and q_5 by 2 for proper weighting during polynomial
  * evaluation.
  */
    Lpq[LpcOrder+0] = L_shr( Lpq[LpcOrder+0], (Word16) 1 ) ;
    Lpq[LpcOrder+1] = L_shr( Lpq[LpcOrder+1], (Word16) 1 ) ;
 /*
  * Normalize the polynomial coefficients and convert to shorts
  */
    /* Find the maximum */
    Acc1 = L_abs( Lpq[0] ) ;
    for ( i = 1 ; i < LpcOrder+2 ; i ++ ) {
        Acc0 = L_abs( Lpq[i] ) ;
        if ( Acc0 > Acc1 )
            Acc1 = Acc0 ;
    }
    /* Compute the normalization factor */
    Exp = norm_l( Acc1 ) ;
    /* Normalize and convert to shorts */
    for ( i = 0 ; i < LpcOrder+2 ; i ++ ) {
        Acc0 = L_shl( Lpq[i], Exp ) ;
        Spq[i] = round( Acc0 ) ;
    }
 /*
  * Initialize the search loop
  */
 /*
  * The variable k is a flag that indicates which polynomial (sum or
  * difference) the algorithm is currently evaluating.  Start with
  * the sum.
  */
    k = 0 ;
    /* Evaluate the sum polynomial at frequency zero */
    PrevVal = (Word32) 0 ;
    for ( j = 0 ; j <= LpcOrder/2 ; j ++ )
        PrevVal = L_mac( PrevVal, Spq[2*j], CosineTable[0] ) ;
 /*
  * Search loop.  Evaluate P(z) and Q(z) at uniform intervals of
  * pi/256 along the unit circle.  Check for zero crossings.  The
  * zeros of P(w) and Q(w) alternate, so only one of them need by
  * evaluated at any given step.
  */
    LspCnt = (Word16) 0 ;
    for ( i = 1 ; i < CosineTableSize/2 ; i ++ ) {
        /* Evaluate the selected polynomial */
        CurrVal = (Word32) 0 ;
        for ( j = 0 ; j <= LpcOrder/2 ; j ++ )
            CurrVal = L_mac( CurrVal, Spq[LpcOrder-2*j+k],
                                    CosineTable[i*j%CosineTableSize] ) ;
        /* Check for a sign change, indicating a zero crossing */
        if ( (CurrVal ^ PrevVal) < (Word32) 0 ) {
 /*
  * Interpolate to find the bottom 7 bits of the
  * zero-crossing frequency
  */
            Acc0 = L_abs( CurrVal ) ;
            Acc1 = L_abs( PrevVal ) ;
            Acc0 = L_add( Acc0, Acc1 ) ;
            /* Normalize the sum */
            Exp = norm_l( Acc0 ) ;
            Acc0 = L_shl( Acc0, Exp ) ;
            Acc1 = L_shl( Acc1, Exp ) ;
            Acc1 = L_shr( Acc1, (Word16) 8 ) ;
            LspVect[LspCnt] = div_l( Acc1, extract_h( Acc0 ) ) ;
 /*
  * Add the upper part of the zero-crossing frequency,
  * i.e. bits 7-15
  */
            Exp = shl( (Word16) (i-1), (Word16) 7 ) ;
            LspVect[LspCnt] = add( LspVect[LspCnt], Exp ) ;
            LspCnt ++ ;
            /* Check if all zeros have been found */
            if ( LspCnt == (Word16) LpcOrder )
                break ;
 /*
  * Switch the pointer between sum and difference polynomials
  */
            k ^= 1 ;
 /*
  * Evaluate the new polynomial at the current frequency
  */
            CurrVal = (Word32) 0 ;
            for ( j = 0 ; j <= LpcOrder/2 ; j ++ )
                CurrVal = L_mac( CurrVal, Spq[LpcOrder-2*j+k],
                                    CosineTable[i*j%CosineTableSize] ) ;
        }
        /* Update the previous value */
        PrevVal = CurrVal ;
    }
 /*
  * Check if all 10 zeros were found.  If not, ignore the results of
  * the search and use the previous frame's LSP frequencies instead.
  */
    if ( LspCnt != (Word16) LpcOrder ) {
        for ( j = 0 ; j < LpcOrder ; j ++ )
            LspVect[j] = PrevLsp[j] ;
    }
    return ;
}
/*
**
** Function:            Lsp_Qnt()
**
** Description:     Vector quantizes the LSP frequencies.  The LSP
**          vector is divided into 3 sub-vectors, or
**          bands, of dimension 3, 3, and 4.  Each band is
**          quantized separately using a different VQ
**          table.  Each table has 256 entries, so the
**          quantization generates three indices of 8 bits
**          each.  (Only the LSP vector for subframe 3 is
**          quantized per frame.)
**
** Links to text:   Section 2.5
**
** Arguments:       
**
**  Word16 CurrLsp[]    Unquantized LSP frequencies for the current frame (10 words)
**  Word16 PrevLsp[]    LSP frequencies from the previous frame (10 words)
**
** Outputs:             Quantized LSP frequencies for the current frame (10 words)
**
** Return value:
**
**  Word32      Long word packed with the 3 VQ indices.  Band 0
**          corresponds to bits [23:16], band 1 corresponds
**          to bits [15:8], and band 2 corresponds to bits [7:0].
**          (Bit 0 is the least significant.)
**
*/
Word32   Lsp_Qnt( Word16 *CurrLsp, Word16 *PrevLsp )
{
    int   i ;
    Word16   Wvect[LpcOrder] ;
    Word16   Tmp0,Tmp1   ;
    Word16   Exp   ;
 /*
  * Compute the VQ weighting vector.  The weights assign greater
  * precision to those frequencies that are closer together.
  */
    /* Compute the end differences */
    Wvect[0] = sub( CurrLsp[1], CurrLsp[0] ) ;
    Wvect[LpcOrder-1] = sub( CurrLsp[LpcOrder-1], CurrLsp[LpcOrder-2] ) ;
    /* Compute the rest of the differences */
    for ( i = 1 ; i < LpcOrder-1 ; i ++ ) {
        Tmp0 = sub( CurrLsp[i+1], CurrLsp[i] ) ;
        Tmp1 = sub( CurrLsp[i], CurrLsp[i-1] ) ;
        if ( Tmp0 > Tmp1 )
            Wvect[i] = Tmp1 ;
        else
            Wvect[i] = Tmp0 ;
    }
    /* Invert the differences */
    Tmp0 = (Word16) 0x0020 ;
    for ( i = 0 ; i < LpcOrder ; i ++ ) {
        if ( Wvect[i] > Tmp0 )
            Wvect[i] = div_s( Tmp0, Wvect[i] ) ;
        else
            Wvect[i] = MAX_16 ;
    }
    /* Normalize the weight vector */
    Tmp0 = (Word16) 0 ;
    for ( i = 0 ; i < LpcOrder ; i ++ )
        if ( Wvect[i] > Tmp0 )
            Tmp0 = Wvect[i] ;
    Exp = norm_s( Tmp0 ) ;
    for ( i = 0 ; i < LpcOrder ; i ++ )
        Wvect[i] = shl( Wvect[i], Exp ) ;
 /*
  * Compute the VQ target vector.  This is the residual that remains
  * after subtracting both the DC and predicted
  * components.
  */
 /*
  * Subtract the DC component from both the current and previous LSP
  * vectors.
  */
    for ( i = 0 ; i < LpcOrder ; i ++ ) {
        CurrLsp[i] = sub( CurrLsp[i], LspDcTable[i] ) ;
        PrevLsp[i] = sub( PrevLsp[i], LspDcTable[i] ) ;
    }
 /*
  * Generate the prediction vector and subtract it.  Use a constant
  * first-order predictor based on the previous (DC-free) LSP
  * vector.
  */
    for ( i = 0 ; i < LpcOrder ; i ++ ) {
        Tmp0 = mult_r( PrevLsp[i], (Word16) LspPrd0 ) ;
        CurrLsp[i] = sub( CurrLsp[i], Tmp0 ) ;
    }
 /*
  * Add the DC component back to the previous LSP vector.  This
  * vector is needed in later routines.
  */
    for ( i = 0 ; i < LpcOrder ; i ++ )
        PrevLsp[i] = add( PrevLsp[i], LspDcTable[i] ) ;
 /*
  * Do the vector quantization for all three bands
  */
    return Lsp_Svq( CurrLsp, Wvect ) ;
}
/*
**
** Function:            Lsp_Svq()
**
** Description:     Performs the search of the VQ tables to find
**          the optimum LSP indices for all three bands.
**          For each band, the search finds the index which 
**          minimizes the weighted squared error between 
**          the table entry and the target.
**
** Links to text:   Section 2.5
**
** Arguments:       
**
**  Word16 Tv[]     VQ target vector (10 words)
**  Word16 Wvect[]      VQ weight vector (10 words)
**
** Outputs:         None
**
** Return value:    
**
**  Word32      Long word packed with the 3 VQ indices.  Band 0
**          corresponds to bits [23:16], band 1 corresponds
**          to bits [15:8], and band 2 corresponds to bits [7:0].
**              
*/
Word32   Lsp_Svq( Word16 *Tv, Word16 *Wvect )
{
    int   i,j,k ;
    Word32   Rez,Indx    ;
    Word32   Acc0,Acc1   ;
    Word16   Tmp[LpcOrder] ;
    Word16  *LspQntPnt  ;
 /*
  * Initialize the return value
  */
    Rez = (Word32) 0 ;
 /*
  * Quantize each band separately
  */
    for ( k = 0 ; k < LspQntBands ; k ++ ) {
 /*
  * Search over the entire VQ table to find the index that
  * minimizes the error.
  */
        /* Initialize the search */
        Acc1 = (Word32) -1 ;
        Indx = (Word32) 0 ;
        LspQntPnt = BandQntTable[k] ;
        for ( i = 0 ; i < LspCbSize ; i ++ ) {
 /*
  * Generate the metric, which is the negative error with the
  * constant component removed.
  */
            for ( j = 0 ; j < BandInfoTable[k][1] ; j ++ )
                Tmp[j] = mult_r( Wvect[BandInfoTable[k][0]+j],
                                                            LspQntPnt[j] ) ;
            Acc0 = (Word32) 0 ;
            for ( j = 0 ; j < BandInfoTable[k][1] ; j ++ )
                Acc0 = L_mac( Acc0, Tv[BandInfoTable[k][0]+j], Tmp[j] ) ;
            Acc0 = L_shl( Acc0, (Word16) 1 ) ;
            for ( j = 0 ; j < BandInfoTable[k][1] ; j ++ )
                Acc0 = L_msu( Acc0, LspQntPnt[j], Tmp[j] ) ;
            LspQntPnt += BandInfoTable[k][1] ;
 /*
  * Compare the metric to the previous maximum and select the
  * new index
  */
            if ( Acc0 > Acc1 ) {
                Acc1 = Acc0 ;
                Indx = (Word32) i ;
            }
        }
 /*
  * Pack the result with the optimum index for this band
  */
        Rez = L_shl( Rez, (Word16) LspCbBits ) ;
        Rez = L_add( Rez, Indx ) ;
    }
    return Rez ;
}
/*
**
** Function:            Lsp_Inq()
**
** Description:     Performs inverse vector quantization of the
**          LSP frequencies.  The LSP vector is divided
**          into 3 sub-vectors, or bands, of dimension 3,
**          3, and 4.  Each band is inverse quantized
**          separately using a different VQ table.  Each
**          table has 256 entries, so each VQ index is 8
**          bits.  (Only the LSP vector for subframe 3 is
**          quantized per frame.)
**
** Links to text:   Sections 2.6, 3.2
**
** Arguments:
**
**  Word16 *Lsp     Empty buffer
**  Word16 PrevLsp[]    Quantized LSP frequencies from the previous frame
**               (10 words)
**  Word32 LspId        Long word packed with the 3 VQ indices.  Band 0
**               corresponds to bits [23:16], band 1 corresponds
**               to bits [15:8], and band 2 corresponds to bits
**               [7:0].
**  Word16 Crc      Frame erasure indicator
**
** Outputs:
**
**  Word16 Lsp[]        Quantized LSP frequencies for current frame (10
**               words)
**
** Return value:         None
**
*/
void Lsp_Inq( Word16 *Lsp, Word16 *PrevLsp, Word32 LspId, Word16 Crc )
{
    int  i,j   ;
    Word16  *LspQntPnt  ;
    Word16   Scon  ;
    Word16   Lprd  ;
    Word16   Tmp   ;
    Flag     Test  ;
 /*
  * Check for frame erasure.  If a frame erasure has occurred, the
  * resulting VQ table entries are zero.  In addition, a different
  * fixed predictor and minimum frequency separation are used.
  */
    if ( Crc == (Word16) 0 ) {
        Scon = (Word16) 0x0100 ;
        Lprd = LspPrd0 ;
    }
    else {
        LspId = (Word32) 0 ;
        Scon = (Word16) 0x0200 ;
        Lprd = LspPrd1 ;
    }
 /*
  * Inverse quantize the 10th-order LSP vector.  Each band is done
  * separately.
  */
    for ( i = LspQntBands-1; i >= 0 ; i -- ) {
 /*
  * Get the VQ table entry corresponding to the transmitted index
  */
        Tmp = (Word16) ( LspId & (Word32) 0x000000ff ) ;
        LspId >>= 8 ;
        LspQntPnt = BandQntTable[i] ;
        for ( j = 0 ; j < BandInfoTable[i][1] ; j ++ )
            Lsp[BandInfoTable[i][0] + j] =
                                LspQntPnt[Tmp*BandInfoTable[i][1] + j] ;
    }
 /*
  * Subtract the DC component from the previous frame's quantized
  * vector
  */
    for ( j = 0 ; j < LpcOrder ; j ++ )
        PrevLsp[j] = sub(PrevLsp[j], LspDcTable[j] ) ;
 /*
  * Generate the prediction vector using a fixed first-order
  * predictor based on the previous frame's (DC-free) quantized
  * vector
  */
    for ( j = 0 ; j < LpcOrder ; j ++ ) {
        Tmp = mult_r( PrevLsp[j], Lprd ) ;
        Lsp[j] = add( Lsp[j], Tmp ) ;
    }
 /*
  * Add the DC component back to the previous quantized vector,
  * which is needed in later routines
  */
    for ( j = 0 ; j < LpcOrder ; j ++ ) {
        PrevLsp[j] = add( PrevLsp[j], LspDcTable[j] ) ;
        Lsp[j] = add( Lsp[j], LspDcTable[j] ) ;
    }
 /*
  * Perform a stability test on the quantized LSP frequencies.  This
  * test checks that the frequencies are ordered, with a minimum
  * separation between each.  If the test fails, the frequencies are
  * iteratively modified until the test passes.  If after 10
  * iterations the test has not passed, the previous frame's
  * quantized LSP vector is used.
  */
    for ( i = 0 ; i < LpcOrder ; i ++ ) {
        /* Check the first frequency */
        if ( Lsp[0] < (Word16) 0x180 )
            Lsp[0] = (Word16) 0x180 ;
        /* Check the last frequency */
        if ( Lsp[LpcOrder-1] > (Word16) 0x7e00 )
            Lsp[LpcOrder-1] = (Word16) 0x7e00 ;
        /* Perform the modification */
        for ( j = 1 ; j < LpcOrder ; j ++ ) {
            Tmp = add( Scon, Lsp[j-1] ) ;
            Tmp = sub( Tmp, Lsp[j] ) ;
            if ( Tmp > (Word16) 0 ) {
                Tmp = shr( Tmp, (Word16) 1 ) ;
                Lsp[j-1] = sub( Lsp[j-1], Tmp ) ;
                Lsp[j] = add( Lsp[j], Tmp ) ;
            }
        }
        Test = False ;
 /*
  * Test the modified frequencies for stability.  Break out of
  * the loop if the frequencies are stable.
  */
        for ( j = 1 ; j < LpcOrder ; j ++ ) {
            Tmp = add( Lsp[j-1], Scon ) ;
            Tmp = sub( Tmp, (Word16) 4 ) ;
            Tmp = sub( Tmp, Lsp[j] ) ;
            if ( Tmp > (Word16) 0 )
                Test = True ;
        }
        if ( Test == False )
            break ;
    }
 /*
  * Return the result of the stability check.  True = not stable,
  * False = stable.
  */
    if ( Test == True) {
        for ( j = 0 ; j < LpcOrder ; j ++ )
            Lsp[j] = PrevLsp[j] ;
    }
    return;
}
/*
**
** Function:            Lsp_Int()
**
** Description:     Computes the quantized LPC coefficients for a
**          frame.  First the quantized LSP frequencies
**          for all subframes are computed by linear
**          interpolation.  These frequencies are then
**          transformed to quantized LPC coefficients.
**
** Links to text:   Sections 2.7, 3.3
**
** Arguments:
**
**  Word16 *QntLpc      Empty buffer
**  Word16 CurrLsp[]    Quantized LSP frequencies for the current frame,
**               subframe 3 (10 words)
**  Word16 PrevLsp[]    Quantized LSP frequencies for the previous frame,
**               subframe 3 (10 words)
**
** Outputs:
**
**  Word16 QntLpc[]     Quantized LPC coefficients for current frame, all
**               subframes (40 words)
**
** Return value:        None
**
*/
void  Lsp_Int( Word16 *QntLpc, Word16 *CurrLsp, Word16 *PrevLsp )
{
    int   i,j   ;
    Word16   Tmp   ;
    Word16  *Dpnt  ;
    Word32   Acc0  ;
 /*
  * Initialize the interpolation factor
  */
    Tmp = (Word16) (MIN_16 / SubFrames ) ;
    Dpnt = QntLpc ;
 /*
  * Do for all subframes
  */
    for ( i = 0 ; i < SubFrames ; i ++ ) {
 /*
  * Compute the quantized LSP frequencies by linear interpolation
  * of the frequencies from subframe 3 of the current and
  * previous frames
  */
        for ( j = 0 ; j < LpcOrder ; j ++ ) {
            Acc0 = L_deposit_h( PrevLsp[j] ) ;
            Acc0 = L_mac( Acc0, Tmp, PrevLsp[j] ) ;
            Acc0 = L_msu( Acc0, Tmp, CurrLsp[j] ) ;
            Dpnt[j] = round( Acc0 ) ;
        }
 /*
  * Convert the quantized LSP frequencies to quantized LPC
  * coefficients
  */
        LsptoA( Dpnt ) ;
        Dpnt += LpcOrder ;
        /* Update the interpolation factor */
        Tmp = add( Tmp, (Word16) (MIN_16 / SubFrames ) ) ;
    }
}
/*
**
** Function:            LsptoA()
**
** Description:     Converts LSP frequencies to LPC coefficients
**          for a subframe.  Sum and difference
**          polynomials are computed from the LSP
**          frequencies (which are the roots of these
**          polynomials).  The LPC coefficients are then
**          computed by adding the sum and difference
**          polynomials.
**          
** Links to text:   Sections 2.7, 3.3
**
** Arguments:       
**
**  Word16 Lsp[]        LSP frequencies (10 words)
**
** Outputs:
**
**  Word16 Lsp[]        LPC coefficients (10 words)
**
** Return value:        None
** 
*/
void  LsptoA( Word16 *Lsp )
{
    int   i,j   ;
    Word32   Acc0,Acc1   ;
    Word16   Tmp ;
    Word32   P[LpcOrder/2+1] ;
    Word32   Q[LpcOrder/2+1] ;
 /*
  * Compute the cosines of the LSP frequencies by table lookup and
  * linear interpolation
  */
    for ( i = 0 ; i < LpcOrder ; i ++ ) {
 /*
  * Do the table lookup using bits [15:7] of the LSP frequency
  */
        j = (int) shr( Lsp[i], (Word16) 7 ) ;
        Acc0 = L_deposit_h( CosineTable[j] ) ;
 /*
  * Do the linear interpolations using bits [6:0] of the LSP
  * frequency
  */
        Tmp = sub(CosineTable[j+1], CosineTable[j] ) ;
        Acc0 = L_mac( Acc0, Tmp, add( shl( (Word16)(Lsp[i] & 0x007f) ,
                                (Word16)8 ), (Word16) 0x0080 ) ) ;
        Acc0 = L_shl( Acc0, (Word16) 1 ) ;
        Lsp[i] = negate( round( Acc0 ) ) ;
    }
 /*
  * Compute the sum and difference polynomials with the real roots
  * removed.  These are computed by polynomial multiplication as
  * follows.  Let the sum polynomial be P(z).  Define the elementary
  * polynomials P_i(z) = 1 - 2cos(w_i) z^{-1} + z^{-2}, for 1<=i<=
  * 5, where {w_i} are the LSP frequencies corresponding to the sum
  * polynomial.  Then P(z) = P_1(z)P_2(z)...P_5(z).  Similarly
  * the difference polynomial Q(z) = Q_1(z)Q_2(z)...Q_5(z).
  */
 /*
  * Initialize the arrays with the coefficients of the product
  * P_1(z)P_2(z) and Q_1(z)Q_2(z).  Scale by 1/8.
  */
    P[0] = (Word32) 0x10000000L ;
    P[1] = L_mult( Lsp[0], (Word16) 0x2000 ) ;
    P[1] = L_mac( P[1], Lsp[2], (Word16) 0x2000 ) ;
    P[2] = L_mult( Lsp[0], Lsp[2] ) ;
    P[2] = L_shr( P[2], (Word16) 1 ) ;
    P[2] = L_add( P[2], (Word32) 0x20000000L ) ;
    Q[0] = (Word32) 0x10000000L ;
    Q[1] = L_mult( Lsp[1], (Word16) 0x2000 ) ;
    Q[1] = L_mac( Q[1], Lsp[3], (Word16) 0x2000 ) ;
    Q[2] = L_mult( Lsp[1], Lsp[3] ) ;
    Q[2] = L_shr( Q[2], (Word16) 1 ) ;
    Q[2] = L_add( Q[2], (Word32) 0x20000000L ) ;
 /*
  * Compute the intermediate polynomials P_1(z)P_2(z)...P_i(z) and
  * Q_1(z)Q_2(z)...Q_i(z), for i = 2, 3, 4.  Each intermediate
  * polynomial is symmetric, so only the coefficients up to i+1 need
  * by computed.  Scale by 1/2 each iteration for a total of 1/8.
  */
    for ( i = 2 ; i < LpcOrder/2 ; i ++ ) {
        /* Compute coefficient (i+1) */
        Acc0 = P[i] ;
        Acc0 = L_mls( Acc0, Lsp[2*i+0] ) ;
        Acc0 = L_add( Acc0, P[i-1] ) ;
        P[i+1] = Acc0 ;
        Acc1 = Q[i] ;
        Acc1 = L_mls( Acc1, Lsp[2*i+1] ) ;
        Acc1 = L_add( Acc1, Q[i-1] ) ;
        Q[i+1] = Acc1 ;
        /* Compute coefficients i, i-1, ..., 2 */
        for ( j = i ; j >= 2 ; j -- ) {
            Acc0 = P[j-1] ;
            Acc0 = L_mls( Acc0, Lsp[2*i+0] ) ;
            Acc0 = L_add( Acc0, L_shr(P[j], (Word16) 1 ) ) ;
            Acc0 = L_add( Acc0, L_shr(P[j-2], (Word16) 1 ) ) ;
            P[j] = Acc0 ;
            Acc1 = Q[j-1] ;
            Acc1 = L_mls( Acc1, Lsp[2*i+1] ) ;
            Acc1 = L_add( Acc1, L_shr(Q[j], (Word16) 1 ) ) ;
            Acc1 = L_add( Acc1, L_shr(Q[j-2], (Word16) 1 ) ) ;
            Q[j] = Acc1 ;
        }
        /* Compute coefficients 1, 0 */
        P[0] = L_shr( P[0], (Word16) 1 ) ;
        Q[0] = L_shr( Q[0], (Word16) 1 ) ;
        Acc0 = L_deposit_h( Lsp[2*i+0] ) ;
        Acc0 = L_shr( Acc0, (Word16) i ) ;
        Acc0 = L_add( Acc0, P[1] ) ;
        Acc0 = L_shr( Acc0, (Word16) 1 ) ;
        P[1] = Acc0 ;
        Acc1 = L_deposit_h( Lsp[2*i+1] ) ;
        Acc1 = L_shr( Acc1, (Word16) i ) ;
        Acc1 = L_add( Acc1, Q[1] ) ;
        Acc1 = L_shr( Acc1, (Word16) 1 ) ;
        Q[1] = Acc1 ;
    }
 /*
  * Convert the sum and difference polynomials to LPC coefficients
  * The LPC polynomial is the sum of the sum and difference
  * polynomials with the real zeros factored in: A(z) = 1/2 {P(z) (1
  * + z^{-1}) + Q(z) (1 - z^{-1})}.  The LPC coefficients are scaled
  * here by 16; the overall scale factor for the LPC coefficients
  * returned by this function is therefore 1/4.
  */
    for ( i = 0 ; i < LpcOrder/2 ; i ++ ) {
        Acc0 = P[i] ;
        Acc0 = L_add( Acc0, P[i+1] ) ;
        Acc0 = L_sub( Acc0, Q[i] ) ;
        Acc0 = L_add( Acc0, Q[i+1] ) ;
        Acc0 = L_shl( Acc0, (Word16) 3 ) ;
        Lsp[i] = negate( round( Acc0 ) ) ;
        Acc1 = P[i] ;
        Acc1 = L_add( Acc1, P[i+1] ) ;
        Acc1 = L_add( Acc1, Q[i] ) ;
        Acc1 = L_sub( Acc1, Q[i+1] ) ;
        Acc1 = L_shl( Acc1, (Word16) 3 ) ;
        Lsp[LpcOrder-1-i] = negate( round( Acc1 ) ) ;
    }
}

						