/*
**
** File:    lpc2.c
**
** Description: Functions that implement linear predictive coding
**      (LPC) operations.
**
** Functions:
**
**  Computing LPC coefficients:
**
**      Comp_Lpc()
**      Durbin()
**
**  Perceptual noise weighting:
**
**      Wght_Lpc()
**      Error_Wght()
**
**  Computing combined impulse response:
**
**      Comp_Ir()
**
**  Computing ringing response:
**
**      Sub_Ring()
**      Upd_Ring()
**
**  Synthesizing speech:
**
**      Synt()
**      Spf()
*/

/*
    ITU-T G.723.1 Floating Point Speech Coder ANSI C Source Code.  Version 5.1F

    Original fixed-point code copyright (c) 1995,
    AudioCodes, DSP Group, France Telecom, Universite de Sherbrooke.
    All rights reserved.

    Floating-point code copyright (c) 1995,
    Intel Corporation and France Telecom (CNET).
    All rights reserved.
*/

#include <stdio.h>
#include <math.h>
#include <float.h>

#include "typedef2.h"
#include "cst2.h"
#include "tab2.h"
#include "lbccode2.h"
#include "coder2.h"
#include "decod2.h"
#include "util2.h"
#include "lpc2.h"
#include "codcng2.h"

/*
**
** Function:        Comp_Lpc()
**
** Description:     Computes the tenth-order LPC filters for an
**          entire frame.  For each subframe, a
**          Hamming-windowed block of 180 samples,
**          centered around the subframe, is used to
**          compute eleven autocorrelation coefficients.
**          The Levinson-Durbin algorithm then generates
**          the LPC coefficients.  This function requires
**          a look-ahead of one subframe, and hence
**          introduces a 7.5 ms encoding delay.
**
** Links to text:   Section 2.4
**
** Arguments:
**
**  FLOAT  *UnqLpc      Empty Buffer
**  FLOAT  PrevDat[]    Previous 2 subframes of samples (120 words)
**  FLOAT  DataBuff[]   Current frame of samples (240 words)
**
** Outputs:
**
**
**  FLOAT  UnqLpc[]     LPC coefficients for entire frame (40 words)
**
** Return value:    None
**
*/
void  Comp_Lpc(FLOAT *UnqLpc, FLOAT *PrevDat, FLOAT *DataBuff)
{
    int   i,j,k;

    FLOAT  Dpnt[Frame+LpcFrame-SubFrLen];
    FLOAT  Vect[LpcFrame];
    FLOAT  Acf_sf[LpcOrderP1*SubFrames];
    FLOAT  *curAcf;
    FLOAT  Pk2;

    /*
     * Generate a buffer of 360 samples.  This consists of 120 samples
     * from the previous frame and 240 samples from the current frame.
     */

    for (i=0; i < LpcFrame-SubFrLen; i++)
        Dpnt[i] = PrevDat[i];
    for (i=0; i < Frame; i++)
        Dpnt[i+LpcFrame-SubFrLen] = DataBuff[i];

    /*
     * Repeat for all subframes
     */

    curAcf = Acf_sf;
    for (k=0; k < SubFrames; k++)
    {
        /*  Apply the Hamming window */

        for (i = 0; i < LpcFrame; i++)
            Vect[i] = Dpnt[k*SubFrLen+i]*HammingWindowTable[i];

        /*  Compute the autocorrelation coefficients  */

        curAcf[0] = DotProd(Vect, Vect, LpcFrame)/(LpcFrame*LpcFrame);

        /*  Do Ridge regression  */

        curAcf[0] *= ((FLOAT)1025.0/(FLOAT)1024.0);
        if (curAcf[0] == (FLOAT)0.0) {
            for (i = 1; i <= LpcOrder; i++)
                curAcf[i] = (FLOAT) 0.0;
        }
        else {
            for (i = 1; i <= LpcOrder; i++)
                curAcf[i] = DotProd(Vect, &Vect[i], LpcFrame-i)
                            / (LpcFrame*LpcFrame) * BinomialWindowTable[i-1];
        }

        /*
         * Apply the Levinson-Durbin algorithm to generate the LPC
         * coefficients
         */

        Durbin(&UnqLpc[k*LpcOrder], &curAcf[1], curAcf[0], &Pk2);
        CodStat.SinDet <<= 1;
        if (Pk2 > (FLOAT) 0.95) {
            CodStat.SinDet++;
        }
        curAcf += LpcOrderP1;
    }

    /* Update sine detector */
    CodStat.SinDet &= 0x7fff ;

    j = CodStat.SinDet ;
    k = 0 ;
    for ( i = 0 ; i < 15 ; i ++ ) {
        k += j & 1 ;
        j >>= 1 ;
    }
    if ( k >= 14 )
        CodStat.SinDet |= 0x8000 ;

    Update_Acf(Acf_sf);
}


/*
**
** Function:        Durbin()
**
** Description:     Implements the Levinson-Durbin algorithm for a
**          subframe.  The Levinson-Durbin algorithm
**          recursively computes the minimum mean-squared
**          error (MMSE) linear prediction filter based on the
**          estimated autocorrelation coefficients.
**
** Links to text:   Section 2.4
**
** Arguments:
**
**  FLOAT *Lpc Empty buffer
**  FLOAT Corr[]   First- through tenth-order autocorrelations (10 words)
**  FLOAT Err  Zeroth-order autocorrelation, or energy
**
** Outputs:
**
**  FLOAT Lpc[]    LPC coefficients (10 words)
**
** Return value:    None
**
*/
FLOAT  Durbin(FLOAT *Lpc, FLOAT *Corr, FLOAT Err, FLOAT *Pk2)
{
    int    i,j;
    FLOAT  Temp[LpcOrder];
    FLOAT  Pk,Tmp0;

    /*  Initialize the LPC vector  */

    for (i=0; i < LpcOrder; i++)
        Lpc[i] = (FLOAT)0.0;

    /*  Compute the partial correlation (parcor) coefficient  */

    for (i=0; i < LpcOrder; i++)
    {
        Tmp0 = Corr[i];
        for (j=0; j<i; j++)
            Tmp0 -= Lpc[j]*Corr[i-j-1];

        if (fabs(Tmp0) >= Err) {
            *Pk2 = (FLOAT)0.99;
            break;
        }

        Lpc[i] = Pk = Tmp0/Err;
        Err -= Tmp0*Pk;

        /*
         * Sine detector
         */
        if ( i == 1 )
            *Pk2 = -Pk;

        for (j=0; j < i; j++)
            Temp[j] = Lpc[j];

        for (j=0; j < i; j++)
            Lpc[j] = Lpc[j] - Pk*Temp[i-j-1];
    }
    return Err;
}


/*
**
** Function:        Wght_Lpc()
**
** Description:     Computes the formant perceptual weighting
**          filter coefficients for a frame.  These
**          coefficients are geometrically scaled versions
**          of the unquantized LPC coefficients.
**
** Links to text:   Section 2.8
**
** Arguments:
**
**  FLOAT  *PerLpc      Empty Buffer
**  FLOAT  UnqLpc[]     Unquantized LPC coefficients (40 words)
**
** Outputs:
**
**
**  FLOAT  PerLpc[]     Perceptual weighting filter coefficients
**              (80 words)
**
** Return value:    None
**
*/
void  Wght_Lpc(FLOAT *PerLpc, FLOAT *UnqLpc)
{
    int  i,j;

    for (i=0; i < SubFrames; i++)
    {
        /*
         * Compute the jth FIR coefficient by multiplying the jth LPC
         * coefficient by (0.9)^j. Compute the jth IIR coefficient by
         * multiplying the jth LPC coefficient by (0.5)^j.
         */

        for (j=0; j < LpcOrder; j++)
        {
            PerLpc[j]          = UnqLpc[j]*PerFiltZeroTable[j];
            PerLpc[j+LpcOrder] = UnqLpc[j]*PerFiltPoleTable[j];
        }
        PerLpc += 2*LpcOrder;
        UnqLpc += LpcOrder;
    }
}

/*
**
** Function:        Error_Wght()
**
** Description:     Implements the formant perceptual weighting
**          filter for a frame. This filter effectively
**          deemphasizes the formant frequencies in the
**          error signal.
**
** Links to text:   Section 2.8
**
** Arguments:
**
**  FLOAT  Dpnt[]       Highpass filtered speech x[n] (240 words)
**  FLOAT  PerLpc[]     Filter coefficients (80 words)
**
** Inputs:
**
**  CodStat.WghtFirDl[] FIR filter memory from previous frame (10 words)
**  CodStat.WghtIirDl[] IIR filter memory from previous frame (10 words)
**
**
** Outputs:
**
**  FLOAT  Dpnt[]       Weighted speech f[n] (240 words)
**
** Return value:    None
**
*/
void  Error_Wght(FLOAT *Dpnt, FLOAT *PerLpc)
{
    int  i,j,k;

    FLOAT Acc0;


    for (k=0; k < SubFrames; k++)
    {
        for (i=0; i < SubFrLen; i++)
        {
            /*  FIR part  */

            Acc0 = *Dpnt - DotProd(PerLpc,CodStat.WghtFirDl,LpcOrder);

            for (j=LpcOrder-1; j > 0; j --)
                CodStat.WghtFirDl[j] = CodStat.WghtFirDl[j-1];

            CodStat.WghtFirDl[0] = *Dpnt;

            /*  IIR part  */

            Acc0 += DotProd(&PerLpc[LpcOrder],CodStat.WghtIirDl,LpcOrder);

            for (j = LpcOrder-1; j > 0; j --)
                CodStat.WghtIirDl[j] = CodStat.WghtIirDl[j-1];

            *Dpnt++ = CodStat.WghtIirDl[0] = Acc0;

        }
        PerLpc += 2*LpcOrder;
    }
}


/*
**
** Function:        Comp_Ir()
**
** Description:     Computes the combined impulse response of the
**          formant perceptual weighting filter, harmonic
**          noise shaping filter, and synthesis filter for
**          a subframe.
**
** Links to text:   Section 2.12
**
** Arguments:
**
**  FLOAT  *ImpResp     Empty Buffer
**  FLOAT  QntLpc[]     Quantized LPC coefficients (10 words)
**  FLOAT  PerLpc[]     Perceptual filter coefficients (20 words)
**  PWDEF  Pw           Harmonic noise shaping filter parameters
**
** Outputs:
**
**  FLOAT  ImpResp[]    Combined impulse response (60 words)
**
** Return value:    None
**
*/
void  Comp_Ir(FLOAT *ImpResp, FLOAT *QntLpc, FLOAT *PerLpc, PWDEF Pw)
{
    int    i,j;

    FLOAT  FirDl[LpcOrder];
    FLOAT  IirDl[LpcOrder];
    FLOAT  Temp[PitchMax+SubFrLen];
    FLOAT  Acc0,Acc1;

    /*