ft_13.cc

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

// file: $isip/class/algo/FourierTransform/ft_13.cc
// version: $Id: ft_13.cc,v 1.5 2001/11/29 21:20:38 gao Exp $
//

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

// method: ditfInit
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// this method creates lookup tables for the decimation in time
// frequency algorithm
//
// this method is based on the DITF algorithm by Ali Saidi and implements the
// butterfly version of the algorithm rather than the recursive version
// for reasons of efficiency.
//
//   Ali Saidi, "Decimation-In-Time-Frequency FFT Algorithm", Proceedings of
//   ICASSP 1994, pp. 453-456, San Francisco, CA, April 1994.
//
boolean FourierTransform::ditfInit(long order_a) {
  
  // if previous implementation and order are the same as
  // the current ones and the current transform object
  // is valid, then we are done
  //
  if (is_valid_d && (prev_implementation_d == DITF) &&
      (prev_order_d == order_a)) {

    // exit gracefully
    //
    return true;
  }

  // refresh the previous implementation type and transform order
  //
  prev_implementation_d = DITF;
  prev_order_d = order_a;
  
  // assign the order
  //
  order_d = order_a;

  // after this method the twiddle factors will be valid
  //
  is_valid_d = true;
  
  // define local variables
  //
  double q = Integral::TWO_PI / order_d;
  double t;
  long i;
  
  // allocate space for lookup tables
  //
  wd_real_d.setLength(order_d);
  wd_imag_d.setLength(order_d);
  
  // compute the lookup table entries
  //
  for (i = 0; i < order_d; i++) {
    t = q * i;
    wd_real_d(i) = cos(t);
    wd_imag_d(i) = sin(t);
  }

  // allocate memory for the temporary workspace used
  //
  tempd_d.setLength(order_d);
  ditf_trans_factor_indices_d.setLength(order_d);
  ditf_indices_d.setLength(order_d >> 1);
  
  // exit gracefully
  //
  return true;     
}

// method: ditfReal
//
// arguments:
//  VectorFloat& output: (output) output data vector
//  const VectorFloat& input: (input) input data vector
//
// return: a boolean value indicating status
//
// this method implements a real decimation in time frequency algorithm
//
// this method is based on the DITF algorithm by Ali Saidi and implements the
// butterfly version of the algorithm rather than the recursive version
// for reasons of efficiency.
//
//   Ali Saidi, "Decimation-In-Time-Frequency FFT Algorithm", Proceedings of
//   ICASSP 1994, pp. 453-456, San Francisco, CA, April 1994.
//
boolean FourierTransform::ditfReal(VectorFloat& output_a,
				   const VectorFloat& input_a) {
  // declare local variables
  //
  long transition_stage = 0;
  long trans_butterfly_size = 0;
  long m, m1, p, l, i1, i2, i, ii, jj, j, ib, jb, kb, k, kk, temp;
  long n1, n2, n3, n4, n5, pn2;
  float real_twid, imag_twid, tempr, tempi;
  
  if (!ditfInit(order_d)) {
    return Error::handle(name(), L"ditfReal", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"ditfReal", ERR_ORDER, __FILE__, __LINE__);
  }
    
  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  // copy input data to temporary memory so that the input data is retained
  // after calculations.
  //
  for (k = 0; k < order_d; ++k) {
    output_a(k) = input_a(k);
    tempd_d(k) = 0;
  }
  
  // initialize indices for the butterflies
  //
  n1 = order_d;   
  n2 = order_d >> 1;
  
  // only do the DIT portion for orders greater than 4
  //
  if (order_d > (long)4) {
    
    // calculate the location of the transition stage
    //
    temp = m;
    transition_stage = temp >> 1;
    
    // do the first decimation-in-time twiddleless butterfly
    //
    for (j = 0; j < n2; ++j) {
      l = j + n2;
      tempr = output_a(j);
      tempi = tempd_d(j);
      output_a(j) = tempr + output_a(l);
      tempd_d(j) = tempi + tempd_d(l);
      output_a(l) = tempr - output_a(l);
      tempd_d(l) = tempi - tempd_d(l);
    }
    
    // initialize a counter for the internal dit bit-reversals
    //
    n3 = 1;
	
    // do the remaining dit butterflies
    //
    for (i = 1; i < transition_stage; ++i) {
      n1 = n2;
      n2 = n2 >> 1;
      
      // procedure for bit reversal
      //
      jb = 0;
      n3 = n3 << 1;
      temp = n3;
      n5 = temp >> 1;
      n4 = n3 - 1;
      
      for (ib = 0; ib < n3; ++ib) {
	ditf_indices_d(ib) = ib;
      }
      for (ib = 0; ib < n4 ; ++ib) {
	
	if (ib < jb) {
	  ditf_indices_d(jb) = ib;
	  ditf_indices_d(ib) = jb;
	}
	kb = n5;
	while (kb <= jb) {
	  jb = jb - kb;
	  kb = kb >> (long)1;
	}
	jb = jb + kb;
      }
	    
      // carry out the dit butterfly
      //
      p = 0;
      for (j = 0; j < n3 ; ++j) {
	if ( j != 0 ) {
	  real_twid = (float)(wd_real_d((long)(n2*ditf_indices_d(j))));
	  imag_twid = (float)(wd_imag_d((long)(n2*ditf_indices_d(j))));
	  
	  for (k = p; k < p + n2; ++k) {
	    l = k + n2;
	    tempr = (real_twid * output_a(l)) +
	      (imag_twid * tempd_d(l));
	    tempi = (real_twid * tempd_d(l)) -
	      (imag_twid * output_a(l));
	    output_a(l) = output_a(k) - tempr;
	    tempd_d(l) = tempd_d(k) - tempi;
	    output_a(k) = output_a(k) + tempr;
	    tempd_d(k) = tempd_d(k) + tempi;
	  }
	}

	else {
	  pn2 = p + n2;
	  for (k = p; k < pn2; ++k) {
	    l = k + n2;
	    tempr = output_a(l);
	    tempi = tempd_d(l);
	    output_a(l) = output_a(k) - tempr;
	    tempd_d(l) = tempd_d(k) - tempi;
	    output_a(k) = output_a(k) + tempr;
	    tempd_d(k) = tempd_d(k) + tempi;
	  }
	}
	p = p + n1;
      }
    }
    
    // find the transition stage factors
    //
	
    // procedure for bit reversal
    //
    n3 = n3 << 1;
    long jb = 0;
    n4 = n3 - 1;
    temp = n3;
    n5 = temp >> 1;
    for (ib = 0; ib < n3; ++ib) {
      ditf_indices_d(ib) = ib;
    }
    for (ib = 0; ib < n4 ; ++ib) {
	    
      if (ib < jb) {
	ditf_indices_d(jb) = ib;
	ditf_indices_d(ib) = jb;
      }
      
      kb = n5;
      while (kb <= jb) {
	jb = jb - kb;
	kb = kb >> 1;
      }
      jb = jb + kb;
    }
    
	
    n5 = 0;
    trans_butterfly_size = order_d >> transition_stage;
    for (jj = 0; jj < n3; ++jj) {
      for (ii = 0; ii < trans_butterfly_size; ++ii) {
	ditf_trans_factor_indices_d(n5) = ii * ditf_indices_d(jj);
	++n5;
      }
    }
	
    // multiply by the transition stage factors.
    //
    for (i = 0; i < order_d; ++i) {
      real_twid = wd_real_d(ditf_trans_factor_indices_d(i));
      imag_twid = wd_imag_d(ditf_trans_factor_indices_d(i));
      tempr = output_a(i);
      tempi = tempd_d(i);
      output_a(i) = (real_twid * tempr) + (imag_twid * tempi);
      tempd_d(i) = (real_twid * tempi) - (imag_twid * tempr);
    }
  }
  
  // reset the stage counter
  //
  n2 = order_d >> transition_stage;
  
  // carry out the decimation-in-frequency portion
  //
  m1 = m - 1;
  for (i = transition_stage; i < m1; ++i) {
    n1 = n2;
    n2 = n2 >> 1;
    i1 = 0;
    i2 = (long)(order_d / n1);
    for (j = 0; j < n2; ++j) {
      if ( j != 0 ) {
	real_twid = wd_real_d(i1);
	imag_twid = wd_imag_d(i1);
	i1 = i1 + i2;
	
	for (k = j; k < order_d; k += n1) {
	  l = k + n2;
	  tempr = output_a(k) - output_a(l);
	  output_a(k) = output_a(k) + output_a(l);
	  tempi = tempd_d(k) - tempd_d(l);
	  tempd_d(k) = tempd_d(k) + tempd_d(l);
	  output_a(l) = (real_twid * tempr) + (imag_twid * tempi);
	  tempd_d(l) =  (real_twid * tempi) - (imag_twid * tempr);
	}
      }
      else {
	i1 = i1 + i2;
	for (k = j; k < order_d; k += n1) {
	  l = k + n2;
	  tempr = output_a(k) - output_a(l);
	  output_a(k) = output_a(k) + output_a(l);
	  tempi = tempd_d(k) - tempd_d(l);
	  tempd_d(k) = tempd_d(k) + tempd_d(l);
	  output_a(l) = tempr;
	  tempd_d(l) = tempi;
	}
      }
    }
  }
    
  // Calculate last twiddleless stage of the DIF portion
  //
  for (j = 0; j < order_d; j += 2) {
    l = j + 1;
    tempr = output_a(j);
    tempi = tempd_d(j);
    output_a(j) = tempr + output_a(l);