ft_11.cc

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

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

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

// method: fhInit
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// this method creates lookup tables for the fast hartley transform
//
// this code is based on the algorithm as described in the work by
// Ron Mayer whose code is available as free-ware. the algorithm is
// described in the paper:
//   Hsieh S. Hou,"The Fast Hartley Transform Algorithm",
//   IEEE Transactions on Computers, pp. 147-153, February 1987.
//
boolean FourierTransform::fhInit(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 == FAST_HARTLEY) &&
      (prev_order_d == order_a)) {

    // exit gracefully
    //
    return true;
  }

  // refresh the previous implementation type and transform order
  //
  prev_implementation_d = FAST_HARTLEY;
  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;
  
  // allocate memory for the temporary workspace used
  //
  tempd_d.setLength(order_d);
  wd_real_d.setLength(order_d);
  wd_imag_d.setLength(order_d);
    
  // exit gracefully
  //
  return true;     
}

// method: fhReal
//
// arguments:
//  VectorFloat& output: (output) output data vector
//  const VectorFloat& input: (input) input data vector
//  boolean isReal: (input) is the input data real?
//
// return: a boolean value indicating status
//
// this method implements a real fast hartley transform
//
// this code is based on the algorithm as described in the work by
// Ron Mayer whose code is available as free-ware. the algorithm is
// described in the paper:
//   Hsieh S. Hou,"The Fast Hartley Transform Algorithm",
//   IEEE Transactions on Computers, pp. 147-153, February 1987.
//
boolean FourierTransform::fhReal(VectorFloat& output_a,
				 const VectorFloat& input_a,
				 boolean isReal_a) {

  // declare local variables
  //
  float a, b;
  float c1, s1, s2, c2, s3, c3, s4, c4;
  float f0, g0, f1, g1, f2, g2, f3, g3;
  long i, j, kk, k, k1, k2, k3, k4, kx, m, temp;
  long fi, fn, gi;
  long N_d_2 = order_d >> 1;
  
  // initializes the variable in the lookup table
  //
  TRIG_VARS;
      
  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  if (!fhInit(order_d)) {
    return Error::handle(name(), L"fhReal", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"fhReal", 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) = 0;
    tempd_d(k) = input_a(k);
  }
  
  // perform bit reversal
  //
  for (k1 = 1, k2 = 0; k1 < order_d; ++k1) {
    for (k = order_d >> 1; (!((k2 ^= k) & k)); k >>= 1) {
      ;
    }
    if (k1 > k2) {
      a = tempd_d(k1);
      tempd_d(k1) = tempd_d(k2);
      tempd_d(k2) = a;
    }
  }
  
  // get the exponent of two
  //
  for (k = 0; (1 << k) < order_d; ++k);
  k &= (long)1;

  // twiddleless radix-2 computations when order is an even power of 2
  // completes the stage before the g's in fig.2
  //
  if (k == 0) {
    for (fi = 0, fn = order_d; fi < fn; fi += 4) {
      f1 = tempd_d(fi + 0) - tempd_d(fi + 1);
      f0 = tempd_d(fi + 0) + tempd_d(fi + 1);
      f3 = tempd_d(fi + 2) - tempd_d(fi + 3);
      f2 = tempd_d(fi + 2) + tempd_d(fi + 3);
      tempd_d(fi + 2) = f0 - f2;	
      tempd_d(fi + 0) = f0 + f2;
      tempd_d(fi + 3) = f1 - f3;	
      tempd_d(fi + 1) = f1 + f3;
    }
  }

  // twiddleless radix-2 computations when order is an odd power of 2
  // completes the stages till g's in fig. 1
  //
  else {
    for (fi = 0, fn = order_d, gi = fi + 1;
	 fi < fn; fi += 8, gi += 8) {
      c1 = tempd_d(fi + 0) - tempd_d(gi + 0);
      s1 = tempd_d(fi + 0) + tempd_d(gi + 0);
      c2 = tempd_d(fi + 2) - tempd_d(gi + 2);
      s2 = tempd_d(fi + 2) + tempd_d(gi + 2);
      c3 = tempd_d(fi + 4) - tempd_d(gi + 4);
      s3 = tempd_d(fi + 4) + tempd_d(gi + 4);
      c4 = tempd_d(fi + 6) - tempd_d(gi + 6);
      s4 = tempd_d(fi + 6) + tempd_d(gi + 6);
      f1 = s1 - s2;	
      f0 = s1 + s2;
      g1 = c1 - c2;	
      g0 = c1 + c2;
      f3 = s3 - s4;	
      f2 = s3 + s4;
      g3 = c4 * Integral::SQRT2;		
      g2 = c3 * Integral::SQRT2;
      tempd_d(fi + 4) = f0 - f2;
      tempd_d(fi + 0) = f0 + f2;
      tempd_d(fi + 6) = f1 - f3;
      tempd_d(fi + 2) = f1 + f3;
      tempd_d(gi + 4) = g0 - g2;
      tempd_d(gi + 0) = g0 + g2;
      tempd_d(gi + 6) = g1 - g3;
      tempd_d(gi + 2) = g1 + g3;
    }
  }
   
  // stop computations here if order is less than 16
  //
  if (order_d < (long)16) { 
        
    if (isReal_a) {

      // get fourier coefficients from hartley coefficients
      //
      for (i = 1, j = ((long)order_d - 1), k = N_d_2; i < k; ++i, --j) {
	a = tempd_d(i);
	b = tempd_d(j);
	tempd_d(j) = (a - b) * (double)0.5;
	tempd_d(i) = (a + b) * (double)0.5;
      }
    
      //  interlace the output and put into output_a
      //
      output_a(0) = tempd_d(0);
      output_a(1) = 0.0;

      for ( kk = 1; kk <= N_d_2; ++kk ) {
	output_a(2 * kk) = output_a(2 * order_d - 2 * kk ) = tempd_d(kk);
	output_a(2 * order_d - 2 * kk + 1 ) = tempd_d(order_d - kk);
	output_a(2 * kk + 1) = -tempd_d(order_d - kk);
      }
      output_a((long)order_d + 1) = 0.0;
    }
    
    if (!isReal_a) {

      // get the output data
      //
      for (kk = 0; kk < order_d; ++kk) {
	output_a(kk) = tempd_d(kk);
      }
    }
    
    // exit gracefully
    //
    return true;
  }

  // if order is greater than 16, continue twiddled computations
  //

  // does the computations for the even indices
  //
  do {
    k  += 2;
    k1  = 1 << k;
    temp = k1;
    k2  = temp << 1;
    temp = k2;
    k4  = temp << 1;
    k3  = k1 + k2;
    temp = k1;
    kx  = temp >> 1;
    fi  = 0;
    gi  = fi + kx;
    fn  = order_d;
    do {
      f1 = tempd_d(fi)      - tempd_d(fi + k1);
      f0 = tempd_d(fi)      + tempd_d(fi + k1);
      f3 = tempd_d(fi + k2) - tempd_d(fi + k3);
      f2 = tempd_d(fi + k2) + tempd_d(fi + k3);
      tempd_d(fi + k2) = f0 - f2;
      tempd_d(fi)      = f0 + f2;
      tempd_d(fi + k3) = f1 - f3;
      tempd_d(fi + k1) = f1 + f3;
      g1 = tempd_d(gi) - tempd_d(gi + k1);
      g0 = tempd_d(gi) + tempd_d(gi + k1);
      g3 = (float)tempd_d(gi + k3) * Integral::SQRT2;
      g2 = (float)tempd_d(gi + k2) * Integral::SQRT2;
      tempd_d(gi + k2) = g0 - g2;
      tempd_d(gi)      = g0 + g2;
      tempd_d(gi + k3) = g1 - g3;			
      tempd_d(gi + k1) = g1 + g3;
      gi = gi + k4;
      fi = fi + k4;
    } while (fi < fn);
    
    TRIG_INIT(k, c1, s1);

    // does computations for odd indices, this includes multiplication by
    // the twiddles. note that these computations use some sort of look
    // ahead and avoids redundancy of computations.
    //
    for (i = 1; (long)i < (long)kx; ++i) {

      // uses the euler form of approximation to compute the cosines
      //
      TRIG_NEXT(k, c1, s1);
      c2 = (c1 * c1) - (s1 * s1);
      s2 = (c1 * s1) * 2;
      fn = order_d;
      fi = i;
      gi = k1 - i;
      do {
	b = s2 * tempd_d(fi + k1) - c2 * tempd_d(gi + k1);
	a = c2 * tempd_d(fi + k1) + s2 * tempd_d(gi + k1);
	f1 = tempd_d(fi) - a;
	f0 = tempd_d(fi) + a;
	g1 = tempd_d(gi) - b;
	g0 = tempd_d(gi) + b;
	b = s2 * tempd_d(fi + k3) - c2 * tempd_d(gi + k3);
	a = c2 * tempd_d(fi + k3) + s2 * tempd_d(gi + k3);
	f3 = tempd_d(fi + k2) - a;
	f2 = tempd_d(fi + k2) + a;
	g3 = tempd_d(gi + k2) - b;
	g2 = tempd_d(gi + k2) + b;
	b = s1 * f2  - c1 * g3;
	a = c1 * f2  + s1 * g3;
	tempd_d(fi + k2) = f0 - a;
	tempd_d(fi)      = f0 + a;
	tempd_d(gi + k3) = g1 - b;
	tempd_d(gi + k1) = g1 + b;
	b = c1 * g2 - s1 * f3;
	a = s1 * g2 + c1 * f3;
	tempd_d(gi + k2) = g0 - a;
	tempd_d(gi)      = g0 + a;
	tempd_d(fi + k3) = f1 - b;
	tempd_d(fi + k1) = f1 + b;
	gi = gi + k4;
	fi = fi + k4;
      } while (fi < fn);
    }
    TRIG_RESET(k, c1, s1);
  } while (k4 < order_d);
  
  a = b = (double)0;  

  if (!isReal_a) {
    
    // get the output data
    //
    for (kk = 0; kk < order_d; ++kk) {
      output_a(kk) = tempd_d(kk);
    }
  }

  if (isReal_a) {
    
    // get fourier coefficients from hartley coefficients
    //
    for (i = 1, j = ((long)order_d - 1), k = N_d_2; i < k; ++i, --j) {
      a = tempd_d(i);
      b = tempd_d(j);
      tempd_d(j) = (a - b) * (double)0.5;
      tempd_d(i) = (a + b) * (double)0.5;
    }
    
    //  interlace the output and put into output_a
    //
    output_a(0) = tempd_d(0);
    output_a(1) = (double)0;

    for (kk = 1; kk <= N_d_2; ++kk) {
      output_a(kk * 2) = output_a(2 * order_d - kk * 2 ) = tempd_d(kk);
      output_a(kk * 2 + 1) = -tempd_d(order_d - kk);
      output_a(2 * order_d - kk * 2 + 1 ) = tempd_d(order_d - kk);
    }
    output_a((long)order_d + 1) = 0.0;
  }

  // exit gracefully
  //
  return true;
}

// method: fhComplex
//
// 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 fast hartley transform
//
// this code is based on the algorithm as described in the work by
// Ron Mayer whose code is available as free-ware. the algorithm is
// described in the paper:
//   Hsieh S. Hou,"The Fast Hartley Transform Algorithm",
//   IEEE Transactions on Computers, pp. 147-153, February 1987.
//
boolean FourierTransform::fhComplex(VectorFloat& output_a,
				    const VectorFloat& input_a) {

  // declare local variables
  //
  float a, b, c, d;
  float q, r, s, t;
  long i, j, k, m, N_d_2;

  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  if (!fhInit(order_d)) {
    return Error::handle(name(), L"fhComplex", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"fhComplex", ERR_ORDER, __FILE__, __LINE__);
  }

  // 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) {
    wd_real_d(k) = input_a(k*2);
    wd_imag_d(k) = input_a(k*2 + 1);
  }

  // prepare data for FHT computations
  //
  N_d_2 = order_d >> 1;
  for (i = 1, j = (long)order_d - 1, k = N_d_2; i < k; ++i, --j) {
    a = wd_real_d(i);
    b = wd_real_d(j);
    q = a + b;
    r = a - b;
    c = wd_imag_d(i);
    d = wd_imag_d(j);
    s = c + d;
    t = c - d;
    wd_real_d(i) = (q + t) * 0.5;
    wd_real_d(j) = (q - t) * 0.5;
    wd_imag_d(i) = (s - r) * 0.5;
    wd_imag_d(j) = (s + r) * 0.5;
  }
  
  // compute hartley transform of the real and imaginary parts separately
  //
  fhReal(output_a, wd_real_d, false);
  fhReal(wd_real_d, wd_imag_d, false);

  // interlace output and put into output_d array
  //
  for (k = (long)order_d - 1; k >= 0; --k) {
    output_a(k * 2) = output_a(k);
    output_a(k * 2 + 1) = wd_real_d(k);
  }
  
  // exit gracefully
  //
  return true;
}