ft_12.cc

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

// file: $isip/class/algo/ourierTransform/ft_12.cc
// version: $Id: ft_12.cc,v 1.7 2002/07/08 20:41:34 parihar Exp $
//

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

// method: qfInit
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// this method creates lookup tables for the quick fourier transform
//
// this code very closely follows the algorithm described in:
//    H.Guo et.al,"The Quick Discrete Fourier Transform",  Proceedings of
//    ICASSP 1994, pp. 453-456, SanFrancisco, CA, April 1994.
//
// this implementation is adapted from the code written by G.A.Sitton at
// Rice University
//
boolean FourierTransform::qfInit(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 == QF) &&
      (prev_order_d == order_a)) {

    // exit gracefully
    //
    return true;
  }

  // refresh the previous implementation type and transform order
  //
  prev_implementation_d = QF;
  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;

  double q;
  long i, m, n2, n2m1, temp;

  n2 = order_d >> 1;
  isPower((long&)m, 2, order_d);
  
  // allocate the QFT workspace
  //
  qf_ws_d.setLength(sizeof(long) * m + sizeof(long) * 6 * m +
		   sizeof(float) * 3 * (order_d + m));

  // seperate memory is allocated for real and complex computations
  // in this case since the qft real routine is called from the complex
  // routine
  //
  wd_real_d.setLength(order_d);
  wd_imag_d.setLength(order_d);
  tempd_d.setLength(order_d);

  // memory space specifically used only for the complex routine
  //
  qf_comp_real_coeff_d.setLength(2 * order_d);
  qf_comp_imag_coeff_d.setLength(2 * order_d);
  qf_comp_temp_d.setLength(order_d);
    
  // the following variables are used to do the indexing, recursion and
  // workspace traversing
  //
  qf_nc_d = m;
  qf_ic_d = qf_nc_d + m;
  qf_mc_d = qf_ic_d + m;
  qf_sc_d = qf_mc_d + m;
  qf_ws_d(0) = qf_sc_d + (n2 - 1);
  
  // Initialize tables	
  //
  qf_ws_d(qf_nc_d + 0) = order_d;
  qf_ws_d(qf_ic_d + 0) = 1;
  qf_ws_d(qf_mc_d + 0) = m;
  
  for (i = 1; i < m; ++i) {
    temp = qf_ws_d(qf_nc_d + (i - 1));
    qf_ws_d(qf_nc_d + i) = temp >>  (long)1;
    temp = qf_ws_d(qf_ic_d + (i - 1));
    qf_ws_d(qf_ic_d + i) = temp << (long)1;
    
    if (i > 1)
      qf_ws_d(i) = (float)qf_ws_d(i - 1) +
	(2 * (long)((float)qf_ws_d(qf_nc_d + (i - 1)) + 2));
    else
      qf_ws_d(i) = qf_ws_d(i - 1) + qf_ws_d(qf_nc_d + i) + 1;	
  }	
  
  // Compute scaled secants to be used to compute the even and odd functions
  // which are used in the DCT and DST computations.
  //
  q = Integral::PI / order_d;
  n2m1 = n2 - 1;
  for (i = 0; i < n2m1; ++i) {
    qf_ws_d(qf_sc_d + i) = 0.5 / cos(q * i);	
  }

  // exit gracefully
  //
  return true;     
}

// method: qfReal
//
// 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 quick fourier transform
//
// this code very closely follows the algorithm described in:
//    H.Guo et.al,"The Quick Discrete Fourier Transform",  Proceedings of
//    ICASSP 1994, pp. 453-456, SanFrancisco, CA, April 1994.
//
// this implementation is adapted from the code written by G.A.Sitton at
// Rice University
//
boolean FourierTransform::qfReal(VectorFloat& output_a,
				    const VectorFloat& input_a) {
  
  // declare local variables
  //
  long xx;
  long	i, k, m, n2, temp;
  
  if (!qfInit(order_d)) {
    return Error::handle(name(), L"qfReal", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"qfReal", ERR_ORDER, __FILE__, __LINE__);
  }

  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);
  
  // copy input 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);
  }

  // Set up constants	
  //
  if (order_d < (long)MM)
    qf_nn_d = MM;
  else
    qf_nn_d = order_d;
  
  qf_mm_d = order_d;
  n2 = order_d >> 1;
  qf_ii_d = 0;
  xx = qf_ws_d(qf_ii_d);
  
  for (i = 1; i < qf_ws_d(qf_mc_d + 0); ++i) {
    qf_ws_d(qf_mc_d + i) = order_d;
  }

  // form even function; eqn (15)
  //
  if (qf_nn_d <= n2) {

    // Do a real DCT
    //
    qf_input_d = 0;
    qf_output_d = 0;
    dct_real_cc(wd_real_d, tempd_d);
  }
  else {
    qf_ws_d(xx + 0) = tempd_d(0);
    temp = -qf_nn_d + order_d;
    for (i = 1; i <= temp; i += 2) {
      qf_ws_d(xx + i) = tempd_d(i);
      qf_ws_d(xx + (i + 1)) = tempd_d(i + 1);
    }
    
    qf_ws_d(xx + i) = tempd_d(i) + tempd_d(order_d - i);
    ++i;
    
    for (; i < n2; i += 2) {
      qf_ws_d(xx + (i)) = tempd_d(i) + tempd_d(order_d - i);
      qf_ws_d(xx + (i + 1)) = tempd_d(i + 1) +  tempd_d(order_d - i - 1);
    }
    
    qf_ws_d(xx + i) = tempd_d(i);
    i++;
    
    for (; i < MM; ++i) {
      qf_ws_d(xx + i) = 0.0;
    }
    
    // Do a real DCT
    //
    qf_input_d = xx;
    qf_output_d = 0;
    dct_real_cc(wd_real_d, qf_ws_d);
  }
  
  // Form negative odd function; eqn (16)
  //
  if (qf_nn_d <= n2) {
    qf_ws_d(xx + 1) = -tempd_d(1);
    
    for (i = 2; i < qf_nn_d; i += 2) {
      qf_ws_d(xx + i) = -tempd_d(i);
      qf_ws_d(xx + (i + 1)) = -tempd_d(i + 1);
    }
  }
  else {
    temp = -qf_nn_d + order_d;
    for (i = 1; i <= temp; i += 2) {
      qf_ws_d(xx + i) = -tempd_d(i);
      qf_ws_d(xx + (i + 1)) = -tempd_d(i + 1);
    }
    
    qf_ws_d(xx + i) = tempd_d(-i + order_d) - tempd_d(i);
    i++;
    
    for (; i < n2; i += 2) {
      qf_ws_d(xx + i) = tempd_d(-i + order_d) - tempd_d(i);
      qf_ws_d(xx + (i + 1)) = tempd_d(-(i + 1) + order_d) -
	tempd_d(i + 1);
    }
  }
  
  for (; i < MM; ++i)
    qf_ws_d(xx + i) = 0.0;
  
  // Do a real DST
  //
  qf_input_d = xx + 1;
  qf_output_d = 1;
  dst_real_cc(wd_imag_d, qf_ws_d);
  wd_imag_d(0) = wd_imag_d(n2) = (float)0.0;

  // create symmetric part of the spectrum
  //
  for (k = 1; k <= (long)order_d / 2; ++k) {
    wd_imag_d(-k + order_d ) = -wd_imag_d(k);
    wd_real_d(- k + order_d) = wd_real_d(k);
  }
  wd_imag_d((long)order_d / 2) = 0.0;
  
  // interlace output into output_a
  //
  for (k = (long)order_d - 1; k >= 0; --k) {
    output_a(k * 2) = wd_real_d(k);
    output_a(k * 2 + 1) = wd_imag_d(k);
  }

  // exit gracefully
  //
  return true;
}

// method: qfComplex
//
// 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 quick fourier transform
//
// this code very closely follows the algorithm described in:
//    H.Guo et.al,"The Quick Discrete Fourier Transform",  Proceedings of
//    ICASSP 1994, pp. 453-456, SanFrancisco, CA, April 1994.
//
// this implementation is adapted from the code written by G.A.Sitton at
// Rice University
//
boolean FourierTransform::qfComplex(VectorFloat& output_a,
				    const VectorFloat& input_a) {

  // declare local variables
  //
  long m;

  // allocate memory for output vector
  //
  output_a.setLength(2 * order_d);

  if (!qfInit(order_d)) {
    return Error::handle(name(), L"qfComplex", ERR_INIT, __FILE__, __LINE__);
  }
  
  if (!isPower((long&)m, 2, order_d)) {
    return Error::handle(name(), L"qfComplex", ERR_ORDER, __FILE__, __LINE__);
  }

  // copy input real data to temporary memory so that the input data is
  // retained after calculations.
  //
  for (long i = 0; i < order_d; i++) {
    qf_comp_temp_d(i) = input_a(i * 2);
  }

  // compute packed QFT of the real parts; as stated on page 447
  //
  qfReal(qf_comp_real_coeff_d, qf_comp_temp_d);

  // copy input imaginary data to temporary memory so that the input data is
  // retained after calculations.
  //
  for (long i = 0; i < order_d; i++) {
    qf_comp_temp_d(i) = input_a((i * 2) + 1);
  }

  // compute packed QFT of the imaginary parts; as stated on page 447
  //
  qfReal(qf_comp_imag_coeff_d, qf_comp_temp_d);

  // interlace the output data into output_a
  //
  for (long i = 0; i < order_d; ++i) {
    output_a(i*2) = qf_comp_real_coeff_d(i*2) - qf_comp_imag_coeff_d(i*2+1);
    output_a(i*2+1) = qf_comp_real_coeff_d(i*2+1) + qf_comp_imag_coeff_d(i*2);
  }

  // exit gracefully
  //
  return true;
}