corr_05.cc

这是一个从音频信号里提取特征参量的程序

// file: $isip/class/algo/Correlation/corr_05.cc
// version: $Id: corr_05.cc,v 1.30 2002/05/31 21:57:24 picone Exp $
//

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

// method: apply
//
// arguments:
//  Vector<AlgorithmData>& output: (output) output data
//  const Vector< CircularBuffer<AlgorithmData> >& input: (input) input data
//
// return: a boolean value indicating status
//
// this method selects the appropriate computation methods according
// to the input data type
//
boolean Correlation::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
  //
  long len = input_a.length();
  output_a.setLength(len);
  boolean res = true;

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

  // loop over the channels and call the compute method
  //
  for (long c = 0; c < len; c++) {

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

    // branch on algorithm: autocorrelation only has one input,
    // crosscorrelation has two inputs
    //
    if (input_a(c)(0).getDataType() == AlgorithmData::VECTOR_FLOAT) {

      // 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),
		     input_a(c)(0).getCoefType(),
		     c);
    }

    // for crosscorrelation the input should be a combination of two
    // VectorFloat's -- one for each signal to obtain the correlation
    // between.
    //
    else if (input_a(c)(0).getDataType() == AlgorithmData::COMBINATION) {

      // two input VectorFloat signals are expected here
      //
      if (input_a(c)(0).getCombination().length() != 2) {
	return Error::handle(name(), L"apply", ERR_INPUT, __FILE__, __LINE__);
      }

      // force the output to be a vector through
      // AlgorithmData::makeVectorFloat. for each of the two input signals
      // call getVectorFloat to obtain the VectorFloat signal.
      //
      res &= compute(output_a(c).makeVectorFloat(),
		     input_a(c)(0).getCombination()(0).getVectorFloat(),
		     input_a(c)(0).getCombination()(1).getVectorFloat(),
		     input_a(c)(0).getCoefType(),
		     c);
    }

    // error condition
    //
    else {
      return Error::handle(name(), L"apply", ERR_UNKTYP, __FILE__, __LINE__);
    }

    // set the coefficient type for the output
    //
    // branch on algorithm:
    //  Algorithm: AUTO or AUTO_TAPERED
    //
    if ((algorithm_d == AUTO) || (algorithm_d == AUTO_TAPERED)) {
      output_a(c).setCoefType(AlgorithmData::CORRELATION);
    }

    //  Algorithm: CROSS or CONV
    //
    else {
      output_a(c).setCoefType(AlgorithmData::SIGNAL);
    }
  }

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

  // exit gracefully
  //
  return res;
}

// method: compute
//
// arguments:
//  VectorFloat& output: (output) output data
//  const CircularBuffer<AlgorithmData>& input: (input) input data
//  AlgorithmData::COEF_TYPE coef_type: (input) coeff type of input signal
//  long channel_index: (input) the channel index of input signal
//
// return: a boolean value containing status
//
// this method applies the specified window on the input channel
//
boolean Correlation::compute(VectorFloat& output_a,
			const CircularBuffer<AlgorithmData>& input_a,
			AlgorithmData::COEF_TYPE coef_type_a,
			long channel_index_a) {

  // for FRAME_INTERNAL mode, we just process what data we have
  //
  if (cmode_d == FRAME_INTERNAL) {
    return compute(output_a, input_a(0).getVectorFloat(), coef_type_a,
		   channel_index_a);
  }
  
  // for other modes, we error
  //
  else if (cmode_d != CROSS_FRAME) {
    return Error::handle(name(), L"compute",
			 ERR_UNSUPM, __FILE__, __LINE__);
  }

  // make sure there is enough data
  //
  if ((input_a.getNumForward() < getLeadingPad()) ||
      (input_a.getNumBackward() < getTrailingPad())) {
    return Error::handle(name(), L"compute", ERR_DATA, __FILE__, __LINE__);
  }

  // precompute some useful constants
  //
  double frame_nsamp = frame_dur_d * sample_freq_d;
  long order = order_d;

  double num_frames = order / frame_nsamp;
  long frame = (long)Integral::ceil(num_frames);
  
  VectorFloat buf;

  // algorithm: AUTO
  //
  if (algorithm_d == AUTO) {
    
    // this is the number of samples not included from the first frame
    //
    long index = frame ;
    long first_offset = frame * (long)Integral::ceil(frame_nsamp) - order;
    long first_nelem = (long)Integral::ceil(frame_nsamp) - first_offset;
    double num_left = order + frame_nsamp;
    
    // loop over the remaining data
    //
    while (num_left > 0) {
      
      // for the first frame, grab only the non-zero samples
      //
      if (index == frame) {
	buf.move(input_a(-index).getVectorFloat(),
		 first_nelem, first_offset, 0);
	num_left -= first_nelem;
      }
      
      // for all other frames, we copy the entire frame
      //
      else {
	buf.concat(input_a(-index).getVectorFloat());
	num_left -= frame_nsamp;
      }

      // decrement the frame counter
      //
      index--;
    }
    
    // apply the correlation function
    //
    if (implementation_d == FACTORED) {
      return computeAutoFactoredFromSignal(output_a, buf);
    }
    else if (implementation_d == UNFACTORED) {
      return computeAutoUnfactoredFromSignal(output_a, buf);  
    }
    else
      return Error::handle(name(), L"compute",
			   ERR_UNSUPI, __FILE__, __LINE__);
  }

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

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

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

// method: compute
//
// arguments:
//  VectorFloat& output: (output) output data
//  const VectorFloat& input: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  long channel_index: (input) channel index
//
// return: a boolean value indicating status
//
// this method selects the appropriate computational method to compute
// autocorrelation.
//
boolean Correlation::compute(VectorFloat& output_a,
			     const VectorFloat& input_a,
			     AlgorithmData::COEF_TYPE input_coef_type_a,
			     long channel_index_a) {

  // declare local variable
  //
  boolean status = false;

  // algorithm: AUTO or AUTO_TAPERED
  //
  if ((algorithm_d == AUTO) || (algorithm_d == AUTO_TAPERED)) {

    // input coefficient type: linear prediction coefficients
    //
    if (input_coef_type_a == AlgorithmData::PREDICTION) {
      status = computeAutoFromPrediction(output_a, input_a);
    }
    else if (input_coef_type_a == AlgorithmData::REFLECTION) {
      status = computeAutoFromReflection(output_a, input_a);
    }

    // else: signals
    //
    else if (input_coef_type_a == AlgorithmData::SIGNAL) {

      // check the implementation
      //
      if (implementation_d == FACTORED) {
	status = computeAutoFactoredFromSignal(output_a, input_a);
      }
      else if (implementation_d == UNFACTORED) {
	status = computeAutoUnfactoredFromSignal(output_a, input_a);
      }
      else {
	return Error::handle(name(), L"compute",
			     ERR_UNKIMP, __FILE__, __LINE__);
      }
    }
      
    // error: unknown coefficient type
    //
    else {
      return Error::handle(name(), L"compute",
			   ERR_UNCTYP, __FILE__, __LINE__);
    }
  }

  // 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:
//  VectorFloat& output: (output) output data
//  const VectorFloat& input1: (input) input data
//  const VectorFloat& input2: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  long channel_index: (input) channel index
//
// return: a boolean value indicating status
//
// this method selects the appropriate computational method to compute
// cross-correlation and convolution
//
boolean Correlation::compute(VectorFloat& output_a,
			     const VectorFloat& input1_a,
			     const VectorFloat& input2_a,
			     AlgorithmData::COEF_TYPE input_coef_type_a,
			     long channel_index_a) {
  
  // declare local variable
  //
  boolean status = false;
  
  // algorithm: CROSS
  //
  if (algorithm_d == CROSS) {

    // implementation: UNFACTORED
    //
    if (implementation_d == UNFACTORED) {
      status = computeCross(output_a, input1_a, input2_a);
    }

    // check known but unsupported implementations
    //
    else if (implementation_d == FACTORED) {
      return Error::handle(name(), L"compute",
			   ERR_UNSUPA, __FILE__, __LINE__);
    }

    // error: unknown implementation
    //
    else {
      return Error::handle(name(), L"compute",
			   ERR_UNKIMP, __FILE__, __LINE__);
    }
  }
  
  // algorithm: CONVOLUTION
  //
  else if (algorithm_d == CONV) {

    // implementation: UNFACTORED
    //
    if (implementation_d == UNFACTORED) {
      status = computeConv(output_a, input1_a, input2_a);
    }

    // check known but unsupported implementations
    //
    else if (implementation_d == CIRCULAR) {
      return Error::handle(name(), L"compute",
			   ERR_UNSUPA, __FILE__, __LINE__);
    }

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

  // possibly display the data. note that we generate an input vector
  // if the debug level warrants that it will actually be output by
  // the AlgorithmBase method -- if not we just pass in a dummy value.
  //
  if (debug_level_d >= Integral::DETAILED) {
    Vector<VectorFloat> input(2);
    input(0).assign(input1_a);
    input(1).assign(input2_a);
    display(output_a, input, name());
  }
  else {
    display(output_a, input1_a, name());
  }
  
  // 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
//  long channel_index: (input) channel index
//
// return: a boolean value indicating status
//
// this method selects the appropriate computational method for
// the complex data type.
//
boolean Correlation::compute(VectorComplexFloat& output_a,
			     const VectorComplexFloat& input_a,
			     AlgorithmData::COEF_TYPE input_coef_type_a,
			     long channel_index_a) {
  
  // this method is not implemented yet
  //
  return Error::handle(name(), L"compute", Error::NOT_IMPLEM,
		       __FILE__, __LINE__);      
}

// method: compute
//
// arguments:
//  VectorComplexFloat& output: (output) output data
//  const VectorComplexFloat& input1: (input) input data
//  const VectorComplexFloat& input2: (input) input data
//  AlgorithmData::COEF_TYPE input_coef_type: (input) type of input
//  long channel_index: (input) channel index
//
// return: a boolean value indicating status
//
// this method selects the appropriate computational method for
// complex data type
//
boolean Correlation::compute(VectorComplexFloat& output_a,
			     const VectorComplexFloat& input1_a,
			     const VectorComplexFloat& input2_a,
			     AlgorithmData::COEF_TYPE input_coef_type_a,
			     long channel_index_a) {

  // this method is not implemented yet
  //
  return Error::handle(name(), L"compute", Error::NOT_IMPLEM,
		       __FILE__, __LINE__);      
}

// method: computeAutoFactoredFromSignal
//
// arguments:
//  VectorFloat& output: (output) Correlation coefficients
//  const VectorFloat& input: (input) input data
//
// return: a boolean value indicating status
//
// this method computes the Autocorrelation function using a factored
// Correlation computation (for sampled data). see:
//
//  J.D. Markel and A.H. Gray, Jr.,
//  Linear Prediction of Speech, Springer-Verlag Berlin Heidelberg,
//  New York, New York, USA, pp. 218, 1976.
//
// see figure 9.2 for a description of the algorithm. note that for
// this algorithm to work, the order must be less than the length/3
//
boolean Correlation::computeAutoFactoredFromSignal(VectorFloat& output_a,
						   const VectorFloat& input_a) {

  // declare local variable
  //
  long order = order_d;
  long N = input_a.length();
  boolean status = false;
  
  // resize the output
  //
  if (!output_a.setLength(order + 1)) {
    return Error::handle(name(), L"computeAutoFactoredFromSignal",
			 ERR, __FILE__, __LINE__);
  }
  
  // algorithm: AUTO
  //  in this approach, we consider the sum from order to N-1
  //
  if (algorithm_d == AUTO) {

    // compute the zeroth term explicitly over the range [order,N]
    //
    double sum = 0.0;
    for (long i = order; i < N; i++) {
      sum += input_a(i) * input_a(i);
    }
    output_a(0) = sum;

    // compute lags on the range [1,order]:
    //  loop over the order
    //
    for (long i = 1; i <= order; i++) {

      // reset the accumulators
      //
      long next = order;
      long last = next;
      sum = 0.0;

      // make sure we don't overrun the end of the vector
      //
      while (last < N) {

	// loop over all partial sums in this block.
	// make sure we don't overrun the edge of the vector
	//
	long np  = next;
	long npm = next - i;
	long npp = next + i;

	for (long j = 0; j < i; j++) {
	  if (npp < N) {
	    sum += input_a(np++) * (input_a(npm++) + input_a(npp++));
	  }
	}

	// increase counters
	//
	next += 2*i;
	last = next + 2*i;
      }
      
      // accumulate the leftover stuff (products that can no longer be grouped
      // because they are too close to the edge of the frame
      //
      long nmi = next - i;
      for (long n = next; n < N; n++) {
	sum += input_a(n) * input_a(nmi++);
      }

      // output the value
      //
      output_a(i) = sum;
    }

    // everything is ok
    //
    status = true;
  }

  // algorithm: AUTO_TAPERED
  //  note that this implementation follows the Markel and Gray reference
  //  above (pg. 217) very closely.
  //
  else {

    // check the order
    //
    if ((order_d < (long)1) ||
	(order_d > (input_a.length() / 3))) {
      return Error::handle(name(), L"computeAutoFactoredFromSignal",
			   ERR_ORDER, __FILE__, __LINE__);
    }  
    
    // declare some local variables
    //
    long mp = (long)order_d + 1;
    long n = input_a.length();
    double rr;
    double tmp;
    
    long i, j, k, im1, mod, nterm, kstrt, kstop, kk;
    long kl = (long)0;
    long kr = (long)0;
    long j_start;
    
    // initialize
    //
    rr = 0.0;
    for (i = 0; i < n; i++) {
      tmp = (double)input_a(i);
      rr += tmp*tmp;
    }
    output_a((long)0) = rr;
    
    // loop over all lags
    //
    for (i = 1; i < mp; i++) {
      
      // compute a new value for rr
      //
      im1 = i;
      mod = im1 + im1;
      nterm = n - im1;