ft_10.cc

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

// file: $isip/class/algo/FourierTransform/ft_10.cc
// version: $Id: ft_10.cc,v 1.6 2002/01/31 21:41:46 zheng Exp $
//
 
// isip include files
//
#include "FourierTransform.h"

// method: rad4Init
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// creates lookup tables for the radix-4 transform
//
// this method performs a standard radix-4 decimation-in-frequency algorithm
// the code follows the exact structure of that in the book :
//
//  C.S. Burrus and T.W. Parks, "DFT/FFT and Convolution Algorithms - Theory
//  and Implementation", 2nd Edition, John Wiley and Sons, New York, NY, 1985.
// 
boolean FourierTransform::rad4Init(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_4) &&
      (prev_order_d == order_a)) {

    // exit gracefully
    //
    return true;
  }

  // refresh the previous implementation type and transform order
  //
  prev_implementation_d = RADIX_4;
  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 = (double)(Integral::TWO_PI / order_d);
  double t;
  
  // 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++) {
    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);
  
  // exit gracefully
  //
  return true;
}

// method: rad4Real
//
// 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-4 transform
//
// this method performs a standard radix-4 decimation-in-frequency algorithm
// the code follows the exact structure of that in the book :
//
//  C.S. Burrus and T.W. Parks, "DFT/FFT and Convolution Algorithms - Theory
//  and Implementation", 2nd Edition, John Wiley and Sons, New York, NY, 1985.
// 
boolean FourierTransform::rad4Real(VectorFloat& output_a,
				   const VectorFloat& input_a) {

  // declare local variables
  //
  long i, j, k, m, i1, i2, i3;
  long n1, n2, N_d_1, N_d_2;
  long a, b, c;
  float r1, r2, r3, r4, s1, s2, s3, s4, co1, co2, co3, si1, si2, si3, temp;

  // check the order
  //
  if (!isPower((long&)m, 4, order_d)) {
    return Error::handle(name(), L"rad4Real", ERR_ORDER, __FILE__, __LINE__);  
  }
  n2 = order_d;
 
  // create lookup table ( we use the same lookup table generator as
  // we used for rad2 computations )
  //
  if (!rad4Init(order_d)) {
      return Error::handle(name(), L"rad4Real", ERR_INIT, __FILE__, __LINE__);
  }

  // allocate memory for output vector
  //
  output_a.setLength((long)order_d * 2);
  
  // copy input data to  temporary memory and zero the output 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.0;
  }

  // main computation loop; looping through all the stages; there are
  // log4(order_d) stages.
  //
  
  // go through all the stages of the radix-4 butterfly.
  //
  for (k = 0; k < m; ++k) {
    n1 = n2;
    n2 = n2 >> 2;
    a = 0;

    // go through each butterfly in a stage; there are order_d/4 butterflies
    // in each stage.
    //
    for (j = 0; j < n2; j++) {
      b = a + a;
      c = a + b;
      co1 = wd_real_d(a);
      co2 = wd_real_d(b);
      co3 = wd_real_d(c);
      si1 = wd_imag_d(a);
      si2 = wd_imag_d(b);
      si3 = wd_imag_d(c);
      a = (long)((j + 1) * pow(4, k));
      
      // compute each radix-4 butterfly;
      //
      for (i = j; i < order_d; i += n1) {
	i1 = i + n2;
	i2 = i1 + n2;
	i3 = i2 + n2;
	r1 = output_a(i) + output_a(i2);
	r3 = output_a(i) - output_a(i2);
	s1 = tempd_d(i) + tempd_d(i2);
	s3 = tempd_d(i) - tempd_d(i2);
	r2 = output_a(i1) + output_a(i3);
	r4 = output_a(i1) - output_a(i3);
	s2 = tempd_d(i1) + tempd_d(i3);
	s4 = tempd_d(i1) - tempd_d(i3);
	output_a(i) = r1 + r2;
	r2 = r1 - r2;      
	r1 = r3 - s4;
	r3 = r3 + s4;
	tempd_d(i) = s1 + s2;
        s2 = s1 - s2;      
	s1 = s3 + r4;
	s3 = s3 - r4 ;
	output_a(i1) = (r3 * co1) + (s3 * si1);
	tempd_d(i1) = (s3 * co1) - (r3 * si1);
	output_a(i2) = (r2 * co2) + (s2 * si2);
	tempd_d(i2) = (s2 * co2) - (r2 * si2);
	output_a(i3) = (r1 * co3) + (s1 * si3);
	tempd_d(i3) = (s1 * co3) - (r1 * si3);
      }
    }
  }

  // procedure for bit reversal
  //
  temp = 0;
  j = (long)0;
  N_d_1 = (long)order_d - 1;
  N_d_2 = order_d >> 2;
  for (i = 0; i < N_d_1; ++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 = order_d >> 2;
    while (k * 3 <= j) {
      j -= k * 3;
      k = k >> 2;
    }
    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;
}

// method: rad4Complex
//
// 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-4 transform
//
// this method performs a standard radix-4 decimation-in-frequency algorithm
// the code follows the exact structure of that in the book :
//
//  C.S. Burrus and T.W. Parks, "DFT/FFT and Convolution Algorithms - Theory
//  and Implementation", 2nd Edition, John Wiley and Sons, New York, NY, 1985.
// 
boolean FourierTransform::rad4Complex(VectorFloat& output_a,
				      const VectorFloat& input_a) {
  
  // declare local variables
  //
  long i, j, k, m, i1, i2, i3;
  long n1, n2, N_d_1;
  long a, b, c;
  float r1, r2, r3, r4, s1, s2, s3, s4, co1, co2, co3, si1, si2, si3;
  
  // check the order
  //
  if (!isPower((long&)m, 4, order_d)) {
    return Error::handle(name(), L"rad4Complex", ERR_ORDER,
			 __FILE__, __LINE__);  
  }
  n2 = order_d;

  // create lookup table ( we use the same lookup table generator as
  // we used for rad2 computations )
  //
  if (!rad4Init(order_d)) {
      return Error::handle(name(), L"rad4Complex", ERR_INIT,
			   __FILE__, __LINE__);
  }

  // allocate memory for output vector
  //
  output_a.setLength((long)order_d * 2);
  
  // copy input data to temporary memory and zero the output so that
  // the input data is retained after calculations.
  //
  for (k = 0; k < order_d; ++k) {
    output_a(k) = input_a(k * 2);
    tempd_d(k) = input_a(k * 2 + 1);
  }
  
  // main computation loop; looping through all the stages; there are
  // log4(order_d) stages.
  //
  for (k = 0; k < m; ++k) {
    n1 = n2;
    n2  = n2 >> (long)2;
    a = 0;

    // go through each butterfly in a stage; there are order_d/4 butterflies
    // in each stage.
    //
    for (j = 0; j < n2; ++j) {
      b = a + a;
      c = a + b;
      co1 = wd_real_d(a);
      co2 = wd_real_d(b);
      co3 = wd_real_d(c);
      si1 = wd_imag_d(a);
      si2 = wd_imag_d(b);
      si3 = wd_imag_d(c);
      a = (long)((j + 1) * pow(4, k));
      
      // compute each radix-4 butterfly;
      //
      for (i = j; i < order_d; i += n1) {
	i1 = i + n2;
	i2 = i1 + n2;
	i3 = i2 + n2;
	r1 = output_a(i) + output_a(i2);
	r3 = output_a(i) - output_a(i2);
	s1 = tempd_d(i) + tempd_d(i2);
        s3 = tempd_d(i) - tempd_d(i2);
        r2 = output_a(i1) + output_a(i3);
        r4 = output_a(i1) - output_a(i3);
        s2 = tempd_d(i1) + tempd_d(i3);
        s4 = tempd_d(i1) - tempd_d(i3);
        output_a(i) = r1 + r2;
        r2 = r1 - r2;      
        r1 = r3 - s4;
        r3 = r3 + s4;     
        tempd_d(i) = s1 + s2;
        s2 = s1 - s2;      
        s1 = s3 + r4;
        s3 = s3 - r4;
	output_a(i1) = (r3 * co1) + (s3 * si1);
	tempd_d(i1) = (s3 * co1) - (r3 * si1);
	output_a(i2) = (r2 * co2) + (s2 * si2);
	tempd_d(i2) = (s2 * co2) - (r2 * si2);
	output_a(i3) = (r1 * co3) + (s1 * si3);
	tempd_d(i3) = (s1 * co3) - (r1 * si3);
      }
    }
  }
  
  // procedure for bit reversal
  //
  double temp = 0;
  j = 0;
  N_d_1 = (long)order_d - 1;
  for (i = 0; i < N_d_1; 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 = order_d >> 2;
    while (k * 3 <= j) {
      j -= k * 3;
      k  = k >> 2;
    }
    j += k;
  }

  N_d_1 = (long)order_d - 1;
  for (k = N_d_1; k >= 0; --k) {
    output_a(k * 2 + 1) = tempd_d(k);
    output_a(k * 2) = output_a(k);    
  }

  // exit gracefully
  //
  return true;
}