// file: $isip/class/algo/Calculus/calc_05.cc
// version: $Id: calc_05.cc 8226 2002-06-25 01:42:42Z picone $
//

// isip include files
//
#include "Calculus.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 method
//
bool8 Calculus::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;
  frame_index_d++;
  
  // 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);
    
    if (input_a(c)(0).getDataType() != AlgorithmData::VECTOR_FLOAT) {
      return Error::handle(name(), L"apply", ERR_UNSUPA, __FILE__, __LINE__);
    }

    // cmode_d: FRAME_INTERNAL
    //
    if (cmode_d == FRAME_INTERNAL) {
      
      // call AlgorithmData::makeVectorFloat to force the output vector for
      // this channel to be a VectorFloat, use the
      // CircularBuffer<AlgorithmData> directly on the input for this
      // channel
      //          
      res &= computeInternal(output_a(c).makeVectorFloat(), input_a(c));
    }
    
    // cmode_d: CROSS_FRAME
    //
    else if (cmode_d == CROSS_FRAME) {
      
      // call AlgorithmData::makeVector to force the output vector for
      // this channel to be a VectorFloat, use the
      // CircularBuffer<AlgorithmData> directly on the input for this
      // channel
      //          
      res &= computeCross(output_a(c).makeVectorFloat(), input_a(c));
    }
    
    // mode_d: ACCUMULATE
    //
    else {
      return Error::handle(name(), L"apply",
			   ERR_UNSUPM, __FILE__, __LINE__);
    }    
    // set the coefficient type for the output
    //    
    output_a(c).setCoefType(input_a(c)(0).getCoefType());
  }
  
  // finish the debugging output
  //
  displayFinish(this);
  
  // exit gracefully
  //
  return res;
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) output data
//  const VectorFloat& input: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel index
//
// return: a bool8 value indicating status
//
// this method performs appropriate calculus operation on input vector.
// note that none of the algorithms or implementations currently
// require memory to be maintained between calls, so channel_index is
// not used.
//
  bool8 Calculus::compute(VectorFloat& output_a, const VectorFloat& input_a,
			  AlgorithmData::COEF_TYPE input_coef_type_a,
			  int32 channel_index_a) {

  // declare local variables
  //
  bool8 status = false;
  
  // check valid flag
  //
  if (!is_valid_d) {
    init();
  }

  // branch on algorithm:
  //  Algorithm: DIFFERENTIATION
  //
  if (algorithm_d == DIFFERENTIATION) {

    // branch on implementation:
    //
    //  Implementation: REGRESSION
    //
    if (implementation_d == REGRESSION) {
      status = computeDifRegression(output_a, input_a);
    }

    //  Implementation: CENTRAL_DEFFERENCE
    //
    else if (implementation_d == CENTRAL_DIFFERENCE) {
      status = computeDifCentral(output_a, input_a);
    }

    //  Implementation: BACKWARD_DIFFERENCE
    //
    else if (implementation_d == BACKWARD_DIFFERENCE) {
      status = computeDifBackward(output_a, input_a);
    }

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

  //  Algorithm: INTEGRATION
  //
  else if (algorithm_d == INTEGRATION) {
    return Error::handle(name(), L"compute",
			 ERR_UNSUPA, __FILE__, __LINE__);
  }  

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

  // exit gracefully
  //
  return status;
}

// method: compute
//
// arguments:
//  VectorComplexFloat& output: (output) output data
//  const VectorComplexFloat& input: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 index: (input) channel index
//
// return: a bool8 value indicating status
//
// this method performs the appropriate calculus operation on the input vector
//
bool8 Calculus::compute(VectorComplexFloat& output_a,
			  const VectorComplexFloat& input_a,
			  AlgorithmData::COEF_TYPE input_coef_type_a,
			  int32 index_a) {
  return Error::handle(name(), L"compute", ERR_UNSUPA, __FILE__, __LINE__);
}  

// method: computeInternal
//
// arguments:
//  VectorFloat& output: (output) output vector
//  const CircularBuffer<AlgorithmData>& input: (input) input vectors
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel index
//
// return: a bool8 value indicating status
//
// this method transfers data from the circular buffer to the
// lower level compute method
//
bool8 Calculus::computeInternal(VectorFloat& output_a,
				  const CircularBuffer<AlgorithmData>& input_a,
				  AlgorithmData::COEF_TYPE input_coef_type_a,
				  int32 channel_index_a) {
  
  // check that the type of data is correct
  //
  if (input_a(0).getDataType() != AlgorithmData::VECTOR_FLOAT) {
    return Error::handle(name(), L"computeInternal", ERR_UNKTYP, __FILE__, __LINE__);
  }
    
  // create space for the output
  //
  int32 num_coeffs = input_a(0).getVectorFloat().length();
  output_a.setLength(num_coeffs);
    
  // check how much past and future data is in the buffer, and whether
  // this is adequate. note that the data must span from the interval:
  //
  //  [offset_d, offset_d+num_elements-1]
  //
  // getNumBackward and getNumForward retrieve this information from
  // the circular buffer. the calling class, typically Algorithm,
  // has the responsibility of guaranteeing there is enough data in
  // the circular buffer
  //
  if (offset_d < (-1 * input_a.getNumBackward())) {
    return Error::handle(name(), L"computeInternal", ERR_INSUFD,
			 __FILE__, __LINE__);
  }
  else if ((offset_d + num_elements_d - 1) > input_a.getNumForward()) {
    return Error::handle(name(), L"computeInternal", ERR_INSUFD,
			 __FILE__, __LINE__);
  }
    
  // declare some scratch buffers:
  //  a vector to hold frame-spanning data
  //  a vector to accept output from compute()
  //
  VectorFloat data(num_elements_d);
  VectorFloat result(1);

  // loop over the number of coefficients
  //
  for (int32 i = 0; i < num_coeffs; i++) {
      
    // copy data onto the data buffer
    //
    for (int32 k = 0; k < num_elements_d; k++) {
	
      // fetch the data backward in time
      //
      if (i +  k < delta_win_d) {
	data(k) = (input_a(-1).getVectorFloat())(num_coeffs -
						 delta_win_d + i + k);
      }
      else if ( (i > num_coeffs - 1 - delta_win_d) &&
		(i + k - delta_win_d > num_coeffs - 1) &&
		(k > delta_win_d)) {
	data(k) = (input_a(1).getVectorFloat())(i + k -
						num_coeffs - delta_win_d);
      }
      else {
	data(k) = (input_a(0).getVectorFloat())(i + k - delta_win_d);
      }
    }

    // compute the algorithm: we now have the temporal data in a vector
    //  and can process it
    //
    if (!compute(result, data, input_coef_type_a)) {
      return Error::handle(name(), L"computeInternal",
			   ERR_COMPF, __FILE__, __LINE__, Error::WARNING);
    }

    // save the result into the output buffer and move to the next element
    //
    output_a(i) = result(0);
  }

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

  // exit gracefully
  //
  return true;
}

// method: computeCross
//
// arguments:
//  VectorFloat& output: (output) output vector
//  const CircularBuffer<AlgorithmData>& input: (input) input vectors
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  int32 channel_index: (input) channel index
//
// return: a bool8 value indicating status
//
// this method transfers data from the circular buffer to the
// lower level compute method
//
bool8 Calculus::computeCross(VectorFloat& output_a,
			       const CircularBuffer<AlgorithmData>& input_a,
			       AlgorithmData::COEF_TYPE input_coef_type_a,
			       int32 channel_index_a) {
  
  // check that the type of data is correct
  //
  if (input_a(0).getDataType() != AlgorithmData::VECTOR_FLOAT) {
    return Error::handle(name(), L"computeCross",
			 ERR_UNKTYP, __FILE__, __LINE__);
  }
    
  // create space for the output
  //
  int32 num_coeffs = input_a(0).getVectorFloat().length();
  output_a.setLength(num_coeffs);
    
  // check how much past and future data is in the buffer, and whether
  // this is adequate. note that the data must span from the interval:
  //
  //  [offset_d, offset_d+num_elements-1]
  //
  // getNumBackward and getNumForward retrieve this information from
  // the circular buffer. the calling class, typically Algorithm,
  // has the responsibility of guaranteeing there is enough data in
  // the circular buffer
  //
  if (offset_d < (-1 * input_a.getNumBackward())) {
    return Error::handle(name(), L"computeCross",
			 ERR_INSUFD, __FILE__, __LINE__);
  }
  else if ((offset_d + num_elements_d - 1) > input_a.getNumForward()) {
    return Error::handle(name(), L"computeCross",
			 ERR_INSUFD, __FILE__, __LINE__);
  }
    
  // declare some scratch buffers:
  //  a vector to hold frame-spanning data
  //  a vector to accept output from compute()
  //
  VectorFloat data(num_elements_d);
  VectorFloat result(1);    
  int32 number_frame = (int32)(Integral::round(signal_duration_d / frame_dur_d));

  // loop over the number of coefficients
  //
  for (int32 i = 0; i < num_coeffs; i++) {

    int32 j = offset_d;

    // copy data onto the data buffer
    //
    for (int32 k = 0; k < num_elements_d; k++) {

      // fetch the data backward in time
      //
      if ((k + frame_index_d + offset_d) < 0) {
	data(k) = (input_a(0 - frame_index_d).getVectorFloat())(i);
	j++;
      }
      else if ((k + frame_index_d + offset_d) > number_frame - 1) {
	data(k) =
	  (input_a(number_frame - 1 - frame_index_d ).getVectorFloat())(i);
      }
      else {	
	data(k) = (input_a(j++).getVectorFloat())(i);
      }
    }
    
    // compute the algorithm: we now have the temporal data in a vector
    //  and can process it
    //
    if (!compute(result, data, input_coef_type_a)) {
      return Error::handle(name(), L"computeCross",
			   ERR_COMPF, __FILE__, __LINE__, Error::WARNING);
    }

    // save the result into the output buffer and move to the next element
    //
    output_a(i) = result(0);
  }

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

  // exit gracefully
  //
  return true;
}


// method: computeDifCentral
//
// arguments:
//  VectorFloat& output: (output) output data
//  const VectorFloat& input: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes a temporal derivative using a central difference.
// two pertinent references are:
//
//  J.G. Proakis and D.G. Manolakis, Digital Signal Processing (Third Edition)
//  Prentice-Hall, Upper Saddle River, New Jersey, USA, 1996.
//
// and:
//
//  J. Picone, "Signal Processing in Speech Recognition," Proceedings
//  of the Speech Recognition System Training Workshop, Institute
//  for Signal and Information Processing, Mississippi State University,
//  May 2000 (see http://www.isip.msstate.edu/conferences/srstw00/program/
//  session_04/signal_processing/index.html).
//
bool8 Calculus::computeDifCentral(VectorFloat& output_a,
				    const VectorFloat& input_a) {

  // check the order: we currently only support a first-order difference
  //
  if (order_d != (int32)1) {
    return Error::handle(name(), L"computeDifCentral",
			 ERR_ORDER, __FILE__, __LINE__);
  }

  // check the length
  //
  else if (input_a.length() != (2 * delta_win_d + 1)) {
    return Error::handle(name(), L"computeDifCentral",
			 ERR_INSUFD, __FILE__, __LINE__);
  }
  
  // compute the numerator of the difference
  //
  int32 N = (int32)delta_win_d;
  float64 num = input_a(N + delta_win_d) - input_a(N - delta_win_d);

  // exit gracefully: return the value
  //
  output_a.setLength(1);

  if (denom_d != (float64)0.0) {
    output_a(0) = num * denom_d;
    return true;
  }    
  else {
    return Error::handle(name(), L"computeDifCentral",
			 ERR_ZERODENOM, __FILE__, __LINE__, Error::WARNING);
  }
}

// method: computeDifBackward
//
// arguments:
//  VectorFloat& output: (output) output data
//  const VectorFloat& input: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes a temporal derivative using a simple
// backward difference. two pertinent references are:
//
//  J.G. Proakis and D.G. Manolakis, Digital Signal Processing (Third Edition)
//  Prentice-Hall, Upper Saddle River, New Jersey, USA, 1996.
//
// and:
//
//  J. Picone, "Signal Processing in Speech Recognition," Proceedings
//  of the Speech Recognition System Training Workshop, Institute
//  for Signal and Information Processing, Mississippi State University,
//  May 2000 (see http://www.isip.msstate.edu/conferences/srstw00/program/
//  session_04/signal_processing/index.html).
//
bool8 Calculus::computeDifBackward(VectorFloat& output_a,
				     const VectorFloat& input_a) {

  // check the order: we currently only support a first-order difference
  //
  if (order_d != (int32)1) {
    return Error::handle(name(), L"computeDifBackward",
			 ERR_ORDER, __FILE__, __LINE__);
  }

  // check the length
  //
  else if (input_a.length() != (delta_win_d + (int32)1)) {
    return Error::handle(name(), L"computeDifBackward",
			 ERR_INSUFD, __FILE__, __LINE__);
  }
  
  // compute the numerator of the difference
  //
  int32 N = (int32)delta_win_d;
  float64 num = input_a(N) - input_a(0);

  // exit gracefully: return the value
  //
  output_a.setLength(1);

  if (denom_d != (float64)0.0) {
    output_a(0) = num * denom_d;
    return true;
  }    
  else {
    return Error::handle(name(), L"computeDifBackward",
			 ERR_ZERODENOM, __FILE__, __LINE__, Error::WARNING);
  }
}

// method: computeDifRegression
//
// arguments:
//  VectorFloat& output: (output) output data
//  const VectorFloat& input: (input) input data
//
// return: a bool8 value indicating status
//
// this method computes a temporal derivative using a regression
// approximation. Two pertinent references are:
//
//  F.K. Soong and A.E. Rosenburg, "On the Use of Instantaneous and
//  Transitional Spectral Information in Speaker Recognition," Proceedings
//  of the International Conference on Acoustics, Speech, and Signal
//  Processing, Tokyo, Japan, pp. 877-880, April 1986.
//
// and:
//
//  J. Picone, "Signal Processing in Speech Recognition," Proceedings
//  of the Speech Recognition System Training Workshop, Institute
//  for Signal and Information Processing, Mississippi State University,
//  May 2000 (see http://www.isip.msstate.edu/conferences/srstw00/program/
//  session_04/signal_processing/index.html).
//
bool8 Calculus::computeDifRegression(VectorFloat& output_a,
				       const VectorFloat& input_a) {
  
  // check the order: we currently only support a first-order regression
  //
  if (order_d != (int32)1) {
    return Error::handle(name(), L"computeDifRegression",
			 ERR_ORDER, __FILE__, __LINE__);
  }

  // check the length
  //
  else if (input_a.length() != (2 * delta_win_d + 1)) {
    return Error::handle(name(), L"computeDifRegression",
			 ERR_INSUFD, __FILE__, __LINE__);
  }

  // compute the numerator of the regression formula 
  //
  float64 num = 0.0;
  int32 N = (int32)delta_win_d;
  
  for (int32 i = 1; i <= N; i++) {
    num += (float64)i * (input_a(N + i) - input_a(N - i));
  }

  // exit gracefully: return the value
  //
  output_a.setLength(1);

  if (denom_d != (float64)0.0) {
    output_a(0) = num * denom_d;
    return true;
  }    
  else {
    return Error::handle(name(), L"computeDifRegression",
			 ERR_ZERODENOM, __FILE__, __LINE__, Error::WARNING);
  }
}