// file: $ECE_8993/util/calculate_lpc/v1.0/clpc_comp_6.cc
//
// computes the lp predictor coefficients
//

// system include files
//
#include <string.h>

// isip include files
//
#include <integral.h>
#include <integral_constants.h>

// local include files
//
#include "calculate_lpc.h"

#include "calculate_lpc_constants.h"

// function: compute_lp_coeffs_cc
//
// arguments:
//
//  float_4* pred_coeffs: (output) predictor coefficients
//  int_4 lp_order: (input) order of the lp model
//  float_4* autocor_coeffs: (input) autocorrelation coefficients
//  float_4& gain: (output) gain of the lp model
//  
// return value: a logical_1 indicating status
//
// this function computes the coefficients for the lp model using the
// levinson-durbin recursion
//
logical_1 compute_lp_coeffs_cc(float_4* pred_coeffs_a, int_4 lp_order_a,
			       float_4* autocor_coeffs_a, float_4& gain_a) {

  // allocate temporary memory for the reflection coefficients and
  // the intermediate predictor coefficients
  //
  float_4* refl_coeffs = new float_4[lp_order_a + 1];
  float_4* temp_pred_coeffs = new float_4[lp_order_a + 1];

  // initialize the each of the coefficients
  //
  memset((void*)pred_coeffs_a, 0, (lp_order_a + 1) * sizeof(float_4));
  memset((void*)temp_pred_coeffs, 0, (lp_order_a + 1) * sizeof(float_4));
  memset((void*)refl_coeffs, 0, (lp_order_a + 1) * sizeof(float_4));
  
  // initialize the error term
  //
  float_4 lp_error = autocor_coeffs_a[0];

  // this is the core of the levinson-durbin recursion
  //
  for (int_4 iter = (int_4)1; iter <= lp_order_a; iter++) {

    // initialize the reflection coefficient
    //
    refl_coeffs[iter] = (float_4)0.0;

    // update all of the reflection coefficients
    //
    for (int_4 j = (int_4)1; j < iter; j++) {

      refl_coeffs[iter] -= temp_pred_coeffs[j] * autocor_coeffs_a[iter - j];
    }

    // add the autocorrelation coefficient to the reflection coeff
    // and divide by the previous error term
    //
    refl_coeffs[iter] += autocor_coeffs_a[iter];
    refl_coeffs[iter] /= lp_error;

    // compute the new set of predictor coefficients
    //
    pred_coeffs_a[iter] = refl_coeffs[iter];

    // update the predictor coefficients
    //
    for (int_4 j = (int_4)1; j < iter; j++) {

      pred_coeffs_a[j] = pred_coeffs_a[j] - (refl_coeffs[iter] *
					     temp_pred_coeffs[iter - j]);
    }

    // compute the new error term
    //
    lp_error = lp_error * ((float_4)1.0 - (refl_coeffs[iter] *
					   refl_coeffs[iter]));

    // update the temporary predictor coeffs
    //
    for (int_4 j = 0; j <= iter; j++) {
      temp_pred_coeffs[j] = pred_coeffs_a[j];
    }
  }

  // set the gain value
  //
  gain_a = lp_error;
  
  // delete memory
  //
  delete [] refl_coeffs;
  delete [] temp_pred_coeffs;

  // exit gracefully
  //
  return ISIP_TRUE;
}