ft_13.cc

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

    tempd_d(j) = tempi + tempd_d(l);
    output_a(l) = tempr - output_a(l);
    tempd_d(l) = tempi - tempd_d(l);
  }
  
  // procedure for bit reversal
  //
  temp = 0;
  j = 0;
  n1 = (long)order_d - 1;
  n2 = order_d >> (long)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 = j - k;
      k = k >> (long)1;
    }
    j = j + k;
  }
  
  // interlace the output and store in the output_a array
  //
  for (kk = (long)order_d - 1; kk >= 0; --kk) {
    output_a(2 * kk + 1) = tempd_d(kk);
    output_a(2 * kk) = output_a(kk);    
  }
  
  // exit gracefully
  //
  return true;  
}

// method: ditfComplex
//
// 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 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::ditfComplex(VectorFloat& output_a,
				      const VectorFloat& input_a) {
  // declare local variables
  //
  long transition_stage = 0;
  long trans_butterfly_size = 0;
  long m, m1, 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"ditfComplex", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"ditfComplex", 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(2*k);
    tempd_d(k) = input_a(2*k + 1);
  }
  
  // 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 >> 1;
	}
	jb = jb + kb;
      }
	    
      // carry out the dit butterfly
      //
      long 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 << (long)1;
    long jb = 0;
    n4 = n3 - 1;
    temp = n3;
    n5 = temp >> (long)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);
    tempd_d(j) = tempi + tempd_d(l);
    output_a(l) = tempr - output_a(l);
    tempd_d(l) = tempi - tempd_d(l);
  }
  
  // procedure for bit reversal
  //
  temp = 0;
  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 = j - k;
      k = k >> 1;
    }
    j = j + k;
  }
  
  // interlace the output and store in the output_a array
  //
  for (kk = (long)order_d - 1; kk >= 0; --kk) {
    output_a(2 * kk + 1) = tempd_d(kk);
    output_a(2 * kk) = output_a(kk);    
  }
  
  // exit gracefully
  //
  return true; 
}