ft_08.cc

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

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

// method: srInit
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// this method creates lookup tables for sin and cosine terms
//
// note that lookup table computation is done in this method only when
// the order changes
//
boolean FourierTransform::srInit(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 == SPLIT_RADIX) && (prev_order_d == order_a)) {

    // exit gracefully
    //
    return true;
  }

  // refresh the previous implementation type and transform order
  //
  prev_implementation_d = SPLIT_RADIX;
  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;
  
  // 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: srReal
//
// arguments:
//  VectorFloat& output: (output) output data vector
//  const VectorFloat& input: (input) input data vector
//
// return: a boolean value indicating status
//
// this method is used to compute the real split-radix fft
//
// this code is based on the code by G.A.Sitton of the Rice
// University. the theory behind the implementation closely follows
// the description of the split-radix algorithm in:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 727-730, 1992.
//
boolean FourierTransform::srReal(VectorFloat& output_a,
				 const VectorFloat& input_a) {

  // declare local variables
  //
  float  t1, t2, t3, t4, t5, t6, t;
  float  cc1, ss1, cc3, ss3;
  long  i, j, k, l, m, order_2, inc;
  long  n1, n2, n4, n8;
  long  is, id, i1, i2, i3, i4, i5, i6, i7, i8;

  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  if (!srInit(order_d)) {
    return Error::handle(name(), L"srReal", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"srReal", 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++) {
    tempd_d(k) = input_a(k);
    output_a(k) = 0;
  }

  // bit reversal routine
  //
  n1 = (long)order_d - 1;
  n2 = (long)order_d >> 1;

  for (i = 0, j = 0; i < n1; i++) {
    if (i < j) {
      t = tempd_d(j);
      tempd_d(j) = tempd_d(i);
      tempd_d(i) = t;
    }

    long k = n2;
    while (k <= j) {
      j -= k;
      k = k >> 1;
    }
    j += k;
  }
  
  // First (twiddleless) L Butterflies; completes stage A in fig. 9.23   
  //
  is = 0;
  id = 4;
  n1 = (long)order_d - 1;
  
  do {
    for (i = is; i < n1; i += id) {
      t1 = tempd_d(i);
      i1 = i + 1;
      t2 = tempd_d(i1);
      tempd_d(i) = t1 + t2;
      tempd_d(i1) = t1 - t2;
    }
    is = (id << 1) - 2;
    id = id << 2;
  } while (is < n1); 

  //  Other (twiddled) L Butterflies  
  //
  n2 = 2;
  for (k = 2; k <= m; k++) {
    n4 = n2 >> 1;
    n8 = n4 >> 1;
    n2 = n2 << 1;
    is = 0;
    id = n2 << 1; 
    do {
      for (i = is; i < order_d; i += id) {
        i2 = i + n4;
        i3 = i2 + n4;
        i4 = i3 + n4;
        t3 = tempd_d(i3);
        t4 = tempd_d(i4);
        t1 = tempd_d(i);
        t2 = t4 + t3;
        tempd_d(i4) = t4 - t3;
        tempd_d(i3) = t1 - t2;
        tempd_d(i) = t1 + t2;
 
        if (n4 != 1) {
          i1 = i + n8;
          i2 = n8 + i2;
          i3 = n8 + i3;
          i4 = n8 + i4;
          t3 = tempd_d(i3);
          t4 = tempd_d(i4);
          t1 = (t3 + t4) * Integral::INV_SQRT2;
          t2 = (t3 - t4) * Integral::INV_SQRT2;
          t3 = tempd_d(i1);
          t4 = tempd_d(i2);
          tempd_d(i4) = t4 - t1;
          tempd_d(i3) = -(t4 + t1);
          tempd_d(i2) = t3 - t2;
          tempd_d(i1) = t3 + t2;
        }
      }
      
      is = (id << 1) - n2;
      id = id << 2;
    } while (is < n1);
    
    inc = order_d / n2;
    l = inc;
    for (j = 1; j < n8; j++) {
      cc1 = wd_real_d(l);
      ss1 = wd_imag_d(l);

      // i1 gives the twiddle factors which are multiples of 3 and are needed
      // at the B stage in fig. 9.23
      //
      i1 = l * 3;
      cc3 = wd_real_d(i1);
      ss3 = wd_imag_d(i1);
      l = (j + 1) * inc;
      is = 0;
      id = n2 << 1;
      do {
	for (i = is; i < order_d; i += id) {
	  i1 = i + j;
	  i2 = i1 +  n4;
	  i3 = i2 +  n4;
	  i4 = i3 +  n4;
	  i5 = i - j + n4;
	  i6 = i5 +  n4;
	  i7 = i6 +  n4;
	  i8 = i7 +  n4;
	  t3 = tempd_d(i3);
	  t4 = tempd_d(i7);
	  t1 = (t3 * cc1) + (t4 * ss1);
	  t2 = (t4 * cc1) - (t3 * ss1);
	  t5 = tempd_d(i4);
	  t6 = tempd_d(i8);
	  t3 = (t5 * cc3) + (t6 * ss3);
	  t4 = (t6 * cc3) - (t5 * ss3);
	  t5 = t1 + t3;
	  t3 = t1 - t3;
	  t6 = t2 + t4;
	  t4 = t2 - t4;
	  t1 = tempd_d(i6);
	  tempd_d(i8) = t1 + t6;
	  tempd_d(i3) = t6 - t1;
	  t2 = tempd_d(i2);
	  tempd_d(i4) = t2 - t3;
	  tempd_d(i7) = -(t2 + t3);
	  t1 = tempd_d(i1);
	  tempd_d(i1) = t1 + t5;
	  tempd_d(i6) = t1 - t5;
	  t5 = tempd_d(i5);
	  tempd_d(i2) = t5 + t4;
	  tempd_d(i5) = t5 - t4;
	}
	
	is = (id << 1) - n2;
	id = id << 2;
      } while (is < n1);
    }
  }
    
  // interlace the output and store in the output_a array
  //
  order_2 = (long)order_d / 2;
  for (k = 1; k <= order_2; k++) {
    output_a(k * 2) = output_a(2 * order_d - k * 2) = tempd_d(k);
    output_a(2 * order_d - k * 2 + 1 ) = -tempd_d(order_d - k);
    output_a(k * 2 + 1) = tempd_d(order_d - k);
  }
  
  // manually set values to special points in the output data
  // zero imaginary part of dc, and N/2 point
  //
  output_a(0) = tempd_d(0);
  output_a(1) = (double)0.0;
  output_a((long)order_d + 1) = (double)0.0;
  
  // exit gracefully
  //
  return true;
}

// method: srComplex
//
// arguments:
//  VectorFloat& output: (output) output data vector
//  const VectorFloat& input: (input) input data vector
//
// return: boolean indicating status of computation
//
// this method computes the fourier transform using the split-radix
// algorithm
//
// this code is based on the code by G.A.Sitton of the Rice
// University.  the theory behind the implementation closely follows
// the description of the split-radix algorithm in:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 727-730, 1992.
//
boolean FourierTransform::srComplex(VectorFloat& output_a,
                                    const VectorFloat& input_a) {
  
  // declare local variables
  //
  long	i, j, k, m, q, ii, i0, i1, i2, i3, is, id, temp;
  float cc1, ss1, cc3, ss3, r1, r2, r3, s1, s2, t;
  long n1, n2, n4;

  // allocate memory for output vector
  //
  output_a.setLength((long)order_d * 2);
  
  // initialize lookup tables and temporary memory
  //
  if (!srInit(order_d)) {
    return Error::handle(name(), L"srComplex", ERR_INIT, __FILE__, __LINE__);

  }

  // check if inout data is a power of 2
  //
  if (!isPower((long&)m, (long)2, order_d)) {
    return Error::handle(name(), L"srComplex", ERR_ORDER, __FILE__, __LINE__);
  }
  
  // copy input real data to temporary memory and input imaginary data
  // into output so that the input data is retained
  // after calculations.
  //
  for (k = 0; k < order_d; k++) {
    tempd_d(k) = input_a(k * 2);
    output_a(k) = input_a(k * 2 + 1);
  }
  
  // bit reversal routine
  //
  n1 = (long)order_d - 1;
  n2 = order_d >> 1;
  j = 0;
    
  for (i = 0; i < n1; i++) {
    if (i < j) {
      t = tempd_d(j);
      tempd_d(j) = tempd_d(i);
      tempd_d(i) = t;
      t = output_a(j);
      output_a(j) = output_a(i);
      output_a(i) = t;
    }
    k = n2;
    while (k <= j) {
      j -= k;
      k = k >> 1;
    }
    j += k;
  }
  
  // Length two butterflies	
  //
  is = 0;
  id = 1;
  n1 = (long)order_d - 1;
  
  do {
    id = id << 2;      
    for (i0 = is; i0 < order_d; i0 += id) {
      i1 = i0 + 1;
      
      r1 = tempd_d(i0);
      tempd_d(i0) = r1 + tempd_d(i1);
      tempd_d(i1) = r1 - tempd_d(i1);
	
      r1 = output_a(i0);
      output_a(i0) = r1 + output_a(i1);
      output_a(i1) = r1 - output_a(i1); 
      }
    temp = id - 1;
    is = temp << 1;
  } while (is < n1);
  
  // L shaped butterflies	
  //
  n2 = 2;
  for (k = 2; k <= m; k++) {
    temp = n2;
    n4 = temp >> 1;
    n2 = n2 << 1;
    is = 0;
    temp = n2;
    id = temp << 1;

    // Unity twiddle factor loops	
    //
    do {
      for (i0 = is; i0 < n1; i0 += id) {
	i1 = i0 + n4;
	i2 = i1 + n4;
	i3 = i2 + n4;
	r3 = tempd_d(i2) + tempd_d(i3);
	r2 = tempd_d(i2) - tempd_d(i3);
	r1 = output_a(i2) + output_a(i3);
	s2 = output_a(i2) - output_a(i3);
	
	tempd_d(i2) = tempd_d(i0) - r3;
	tempd_d(i0) = tempd_d(i0) + r3;
	tempd_d(i3) = tempd_d(i1) - s2;
	tempd_d(i1) = s2 + tempd_d(i1);
	
	output_a(i2) = output_a(i0) - r1;
	output_a(i0) = output_a(i0) + r1;
	output_a(i3) = output_a(i1) + r2;
	output_a(i1) = output_a(i1) - r2;
	
      } 
      is = (id << 1) - n2;
      id = id << 2; 
      
    } while (is < n1);
    
    q = ii = (long)(order_d / n2);
    
    // Non-trivial twiddle factor loops	
    //
    for (j = 1; j < n4; j++) {
      cc1 = wd_real_d(q);
      ss1 = wd_imag_d(q);
      i3 = q * 3;
      cc3 = wd_real_d(i3);
      ss3 = wd_imag_d(i3);
      
      is = j;
      id = n2 << 1;
      
      do {
	for (i0 = is; i0 < n1; i0 += id) {
	  i1 = i0 + n4;
	  i2 = i1 + n4;
	  i3 = i2 + n4;
	     
	  r1 = (tempd_d(i2) * cc1) + (output_a(i2) * ss1);
	  s1 = (output_a(i2) * cc1) - (tempd_d(i2) * ss1);
	  r2 = (tempd_d(i3) * cc3) + (output_a(i3) * ss3);
	  s2 = (output_a(i3) * cc3) - (tempd_d(i3) * ss3);
	  
	  r3 = r1 + r2;
	  r2 = r1 - r2;
	  r1 = s1 + s2;
	  s2 = s1 - s2;
	  
	  tempd_d(i2) = tempd_d(i0) - r3;
	  tempd_d(i0) = r3 + tempd_d(i0);
	  tempd_d(i3) = tempd_d(i1) - s2;
	  tempd_d(i1) = s2 + tempd_d(i1);
	  
	  output_a(i2) = output_a(i0) - r1;
	  output_a(i0) = r1 + output_a(i0);
	  output_a(i3) = output_a(i1) + r2;
	  output_a(i1) = output_a(i1) - r2;
	}
	is = (id << 1) - n2 + j;
	id = id << 2;
      } while (is < n1);
      q += ii;
    }
  }
  
  // interlace the output and store in the output_a array
  //
  for (k = (long)order_d - 1; k >= 0; k--) {
    output_a(k * 2 + 1) = output_a(k);
    output_a(k * 2) = tempd_d(k);    
  }

  // exit gracefully
  //
  return true;
}