/******************************************************************************
*                                                                             *
*       Copyright (C) 1992-1995 Tony Robinson                                 *
*                                                                             *
*       See the file LICENSE for conditions on distribution and usage         *
*                                                                             *
******************************************************************************/

/*
 * $Id: lpc.c-new,v 1.1 2004/03/17 22:17:54 robux4 Exp $
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "../include/shorten.h"

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#define log2(x)                         ( log (x) * (1./M_LN2) )
#define E_BITS_PER_COEF         ( 2 + LPCQUANT )                // watch out, these are all 0 .. order inclusive arrays


int                                                             // best prediction order model
wav2lpc ( slong*  buf,                  // Samples
                  int     nbuf,                 // Number of samples
                  slong   offset,               //
                  int*    qlpc,                 // quantized prediction coefficients
                  int     nlpc,                 // max. prediction order
                  int     version,              // some strange version information
          float*  psigbit,              // expected number of bits per original signal sample
          float*  presbit )             // expected number of bits per residual signal sample
{
        static double*  fbuf  = NULL;
        static int      nflpc = 0;
        static int      nfbuf = 0;
        int             i;
        int             j;
        int             bestnbit;
        int             bestnlpc;
        double          e;
        double          bestesize;
        double          ci;
        double          esize;
        double          acf [MAX_LPC_ORDER + 1];
        double          ref [MAX_LPC_ORDER + 1];
        double          lpc [MAX_LPC_ORDER + 1];
        double          tmp [MAX_LPC_ORDER + 1];
        double          escale = 0.5 * M_LN2 * M_LN2 / nbuf;
        double          sum;

        if ( nlpc >= nbuf )                                                             // if necessary, limit the LPC order to the number of samples available
                nlpc = nbuf - 1;

        if ( nlpc > nflpc  ||  nbuf > nfbuf ) {                 // grab some space for a 'zero mean' buffer of floats if needed
                if ( fbuf != NULL )
                        free ( fbuf - nflpc );
                fbuf  = nlpc + ((double*) pmalloc ( (nlpc+nbuf) * sizeof (*fbuf) ));
                nfbuf = nbuf;
                nflpc = nlpc;
        }

        e = 0.;
        for ( j = 0; j < nbuf; j++ ) {                                  // zero mean signal and compute energy
                sum = fbuf [j] -= offset;
                e  += sum * sum;
        }

        esize     = e > 0.  ?  0.5 * log2 (escale * e)  :  0.;
        *psigbit  = esize;                                                              // return the expected number of bits per original signal sample

        acf [0]   = e;                                                                  // store the best values so far (the zeroth order predictor)
        bestnlpc  = 0;
        bestnbit  = (int) floor (nbuf * esize);
        bestesize = esize;

        // check all linear predictors up to and including length nlpc if version is 2 or greater, just check two more than bestnlpc
        // AJR:  8 May 1996: the code used to read "version < 12", it should read "bestnlpc + 2" but
        //       changed to "bestnlpc + 3" to be more conservative (and more in line with old behaviour

        for ( i = 1; i <= nlpc  &&  e > 0.  &&  (version < 2  ||  i <= bestnlpc + 3); i++ ) {

                sum = 0.;
                for ( j = i; j < nbuf; j++ )                                                                                    // compute the jth autocorrelation coefficient
                        sum += fbuf [j] * fbuf [j-i];
                acf [i] = sum;

                ci = 0.;                                                                                                                                // compute the reflection and LP coeffients for order j predictor
                for ( j = 1; j < i; j++ )
                        ci += lpc [j] * acf [i-j];
                lpc [i] = ref [i] = ci = (acf [i] - ci) / e;
                for ( j = 1; j < i; j++ )
                        tmp [j] = lpc [j] - ci * lpc [i-j];
                for ( j = 1; j < i; j++ )
                        lpc [j] = tmp [j];

                e    *= 1 - ci*ci;                                                                                                              // compute the new energy in the prediction residual
                esize = e > 0.  ?  0.5 * log2 (escale * e)  :  0.;

                if ( nbuf * esize + i * E_BITS_PER_COEF < bestnbit ) {                                  // store this model if it is the best so far
                        bestnlpc  = i;                                                                                                          // store best model order
                        bestnbit  = (int) floor ( nbuf * esize + i * E_BITS_PER_COEF );
                        bestesize = esize;

                        for ( j = 0; j < bestnlpc; j++ )                                                                        // store the quantised LP coefficients
                                qlpc [j] = (int) floor ( lpc [j+1] * (1 << LPCQUANT) + 0.5 );
                }
        }

        *presbit = bestesize;                                           // return the expected number of bits per residual signal sample
        return bestnlpc;                                                        // return the best model order
}