// file: $isip/class/algo/Prediction/pred_05.cc
// version: $Id: pred_05.cc 8037 2002-04-01 20:04:33Z gao $
//


// isip include files
//
#include "Prediction.h"

// method: apply
//
// arguments:
//  Vector<AlgorithmData>& output: (output) output data
//  const Vector< CircularBuffer<AlgorithmData> >& input: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes the predictor coefficients for input data
//
bool8 Prediction::apply(Vector<AlgorithmData>& output_a,
			  const Vector< CircularBuffer<AlgorithmData> >&
			  input_a) {

  // determine the number of input channels and force the output to be
  // that number
  //
  int32 len = input_a.length();
  output_a.setLength(len);
  bool8 res = true;

  // start the debugging output
  //
  displayStart(this);
  
  // loop over the channels and call the compute method
  //
  for (int32 c = 0; c < len; c++) {
    
    // display the channel number
    //
    displayChannel(c);
    
    // branch on input type - covariance takes matrix as input
    //
    if (input_a(c)(0).getCoefType() == AlgorithmData::COVARIANCE) {
      
      // call AlgorithmData::makeVectorFloat to force the output vector for
      // this channel to be a VectorFloat.
      // call AlgorithmData::getMatrixFloat on the input for this channel to
      // check that the input is already a MatrixFloat and return that
      // vector.
      //      
      res &= compute(output_a(c).makeVectorFloat(),
		     input_a(c)(0).getMatrixFloat(),
		     AlgorithmData::COVARIANCE, c);
    }
  
    else {

      // call AlgorithmData::makeVectorFloat to force the output vector for
      // this channel to be a VectorFloat.
      // call AlgorithmData::getVectorFloat on the input for this channel to
      // check that the input is already a VectorFloat and return that
      // vector.
      // 
      res &= compute(output_a(c).makeVectorFloat(),
		     input_a(c)(0).getVectorFloat(),
		     input_a(c)(0).getCoefType(), c);
    }

    // set the coefficient type for the output
    //    
    output_a(c).setCoefType(AlgorithmData::PREDICTION);    
  }

  // finish the debugging output
  //
  displayFinish(this);

  // exit gracefully
  //
  return res;  
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) predictor or reflection coefficents
//  const VectorFloat& input: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel number
//
// return: a bool8 value indicating status
//
// this method computes the predictor coefficents for the given input data
//
bool8 Prediction::compute(VectorFloat& output_a,
			    const VectorFloat& input_a,
			    AlgorithmData::COEF_TYPE input_coef_type_a,
			    int32 channel_index_a) {
  
  // declare local variables
  //
  float err_energy;
  bool8 status = false;

  // check order
  //
  if (order_d < (int32)1) {
    output_a.setLength(0);
    return Error::handle(name(), L"compute", ERR_ORDER,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // branch on algorithm
  //  Algorithm: AUTOCORRELATION
  //
  if (algorithm_d == AUTOCORRELATION) {
    if (implementation_d == DURBIN) {
      if (input_coef_type_a == AlgorithmData::CORRELATION) {
	status = computeAutoDurbin(output_a, err_energy, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }

  // Algorithm: LATTICE
  //
  else if (algorithm_d == LATTICE) {
    if (implementation_d == BURG) {
      if (input_coef_type_a == AlgorithmData::SIGNAL) {
	status = computeLatticeBurg(output_a, err_energy, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }
  
  // Algorithm: STEP_DOWN
  //
  else if (algorithm_d == REFLECTION) {
    if (implementation_d == STEP_DOWN) {
      if (input_coef_type_a == AlgorithmData::REFLECTION) {
	status = computeReflectionStepDown(output_a, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }

  // Algorithm: LATTICE
  //
  else if (algorithm_d == LOG_AREA_RATIO) {
    if (implementation_d == KELLY_LOCHBAUM) {
      if (input_coef_type_a == AlgorithmData::LOG_AREA_RATIO) {
	status = computeLogAreaKellyLochbaum(output_a, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }
  
  // else: error unknown algorithm
  //
  else {
    return Error::handle(name(), L"compute", ERR_UNKALG,
			 __FILE__, __LINE__);
  }      
    
  // possibly display the data
  //
  display(output_a, input_a, name());
  
  // exit gracefully
  //
  return status;
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) predictor or reflection coefficents
//  const MatrixFloat& input: (input) input data (covariance matrix)
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel number
//
// return: a bool8 value indicating status
//
// this method computes the predictor coefficents from Covariance
// input type
//
bool8 Prediction::compute(VectorFloat& output_a,
			    const MatrixFloat& input_a,
			    AlgorithmData::COEF_TYPE input_coef_type_a,
			    int32 channel_index_a) {

  // declare local variables
  //
  float err_energy;
  bool8 status = false;
  
  // check the order
  //
  if (order_d < (int32)1) {
    output_a.setLength(0);
    return Error::handle(name(), L"compute", ERR_ORDER,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // branch on algorithm
  //
  if (algorithm_d == COVARIANCE) {
    if (implementation_d == CHOLESKY) {
      if (input_coef_type_a == AlgorithmData::COVARIANCE) {
	status = computeCovarCholesky(output_a, err_energy, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }

  // else: error unknown algorithm
  //
  else {
    return Error::handle(name(), L"compute", ERR_UNKALG,
			 __FILE__, __LINE__);
  }      

  // possibly display the data
  //
  display(output_a, input_a, name());
  
  // exit gracefully
  //  
  return status;
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) predictor or reflection coefficents
//  float& err_energy: (output) error energy
//  const VectorFloat& input: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel number
//
// return: a bool8 value indicating status
//
// this method computes the predictor coefficents. it also outputs the
// error energy.
//
bool8 Prediction::compute(VectorFloat& output_a, float& err_energy_a,
			    const VectorFloat& input_a,
			    AlgorithmData::COEF_TYPE input_coef_type_a,
			    int32 channel_index_a) {

  // initialize err_energy a value since in some mode there is no
  // err_energy output
  //
  err_energy_a = 0.0;
  bool8 status = false;
  
  // check the order
  //  
  if (order_d < (int32)1) {
    output_a.setLength(0);
    return Error::handle(name(), L"compute", ERR_ORDER,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // branch on algorithm
  //  Algorithm: AUTOCORRELATION
  //
  if (algorithm_d == AUTOCORRELATION) {
    if (implementation_d == DURBIN) {
      if (input_coef_type_a == AlgorithmData::CORRELATION) {
	status = computeAutoDurbin(output_a, err_energy_a, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }

  // Algorithm: LATTICE
  //
  else if (algorithm_d == LATTICE) {
    if (implementation_d == BURG) {
      if (input_coef_type_a == AlgorithmData::SIGNAL) {
	status = computeLatticeBurg(output_a, err_energy_a, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }
  
  // Algorithm: STEP_DOWN
  //
  else if (algorithm_d == REFLECTION) {
    if (implementation_d == STEP_DOWN) {
      if (input_coef_type_a == AlgorithmData::REFLECTION) {
	status = computeReflectionStepDown(output_a, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }

  // else: error unknown algorithm
  //
  else {
    return Error::handle(name(), L"compute", ERR_UNKALG,
			 __FILE__, __LINE__);
  }      

  // possibly display the data
  //
  display(output_a, input_a, name());
    
  // exit gracefully
  //  
  return status;
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) predictor or reflection coefficents
//  float& err_energy: (output) error energy
//  const MatrixFloat& input: (input) input data (covariance matrix)
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel number
//
// return: a bool8 value indicating status
//
// this method computes the prediction coefficents from Covariance
// input type. it also outputs the error energy.
//
bool8 Prediction::compute(VectorFloat& output_a, float& err_energy_a,
			    const MatrixFloat& input_a,
			    AlgorithmData::COEF_TYPE input_coef_type_a,
			    int32 channel_index_a) {

  // initialize a value for err_energy since in some mode there is no
  // err_energy output
  //
  err_energy_a = 0.0;
  bool8 status = false;
  
  // check the order
  //   
  if (order_d < (int32)1) {
    output_a.setLength(0);
    return Error::handle(name(), L"compute", ERR_ORDER,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // branch on algorithm
  //
  if (algorithm_d == COVARIANCE) {
    if (implementation_d == CHOLESKY) {
      if (input_coef_type_a == AlgorithmData::COVARIANCE) {
	status = computeCovarCholesky(output_a, err_energy_a, input_a);
      }

      // else: error unknown coefficient type
      //
      else {
	return Error::handle(name(), L"compute", ERR_UNCTYP,
			     __FILE__, __LINE__);
      }
    }

    // else: error unknown implementation
    //    
    else {
      return Error::handle(name(), L"compute", ERR_UNKIMP,
			   __FILE__, __LINE__);
    }      
  }

  // else: error unknown algorithm
  //
  else {
    return Error::handle(name(), L"compute", ERR_UNKALG,
			 __FILE__, __LINE__);
  }      

  // possibly display the data
  //
  display(output_a, input_a, name());

  // exit gracefully
  //  
  return status;
}

// method: computeAutoDurbin
//
// arguments:
//  VectorFloat& pred_coef: (output) prediction coefficients
//  float& err_energy: (output) error energy  
//  const VectorFloat& auto: (input) autocorrelation coefficients
//
// return: a bool8 value indicating status
//
// this method computes the prediction coefficients from the
// autocorrelation coefficents. it also outputs the error energy.
//
// Reference:
//  J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
//  Springer-Verlag Berlin Heidelberg, New York, USA, pp. 217-219,
//  1980.
//
bool8 Prediction::computeAutoDurbin(VectorFloat& pred_coef_a,
				      float& err_energy_a,
				      const VectorFloat& auto_a) {

  // apply dynamic range threshold
  //
  float old_auto_0 = auto_a(0);
  if (dyn_range_d < (float)0) {
    float64 scale_factor =
      pow(10.0, (float)dyn_range_d / Integral::DB_POW_SCALE_FACTOR);
    const_cast<VectorFloat&>(auto_a)(0) *= (1 + scale_factor);
  }  
  else {
    return Error::handle(name(), L"computeAutoDurbin", ERR_DYNRANGE,
			 __FILE__, __LINE__);
  }  
  
  // declare local variables
  //
  int32 i, j;                                 // loop counters
  int32 lp_order = auto_a.length() - 1;       // analysis order
  int32 j_bckwd;                              // a backwards index
  int32 middle_index;                         // inner loop limit
  float64 sum;                                // accumulator
  float64 rc_reg;                             // register
  float tmp;
  
  // create space
  //
  if (!pred_coef_a.setLength(lp_order + 1)) {
    return Error::handle(name(), L"computeAutoDurbin", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }
  
  // initialization
  //
  pred_coef_a(0) = 1.0;
  err_energy_a = auto_a(0);

  // if the energy of the signal is zero we set all the predictor
  // coefficients to zero and exit
  //
  if (err_energy_a == (float)0) {
    for (int i = 0; i < pred_coef_a.length(); i++) {
      pred_coef_a(i) = 0;
    }
    pred_coef_a(0) = 1;
    return true;
  }
  
  // do the first step manually
  //
  if (err_energy_a == (float)0) {
    return Error::handle(name(), L"computeAutoDurbin", Error::ARG, __FILE__, __LINE__);
  }
  
  sum = -(float64)auto_a(1) / err_energy_a;
  pred_coef_a(1) = sum;
  tmp = 1 - sum * sum;
  err_energy_a *= tmp;

  // recursion
  //
  for (i = 2; i <= lp_order; i++) {

    // compute the next reflection coefficient
    //
    sum = 0;
    for (j = 1; j < i; j++) {
      sum += pred_coef_a(j) * auto_a(i - j);
    }

    if (err_energy_a == (float)0) {
      return Error::handle(name(), L"computeAutoDurbin", Error::ARG, __FILE__, __LINE__);
    }    
    rc_reg = - (auto_a(i) + (float)sum) / err_energy_a;

    // compute the new prediction coefficients
    //  note: the algorithm used here minimizes the storage.
    //
    // Reference:
    //  J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
    //  Springer-Verlag Berlin Heidelberg, New York, USA, pp. 219,
    //  1980.
    //
    pred_coef_a(i) = rc_reg;
    j_bckwd = i - 1;
    middle_index = i / 2;
    for (j = 1; j <= middle_index; j++) {
      sum = pred_coef_a(j_bckwd) + (float)rc_reg* pred_coef_a(j); 
      pred_coef_a(j) = pred_coef_a(j) + (float)rc_reg * pred_coef_a(i - j);
      pred_coef_a(j_bckwd) = sum;
      j_bckwd--;
    }

    // compute new error
    //
    tmp = 1.0 - rc_reg * rc_reg;
    err_energy_a *= tmp;
  }
  
  // resume the auto_a(0)
  //
  const_cast<VectorFloat&>(auto_a)(0) = old_auto_0;
  
  // exit gracefully
  //
  return true;
}

// method: computeCovarCholesky
//
// arguments:
//  VectorFloat& pred_coef: (output) prediction coefficients
//  float& err_energy: (output) error energy  
//  const MatrixFloat& cor: (input) covariance coefficients
//
// return: a bool8 value indicating status
//
// this method computes the prediction coefficients using the
// covariance analysis method. it also outputs the error energy.
//
// Reference:
//  J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
//  Springer-Verlag Berlin Heidelberg, New York, USA, pp. 51-55, 1980.
//
bool8 Prediction::computeCovarCholesky(VectorFloat& pred_coef_a,
					 float& err_energy_a,
					 const MatrixFloat& cor_a) {

  // apply dynamic range threshold
  //
  VectorFloat temp_vector;
  VectorFloat refl_coef;
  float value;
  
  // add the dynamic range threshold
  //
  if (dyn_range_d < (float)0.0) {
    float64 scale_factor =
      pow(10.0, ((float)dyn_range_d / Integral::DB_POW_SCALE_FACTOR));
    int32 num_rows = cor_a.getNumRows();
    temp_vector.setLength(num_rows);
    for (int32 i = 0; i < num_rows; i++) {
      value = cor_a.getValue(i, i);
      temp_vector(i) = value;
      value *= 1 + scale_factor;
      const_cast<MatrixFloat&>(cor_a).setValue(i, i, value);
    }
  }
  else {
    return Error::handle(name(), L"computeCovarCholesky",
			 ERR_DYNRANGE, __FILE__, __LINE__);
  }
  
  // compute the order
  //
  int32 lp_order = order_d;

  // error checking on the order size
  //
  if (lp_order > cor_a.getNumRows() - 1) {
    return Error::handle(name(), L"computeCovarCholesky", Error::ARG, __FILE__, __LINE__);
  }
  
  // declare some scratch variables
  //
  MatrixFloat b(lp_order + (int32)1, lp_order + (int32)1);
  VectorFloat beta(lp_order + 1);
    
  // create space
  //
  if (!pred_coef_a.setLength(lp_order + 1)) {
    return Error::handle(name(), L"computeCovarCholesky", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }
  if (!refl_coef.setLength(lp_order)) {
    return Error::handle(name(), L"computeCovarCholesky", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // declare local variables
  //
  int32 m  = lp_order;
  int32 mf = lp_order;
  int32 minc;
  int32 ip;
  float64 gam;
  int32 j;
  int32 jp;
  float64 s;
  float tmp;
  
  // Note: the equation numbers indicated here are from the COVAR
  // algorithm in following refence:
  //
  //  J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
  //  Springer-Verlag Berlin Heidelberg, New York, USA, pp. 52, 1980.
  
  // initialization
  //
  b.setValue(1, 1, 1.0);                            // Eq. 3.40
  tmp = cor_a.getValue(1 - 1, 1 - 1);
  err_energy_a = tmp;                               // Eq. 3.41
  beta(1) = cor_a.getValue(2 - 1, 2 - 1);           // Eq. 3.41
  refl_coef(0) = -(float64)cor_a.getValue(2 - 1, 1 - 1) / 
                    cor_a.getValue(2 - 1, 2 - 1);   // Eq. 3.42
  pred_coef_a(0) = 1;                               // Eq. 3.43
  pred_coef_a(1) = refl_coef((int32)0);              // Eq. 3.43
  tmp = tmp - refl_coef(0) * refl_coef(0) * beta(1);
  err_energy_a = tmp;                               // Eq. 3.44

  // iteration
  //
  for (minc = 2; minc <= mf; minc++) {              // do 400 ...
    m = minc - 1;
    b.setValue(minc, minc, 1);                      // Eq. 3.51a

    for (ip = 1; ip <= m; ip++) {                   // do 200 ...
      if ((float64)beta(ip) < 0) {                   // if 600, 200, 130
	return Error::handle(name(), L"computeCovarCholesky", ERR_BETA,
			     __FILE__, __LINE__, Error::WARNING);
      }

      // Eq. 3.50 and 3.51b
      //
      else if ((float64)beta(ip) > 0) {
        gam = 0;
        for (j = 1; j <= ip; j++) {
          gam += cor_a.getValue(m + 2 - 1, j + 1 - 1) * b.getValue(ip, j);
	}
        gam /= (float64)beta(ip);
        for (jp = 1; jp <= ip; jp++) {
	  value = b.getValue(minc, jp) - gam * b.getValue(ip, jp);
	  b.setValue(minc, jp, value);
	}
      }

      // if beta(m) goes to zero, output zero predictor
      // coefficients. this is not part of Markel and Gray's algorithm.
      //
      else if ((float64)beta(m) == 0) {
	err_energy_a = 0.0;
	for (int i = 0; i < pred_coef_a.length(); i++) {
	  pred_coef_a(i) = 0;
	}
	pred_coef_a(0) = 1;      
	break;
      }
    }

    // Eq. 3.52
    //
    beta(minc) = 0;
    for (j = 1; j <= minc; j++) {
      beta(minc) += cor_a.getValue(m + 2 - 1, j + 1 - 1) * b.getValue(minc, j);
    }
    m++;
    
    if ((float64)beta(m) < 0) {                      // If 600,360,260
      return Error::handle(name(), L"computeCovarCholesky", ERR_BETA,
			   __FILE__, __LINE__, Error::WARNING);
    }

    // Eq. 3.55
    //
    else if ((float64)beta(m) > 0) {
      s = 0;
      for (ip = 1; ip <= m; ip++) {
        s += cor_a.getValue(m + 1 - 1, ip - 1) * pred_coef_a(ip - 1);
      }

      // Eq. 3.56
      //
      refl_coef(m - 1) = -s / beta(m);
      for (ip = 2; ip <= m; ip++) {
        pred_coef_a(ip - 1) += refl_coef(m - 1) * b.getValue(m,ip - 1);
      }
      pred_coef_a(m) = refl_coef(m - 1);
    }

    // if beta(m) goes to zero, output zero predictor
    // coefficients. this is not part of Markel and Gray's algorithm.
    //
    else if ((float64)beta(m) == 0) {
      err_energy_a = 0.0;
      for (int i = 0; i < pred_coef_a.length(); i++) {
	pred_coef_a(i) = 0;
      }
      pred_coef_a(0) = 1;      
      break;
    }

    // Eq. 3.57
    //
    tmp -= refl_coef(m - 1) * refl_coef(m - 1) * beta(m);
    err_energy_a = tmp;

    // from page 56:
    //
    //  "If any of the norm squares are zero or negative, the Fortran
    //  subroutines terminate the process, since the results are
    //  either in error (in the case of negative norm squares) or may
    //  be exceeding the accuracy of the computer (in the case of a
    //  zero norm square). With floating-point arithmetic, a zero norm
    //  square is highly unlikely. It is more likely that round-off or
    //  truncation errors would result in a very small norm square
    //  rather than one which is zero. A zero norm can theoretically
    //  occur only in the covariance method."
    //
    // since this is the covariance method, we therefore allow the
    // value of zero since it is theoretically possible but error on
    // less than zero.
    //
    if (err_energy_a == 0) {
      break;
    }
    else if (err_energy_a < 0) {
      return Error::handle(name(), L"computeCovarCholesky", ERR_ENERGY, __FILE__,
			   __LINE__, Error::WARNING);
    }
  }

  // resume the values before applying the dynamic threshold
  //
  if (dyn_range_d < (float)0.0) {
    int32 num_rows =  temp_vector.length();
    for (int32 i = 0; i < num_rows; i++) {
      const_cast<MatrixFloat&>(cor_a).setValue(i, i, temp_vector(i));
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: computeLatticeBurg
//
// arguments:
//  VectorFloat& pred_coef: (output) prediction coefficients
//  float& err_energy: (output) error energy  
//  const VectorFloat& input: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes the prediction coefficients using the burg
// analysis method. it also outputs the error energy.
//
bool8 Prediction::computeLatticeBurg(VectorFloat& pred_coef_a,
				       float& err_energy_a,
				       const VectorFloat& input_a) {
  
  // resize output vectors
  //
  int32 pc_order = order_d + (int32)1;

  if (!pred_coef_a.setLength(pc_order)) {
    return Error::handle(name(), L"computeLatticeBurg", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // declare local variables:
  //  Float Registers/Accumulators
  //
  float64 sum, pe, akk, ai, aj; 
  float64 sn, sd;
  akk = 0;
  float tmp;
  
  // Indices In Autocorrelation
  //
  int32 i, k;                                   
  int32 n = input_a.length();

  // find the signal power
  //
  sum = input_a.sumSquare();

  // if the signal power is zero, output zero predictor coefficients
  //
  if (sum == 0) {
    for (int i = 0; i < pred_coef_a.length(); i++) {
      pred_coef_a(i) = 0;      
    }
    pred_coef_a(0) = 1;    
    return true;
  }

  // compute error energy
  //
  if (n == 0) {
    return Error::handle(name(), L"computeLatticeBurg", Error::ARG, __FILE__, __LINE__);
  }  
  err_energy_a = sum / n;

  // declare some internal scratch space:
  //  note that this data is declared static so that it can be shared across
  //  all objects.
  //
  float64 *tmp_r;
  float64 *tmp_rc;
  float64 *tmp_a;
  float64 *tmp_pef;
  float64 *tmp_per;

  // create new space
  //
  tmp_r = new float64[pc_order];
  tmp_rc = new float64[pc_order];
  tmp_a = new float64[pc_order];
  tmp_pef = new float64[input_a.length() + 1];
  tmp_per = new float64[input_a.length() + 1];

  // initialize
  //
  float64 scale_factor = 1.0 +
    pow(10.0, ((float)dyn_range_d / Integral::DB_POW_SCALE_FACTOR));
  pe = sum * scale_factor;
  tmp_r[0] = sum * scale_factor;
  tmp_a[0] = 1.0;

  // main loop
  //
  for (k = 1; k <= order_d; k++) {

    // compute initial prediction errors
    //
    if (k == 1) {
      for (i = 2; i <= n; i++) {
        tmp_pef[i] = (float64)input_a(i - 1);
        tmp_per[i] = (float64)input_a(i - 2);
      }
    }

    // update prediction errors
    //
    else {
      for (i = n; i >= k + 1; i--) {
        tmp_pef[i] = tmp_pef[i] + akk * tmp_per[i];
        tmp_per[i] = tmp_per[i - 1] + akk * tmp_pef[i - 1];
      }
    }

    // compute new reflection coefficient
    //
    sn = 0;
    sd = 0;
    for (i = k + 1; i <= n; i++) {
      sn = sn - tmp_pef[i] * tmp_per[i];
      sd = sd + tmp_pef[i]*tmp_pef[i] + tmp_per[i]*tmp_per[i];
    }
    
    if (sd == 0) {
      return Error::handle(name(), L"computeLatticeBurg", Error::ARG, __FILE__, __LINE__);
    }      
    akk = 2.0 * sn / sd;
    tmp_rc[k] = akk;

    // convert rc[] to a[] using the levinson recursion
    //
    tmp_a[k] = akk;
    for (i = 1; i <= k / 2; i++) {
      ai = tmp_a[i];
      aj = tmp_a[k - i];
      tmp_a[i] = ai + akk * aj;
      tmp_a[k - i] = aj + akk * ai;
    }

    // compute new autocorrelation estimate
    //
    sum = 0;
    for (i = 1; i <= k; i++) {
      sum = sum + tmp_a[i] * tmp_r[k - i];
    }
    tmp_r[k] = sum;

    // update prediction-error power
    //
    pe *= 1 - akk*akk;
    if (pe < 0) {
      break;
    }
  }

  // copy predictor coefficients
  //
  if (!pred_coef_a.assign(pc_order, tmp_a)) {
    return Error::handle(name(), L"computeLatticeBurg", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // compute the filter gain
  //
  err_energy_a = 1.0;
  for (i = 1; i <= order_d; i++) {
    tmp = 1 - tmp_rc[i] * tmp_rc[i];
    err_energy_a *= tmp;
  }
  err_energy_a = sqrt(err_energy_a);

  // clean up old space
  //
  delete [] tmp_r;
  delete [] tmp_rc;
  delete [] tmp_a;
  delete [] tmp_pef;
  delete [] tmp_per;

  // exit gracefully
  //
  if (pe >= (float64)0) {
    return true;
  }
  else {
    return Error::handle(name(), L"computeLatticeBurg", ERR_PREDERR,
			 __FILE__, __LINE__, Error::WARNING);
  }
}

// method: computeReflectionStepDown
//
// arguments:
//  VectorFloat& output: (output) predictor coefficients
//  const VectorFloat& input: (input) reflection coefficients
//
// return: a bool8 value indicating status
//
// this method computes the prediction coefficients from the
// reflection coefficents.
//
// reference Eq 26-30,
//  J. Picone, "Signal Modelling Techniques in Speech Recognition"
//  Proceedings of IEEE, pp. 1210-1215, 1993.
//
bool8 Prediction::computeReflectionStepDown(VectorFloat& output_a,
					      const VectorFloat& input_a) {

  // a backwards index
  //
  int32 j_bckwd;
  
  // inner loop limit
  //
  int32 middle_index;
  
  // accumulator
  //
  float64 sum;
  
  // register
  //
  float64 rc_reg;                             
  
  // initialization
  //
  int32 order = input_a.length();
  output_a.setLength(order + 1);
  output_a(0) = 1;
  output_a(1) = input_a(0);
  
  for (int32 i = 2; i <= order; i++) {
    rc_reg = input_a(i - 1);
    output_a(i) = rc_reg;
    j_bckwd = i - 1;
    middle_index = i / 2;
    
    for (int32 j = 1; j <= middle_index; j++) {
      sum = output_a(j_bckwd) + (float)rc_reg * output_a(j); 
      output_a(j) = output_a(j) + (float)rc_reg * output_a(i - j);
      output_a(j_bckwd) = sum;
      j_bckwd--;
    }
  }

  // exit gracefully
  //
  return true;
}

// method: computeLogAreaKellyLochbaum
//
// arguments:
//  VectorFloat& output: (output) predictor
//  const VectorFloat& input: (input) log erea ratio coefficients
//
// return: a bool8 value indicating status
//
// this method computes the prediction coefficients from the
// log-area-ratio, using LOG_AREA_RATIO algorithm and KELLY_LOCHBAUM
// implementation
//
bool8 Prediction::computeLogAreaKellyLochbaum(VectorFloat& output_a,
						const VectorFloat& input_a) {

  // declare local variables
  //
  bool8 status = false;
  int32 length = input_a.length();

  // compute the reflection coefficients from log-area-ratio input
  //
  Reflection refl(Reflection::LOG_AREA_RATIO,
		  Reflection::KELLY_LOCHBAUM,
		  length,
		  Reflection::DEF_DYN_RANGE);

  VectorFloat refl_coef;
  status = refl.compute(refl_coef, input_a, AlgorithmData::LOG_AREA_RATIO);

  // compute the prediction coefficients from reflection coefficients
  //
  setAlgorithm(Prediction::REFLECTION);
  setImplementation(Prediction::STEP_DOWN);
  status = compute(output_a, refl_coef, AlgorithmData::REFLECTION);

  // set variable value back
  //
  setAlgorithm(Prediction::LOG_AREA_RATIO);
  setImplementation(Prediction::KELLY_LOCHBAUM);
  
  // exit gracefully
  //
  return status;
}