// file: $isip/class/algo/Spectrum/spec_05.cc
// version: $Id: spec_05.cc 8275 2002-07-02 23:20:05Z picone $
//

// isip include files
//
#include "Spectrum.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 calls the appropriate computation methods
//
bool8 Spectrum::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);

  // start the debugging output
  //
  displayStart(this);

  // accumulate a result status -- if compute method returns false
  // this apply method should return false
  //  
  bool8 res = true;
  
  // loop over the channels and call the compute method
  //
  for (int32 c = 0; c < len; c++) {

    // display the channel number
    //
    displayChannel(c);

    if (implementation_d == MAGNITUDE) {

      // 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);
    }

    else if (implementation_d == COMPLEX) {

      // call AlgorithmData::makeComplexVectorFloat 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).makeVectorComplexFloat(),
		     input_a(c)(0).getVectorFloat(),
		     input_a(c)(0).getCoefType(), c);
    }

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

    // set the coefficient type of output
    //
    output_a(c).setCoefType(AlgorithmData::SPECTRUM);
  }
  
  // finish the debugging output
  //
  displayFinish(this);

  // exit gracefully
  //
  return res;  
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 index: (input) channel number
//
// return: a bool8 value indicating status
//
// this method computes the magnitude spectrum of the input signal
//
bool8 Spectrum::compute(VectorFloat& output_a,
			  const VectorFloat& input_a,
			  AlgorithmData::COEF_TYPE input_coef_type_a,
			  int32 index_a) {

  // declare local variables
  //
  bool8 status = false;
  
  // brach on algorithm
  //
  // algorithm: FOURIER
  //
  if (algorithm_d == FOURIER) {

    // coefficient type: SIGNAL
    //
    if (input_coef_type_a == AlgorithmData::SIGNAL) {

      // calculate corresponding fourier transform compute method
      //
      status = computeFourierMag(output_a, input_a);
    }

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

  // check known but unsupported algorithms
  //
  else if (algorithm_d == MAXIMUM_ENTROPY) {

    // coefficient type: REFLECTION
    //
    if (input_coef_type_a == AlgorithmData::REFLECTION) {

      // calculate corresponding fourier transform compute method
      //
      status = computeReflectionMag(output_a, input_a, dyn_range_d);
    }

    // coefficient type: PREDICTION
    //
    else if (input_coef_type_a == AlgorithmData::PREDICTION) {

      // calculate corresponding fourier transform compute method
      //
      status = computePredictionMag(output_a, input_a);
    }

    // coefficient type: CORRELATION
    //
    else if (input_coef_type_a == AlgorithmData::CORRELATION) {

      // calculate corresponding fourier transform compute method
      //
      status = computeCorrelationMag(output_a, input_a, dyn_range_d);
    }        

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

  // else: error unknown algorithm
  //
  else {
    status = 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:
//  VectorComplexFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 index: (input) channel number
//
// return: a bool8 value indicating status
//
// this method computes the complex spectrum of the input signal
//
bool8 Spectrum::compute(VectorComplexFloat& output_a,
			  const VectorFloat& input_a,
			  AlgorithmData::COEF_TYPE input_coef_type_a,
			  int32 index_a) {

  // declare local variables
  //
  bool8 status = false;
  
  // brach on algorithm
  //
  // algorithm: FOURIER
  //
  if (algorithm_d == FOURIER) {

    // coefficient type: SIGNAL
    //
    if (input_coef_type_a == AlgorithmData::SIGNAL) {

      // calculate corresponding fourier transform compute method
      //
      status = computeFourierComplex(output_a, input_a);
    }

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

  // check known but unsupported algorithms
  //
  else if (algorithm_d == MAXIMUM_ENTROPY) {

    // coefficient type: REFLECTION
    //
    if (input_coef_type_a == AlgorithmData::REFLECTION) {

      // calculate corresponding fourier transform compute method
      //
      status = computeReflectionComplex(output_a, input_a);
    }

    // coefficient type: PREDICTION
    //
    else if (input_coef_type_a == AlgorithmData::PREDICTION) {

      // calculate corresponding fourier transform compute method
      //
      status = computePredictionComplex(output_a, input_a);
    }

    // coefficient type: CORRELATION
    //
    else if (input_coef_type_a == AlgorithmData::CORRELATION) {

      // calculate corresponding fourier transform compute method
      //
      status = computeCorrelationComplex(output_a, input_a);
    }        

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

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

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

  // exit gracefully
  //
  return status;  
}

// method: computeReflectionMag
//
// arguments:
//  VectorFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//  float64 range: (input) dynamic range
//
// return: a bool8 value indicating status
//
// this method computes the magnitude spectrum of the input signal for
// REFLECTION implementation type
//
bool8 Spectrum::computeReflectionMag(VectorFloat& output_a,
				       const VectorFloat& input_a,
				       float64 range_a) {


  // declare local variables
  //
  Prediction predc;
  VectorFloat pred_coeff;
  VectorFloat refl_coeff(input_a);

  // set the order of the prediction coefficients
  //
  int32 order = refl_coeff.length() + 1;
  pred_coeff.setLength(order);
  
  // determine the prediction coefficients form the reflection coefficients
  //
  predc.set(Prediction::REFLECTION, Prediction::STEP_DOWN, order, range_a);
  predc.compute(pred_coeff, refl_coeff, AlgorithmData::REFLECTION);

  // negate the predictor coefficients
  //
  pred_coeff.neg();

  // compute the magnitude spectrum
  //
  computeFourierMag(output_a, pred_coeff);
  
  // exit gracefully
  //
  return true;
}

// method: computePredictionMag
//
// arguments:
//  VectorFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//
// return: a bool8 value indicating status
//
// this method computes the magnitude spectrum of the input signal for
// PREDICTION implementation type
//
bool8 Spectrum::computePredictionMag(VectorFloat& output_a,
				       const VectorFloat& input_a) {

  // declare local variables
  //
  VectorFloat pred_coeff(input_a);

  // negate the predictor coefficients
  //
  pred_coeff.neg();

  // compute the magnitude spectrum
  //
  computeFourierMag(output_a, pred_coeff);
  
  // exit gracefully
  //
  return true;
}

// method: computeCorrelationMag
//
// arguments:
//  VectorFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//  float64 range: (input) dynamic range
//
// return: a bool8 value indicating status
//
// this method computes the magnitude spectrum of the input signal for
// CORRELATION implementation type
//
bool8 Spectrum::computeCorrelationMag(VectorFloat& output_a,
					const VectorFloat& input_a,
					float64 range_a) {

  // declare local variables
  //
  Prediction predc;
  VectorFloat pred_coeff;
  VectorFloat auto_coeff(input_a);

  // set the order of the prediction coefficients
  //
  int32 order = auto_coeff.length();
  pred_coeff.setLength(order);

  // compute the prediction coefficients from the autocorrelation coefficients
  //
  predc.set(Prediction::AUTOCORRELATION, Prediction::DURBIN, order, range_a);
  predc.compute(pred_coeff, auto_coeff, AlgorithmData::CORRELATION);

  // negate the predictor coefficients
  //
  pred_coeff.neg();

  // compute the magnitude spectrum
  //
  computeFourierMag(output_a, pred_coeff);
  
  // exit gracefully
  //
  return true;
}

// method: computeFourierMag
//
// arguments:
//  VectorFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//
// return: a bool8 value indicating status
//
// this method computes the magnitude spectrum of the input signal for
// FOURIER implementation type
//
bool8 Spectrum::computeFourierMag(VectorFloat& output_a,
				    const VectorFloat& input_a) {

  // compute the complex spectrum
  //
  VectorComplexFloat tmp;
  computeFourierComplex(tmp, input_a);
  output_a.setLength(tmp.length());

  // compute the magnitude spectrum:
  //  magnitude spectrum = sqrt(real^2 + imag^2)
  //
  for (int32 i = 0; i < tmp.length(); i++) {
    Float sq_spec;
    Float real(tmp(i).real());
    Float imag(tmp(i).imag());
    sq_spec = real.square() + imag.square();
    output_a(i) = sq_spec.sqrt();
  }
  
  // exit gracefully
  //
  return true;
}

// method: computeReflectionComplex
//
// arguments:
//  VectorComplexFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//  float64 range: (input) dynamic range
//
// return: a bool8 value indicating status
//
// this method computes the complex spectrum of the input signal for
// REFLECTION implementation type
//
bool8 Spectrum::computeReflectionComplex(VectorComplexFloat& output_a,
					   const VectorFloat& input_a,
					   float64 range_a) {

  // declare local variables
  //
  Prediction predc;
  VectorFloat pred_coeff;
  VectorFloat refl_coeff(input_a);

  // set the order of the prediction coefficients
  //
  int32 order = refl_coeff.length() + 1;
  pred_coeff.setLength(order);
  
  // determine the prediction coefficients form the reflection coefficients
  //
  predc.set(Prediction::REFLECTION, Prediction::STEP_DOWN, order, range_a);
  predc.compute(pred_coeff, refl_coeff, AlgorithmData::REFLECTION);

  // negate the predictor coefficients
  //
  pred_coeff.neg();
  
  // compute the magnitude spectrum
  //
  computeFourierComplex(output_a, pred_coeff);

  // exit gracefully
  //
  return true;
}

// method: computePredictionComplex
//
// arguments:
//  VectorComplexFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//
// return: a bool8 value indicating status
//
// this method computes the complex spectrum of the input signal for
// PREDICTION implementation type
//
bool8 Spectrum::computePredictionComplex(VectorComplexFloat& output_a,
					   const VectorFloat& input_a) {

  // declare local variables
  //
  VectorFloat pred_coeff(input_a);

  // negate the predictor coefficients
  //
  pred_coeff.neg();
  
  // compute the magnitude spectrum
  //
  computeFourierComplex(output_a, pred_coeff);
  
  // exit gracefully
  //
  return true;
}

// method: computeCorrelationComplex
//
// arguments:
//  VectorComplexFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//  float64 range: (input) dynamic range
//
// return: a bool8 value indicating status
//
// this method computes the complex spectrum of the input signal for
// CORRELATION implementation type
//
bool8 Spectrum::computeCorrelationComplex(VectorComplexFloat& output_a,
					    const VectorFloat& input_a,
					    float64 range_a) {

  // declare local variables
  //
  Prediction predc;
  VectorFloat pred_coeff;
  VectorFloat auto_coeff(input_a);

  // set the order of the prediction coefficients
  //
  int32 order = auto_coeff.length();
  pred_coeff.setLength(order);

  // compute the prediction coefficients from the autocorrelation coefficients
  //
  predc.set(Prediction::AUTOCORRELATION, Prediction::DURBIN, order, range_a);
  predc.compute(pred_coeff, auto_coeff, AlgorithmData::CORRELATION);

  // negate the predictor coefficients
  //
  pred_coeff.neg();
  
  // compute the magnitude spectrum
  //
  computeFourierComplex(output_a, pred_coeff);
  
  // exit gracefully
  //
  return true;
}

// method: computeFourierComplex
//
// arguments:
//  VectorComplexFloat& output: (output) spectrum
//  const VectorFloat& input: (input) signal
//
// return: a bool8 value indicating status
//
// this method computes the spectrum of the input signal for FOURIER
// algorithm type
//
bool8 Spectrum::computeFourierComplex(VectorComplexFloat& output_a,
					const VectorFloat& input_a) {

  // call FourierTransform compute method with a user-specified order
  //
  return ft_d.compute(output_a, input_a);
}