ft_09.cc

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

// file: $isip/class/algo/FourierTransform/ft_09.cc
// version: $Id: ft_09.cc,v 1.6 2002/01/31 21:41:46 zheng Exp $
//

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

// method: rad2Init
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// this method creates lookup tables for the radix-2 transform
//
// this method performs a radix-2 cooley-tukey decimation-in-frequency
// algorithm. the code follows the exact structure of that in the book:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 731-732, 1992.
//
boolean FourierTransform::rad2Init(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 == RADIX_2) &&
      (prev_order_d == order_a)) {

    // exit gracefully
    //
    return true;
  }

  // refresh the previous implementation type and transform order
  //
  prev_implementation_d = RADIX_2;
  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 wn = Integral::TWO_PI / order_d;
  double phi;
  
  // allocate space for lookup tables
  //
  wd_real_d.setLength(order_d);
  wd_imag_d.setLength(order_d);
  
  // compute the lookup table entries
  //
  for (long i = 0; i < order_d; i++) {
    phi = wn * i;
    wd_real_d(i) = cos(phi);
    wd_imag_d(i) = sin(phi);
  }

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

// method: rad2Real
//
// 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 radix-2 transform
//
// this method performs a radix-2 cooley-tukey decimation-in-frequency
// algorithm.  the code follows the exact structure of that in the book:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 731-732, 1992.
//
boolean FourierTransform::rad2Real(VectorFloat& output_a,
				   const VectorFloat& input_a) {

  // declare local variables
  //
  long i, j, k, l, m, i1, i2;
  long n1, n2;
  float real_twid, imag_twid;
  float tempr, tempi;
  
  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  if (!rad2Init(order_d)) {
    return Error::handle(name(), L"rad2Real", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, (long)2, order_d)) {
    return Error::handle(name(), L"rad2Real", ERR_ORDER, __FILE__, __LINE__);
  }
  
  // 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;
  }
  
  // actual dif radix-2 computation
  //
  n2 = order_d;

  // looping  through the first  stage
  //
  n1 = n2;
  n2 = n2 >> 1;
  i1 = 0;
  i2 = (long)((order_d) / n1);
  
  for (j = 0; j < n2; ++j) {
    
    // loop through each butterfly
    //
    real_twid = wd_real_d(i1);
    imag_twid = wd_imag_d(i1);
    i1 = i1 + i2;
    
    for (k = j; k < order_d; k += n1) {
	
      // perform the butterfly computation
      //
      l = k + n2;
      tempr = output_a(k) - output_a(l);
      output_a(k) = output_a(k) + output_a(l);
      output_a(l) = (real_twid * tempr);
      tempd_d(l) = -(imag_twid * tempr);
    }
  }
  
  for (i = 1; i < m; ++i) {
    
    // looping  through each of the remaining stages
    //
    n1 = n2;
    n2 = n2 >> 1;
    i1 = 0;
    i2 = (long)(order_d / n1);
    
    for (j = 0; j < n2; ++j) {

      // loop through each butterfly
      //
      real_twid = wd_real_d(i1);
      imag_twid = wd_imag_d(i1);
      i1 = i1 + i2;
      
      for (k = j; k < order_d; k += n1) {

	// perform the butterfly computation; computation within square braces
	// in eqn (9.3.34) multiplied by the twiddle factor
	//
	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);
      }
    }
  }
  
  // procedure for bit reversal
  //
  j = 0;
  n1 = (long)order_d - 1;
  n2 = order_d >> 1;
  
  for (i = 0; i < n1 ; ++i) {
    if (i < j) {
      tempr = output_a(j);
      output_a(j) = output_a(i);
      output_a(i) = tempr;
      tempr = tempd_d(j);
      tempd_d(j) = tempd_d(i);
      tempd_d(i) = tempr;
    }
    k = n2;
    while (k <= j) {
      j -= k;
      k = k >> 1;
    }
    j += k;
  }
  
  // interlace the output and store in the output array
  //
  for (k = (long)order_d - 1; k >= 0; --k) {
    output_a(k * 2 + 1) = tempd_d(k);
    output_a(k * 2) = output_a(k);    
  }
 
  // exit gracefully
  //
  return true;
}

// method: rad2Complex
//
// arguments:
//  VectorFloat& output: (output) output data vector
//  const VectorFloat& input: (input) input data vector
//
// return: a boolean value indicating status
//
// this method implements a complex radix-2 transform
//
// this method performs a radix-2 cooley-tukey decimation-in-frequency
// algorithm.  the code follows the exact structure of that in the book:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 731-732, 1992.
//
boolean FourierTransform::rad2Complex(VectorFloat& output_a,
				      const VectorFloat& input_a) {

  // declare local variables
  //
  long i, j, k, l, m, i1, i2;
  long n1, n2;
  float real_twid, imag_twid;
  float tempr, tempi, temp;
  
  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  if (!rad2Init(order_d)) {
    return Error::handle(name(), L"rad2Complex", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"rad2Complex", ERR_ORDER, __FILE__, __LINE__);
  }
 
  // copy input data to temporary memory so that the input data is retained
  // after calculations. note that the complex input data is in an interlaced
  // format. after this operation, output_a will contain the real part of the
  // input and tempd_d will contain the imaginary part.
  //
  for (k = 0; k < order_d; ++k) {
    output_a(k) = input_a(k * 2);
    tempd_d(k) = input_a(k * 2 + 1);
  }
  
  // actual dif radix-2 computation
  //
  n2 = order_d;  
  for (i = 0; i < m; ++i) {

    // looping  through each stage
    //
    n1 = n2;
    n2 = n2 >> 1;
    i1 = 0;
    i2 = (long)(order_d / n1);
    
    for (j = 0; j < n2; ++j) {

      // loop through each butterfly
      //
      real_twid = wd_real_d(i1);
      imag_twid = wd_imag_d(i1);
      i1 = i1 + i2;
      
      for (k = j; k < order_d; k += n1) {

	// perform the butterfly computation; computation within square braces
	// in eqn (9.3.34) multiplied by the twiddle factor
	//
	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);
      }
    }
  }

  // procedure for bit reversal
  //
  j = 0;
  n1 = (long)order_d - 1;
  n2 = order_d >> 1;

  for (i = 0; i < n1 ; ++i) {
    if (i < j) {
      temp = output_a(j);
      output_a(j) = output_a(i);
      output_a(i) = temp;
      temp = tempd_d(j);
      tempd_d(j) = tempd_d(i);
      tempd_d(i) = temp;
    }
    k = n2;
    while (k <= j) {
	j -= k;
	k = k >> 1;
    }
    j += k;
  }

  // interlace the output and store in the output_a array
  //
  for (k = (long)order_d - 1; k >= 0; --k) {
    output_a(k * 2 + 1) = tempd_d(k);
    output_a(k * 2) = output_a(k);
  }
    
  // exit gracefully
  //
  return true;
}