ft_14.cc

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

// file: $isip/class/algo/FourierTransform/ft_14.cc
// version: $Id: ft_14.cc,v 1.9 2002/07/08 20:41:34 parihar Exp $
//

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

// method: dct_real_cc
//
// arguments:
//  VectorFloat& coeff: (output) output dct coefficient array
//  const VectorFloat& data: (input) input data array
//
// return: a logical_1 value indicating status
//
// this method implements a discrete cosine transform
//
// this code very closely follows the algorithm described in:
//    H.Guo et.al,"The Quick Discrete Fourier Transform", IEEE Transactions
//    on Acoustics Speech and Signal Processing 1994, pp. 445-448, 1994.
//
// this implementation is adapted from the code written by G.A.Sitton at
// Rice University
//
boolean FourierTransform::dct_real_cc(VectorFloat& coeff_a,
				      const VectorFloat& data_a) {

  // declare local variables
  //
  float t0, t1, t2, t3, t4;
  long XX;
  long xx;
  long i, j, k, n, nq, n2, n2p1, temp;
  long input_temp, output_temp;
  
  // Increment level & start	
  //
  qf_ii_d++;
  n = (long)qf_ws_d(qf_nc_d + qf_ii_d);
  xx = (long)qf_ws_d(qf_ii_d);

  // Last recursion level?	
  //
  if (n == N2) {
    
    // Form even function
    //
    qf_ws_d(xx + 0) = (float)(data_a(qf_input_d + 0) + data_a(qf_input_d + 8));
    qf_ws_d(xx + 1) = (float)(data_a(qf_input_d + 1) + data_a(qf_input_d + 7));
    qf_ws_d(xx + 2) = (float)(data_a(qf_input_d + 2) + data_a(qf_input_d + 6));
    qf_ws_d(xx + 3) = (float)(data_a(qf_input_d + 3) + data_a(qf_input_d + 5));
  
    // Save even DCT coefficients
    //
    t0 = qf_ws_d(xx + 0) - data_a(qf_input_d + 4);
    t3 = qf_ws_d(xx + 0) + data_a(qf_input_d + 4);
    t2 = qf_ws_d(xx + 1) + qf_ws_d(xx + 3);
    t1 = (qf_ws_d(xx + 1) - qf_ws_d(xx + 3)) * sqrt2;
    t4 = t3 + qf_ws_d(xx + 2);
    
    coeff_a(qf_output_d + 0) = t4 + t2;
    coeff_a(qf_output_d + 2) = t0 + t1;
    coeff_a(qf_output_d + 4) = t3 - qf_ws_d(xx + 2);
    coeff_a(qf_output_d + 6) = t0 - t1;
    coeff_a(qf_output_d + 8) = t4 - t2;
    
    // Form twiddled odd function
    //
    qf_ws_d(xx + 0) = (data_a(qf_input_d + 0) - data_a(qf_input_d + 8)) * sec0;
    qf_ws_d(xx + 1) = (data_a(qf_input_d + 1) - data_a(qf_input_d + 7)) * sec1;
    qf_ws_d(xx + 2) = (data_a(qf_input_d + 2) - data_a(qf_input_d + 6)) * sec2;
    qf_ws_d(xx + 3) = (data_a(qf_input_d + 3) - data_a(qf_input_d + 5)) * sec3;
    
    // Save odd DCT coefficients
    //
    t1 = (qf_ws_d(xx + 1) - qf_ws_d(xx + 3)) * sqrt2;
    t2 = qf_ws_d(xx + 1) + qf_ws_d(xx + 3);
    t4 = qf_ws_d(xx + 0) + qf_ws_d(xx + 2);
    
    coeff_a(qf_output_d + 1) = t4 + (t2 + (t0 = (qf_ws_d(xx + 0) + t1)));
    coeff_a(qf_output_d + 3) = t0 + (t3 = (qf_ws_d(xx + 0) - qf_ws_d(xx + 2)));
    coeff_a(qf_output_d + 5) = t3 + (t1 = (qf_ws_d(xx + 0) - t1));
    coeff_a(qf_output_d + 7) = (t1 + t4) - t2; 
  }
  
  else {
    temp = n;
    n2 = temp >> (long)1;
    n2p1 = n2 + 1;
    nq = (long)qf_ws_d(qf_ic_d + qf_ii_d);  
      
    if ((k = (long)qf_ws_d(qf_mc_d + qf_ii_d)) > n2)
      k = n2;
    
    // Form even-odd functions
    //
    if (qf_nn_d <= n2p1) {
      for (i = j = 0; i < qf_nn_d; ++i, j += nq) {
	qf_ws_d(xx + i) = data_a(qf_input_d + i);
	qf_ws_d(xx + (i + n2p1)) = qf_ws_d(xx + i) * qf_ws_d(qf_sc_d + j);
      }
    }
    else {
      temp = n - qf_nn_d;
      for (i = j = 0; i <= temp; ++i, j += nq) {
	qf_ws_d(xx + i) = data_a(qf_input_d + i);
	qf_ws_d(xx + (i + n2p1)) = qf_ws_d(xx + i) * qf_ws_d(qf_sc_d + j);
      }
      
      for (; i < n2; ++i, j += nq) {
	qf_ws_d(xx + i) = data_a(qf_input_d + i) + data_a(qf_input_d + n - i);
	qf_ws_d(xx + i + n2p1) = (data_a(qf_input_d + i) -
			       data_a(qf_input_d + n - i)) * qf_ws_d(qf_sc_d + j);
      }
      
      qf_ws_d(xx + i) = data_a(qf_input_d + i);
      qf_ws_d(xx + (i + n2p1)) = 0.0;
    }
    
    // Do recursive DCTs
    //
    XX = xx + (n + 2);
    input_temp = qf_input_d;
    output_temp = qf_output_d;
    qf_input_d = xx;
    qf_output_d = XX;
    dct_real_cc(qf_ws_d, qf_ws_d);

    qf_input_d = xx + n2p1;
    qf_output_d = XX + n2p1;
    dct_real_cc(qf_ws_d, qf_ws_d);

    qf_input_d = input_temp;
    qf_output_d = output_temp;
    
    // Save DCT coefficients
    //
    t0 = qf_ws_d(XX + n2p1);
    
    for (i = 0; i < k; i += 2) {
      j = i + i;
      coeff_a(qf_output_d + j) = qf_ws_d(XX + i);
      t1 = qf_ws_d(XX + (i + n2p1) + 1);
      coeff_a(qf_output_d + j + 1) = t0 + t1;
      coeff_a(qf_output_d + j + 2) = qf_ws_d(XX + (i + 1));
      t0 = qf_ws_d(XX + (i + n2p1) + 2);
      coeff_a(qf_output_d + j + 3) = t1 + t0;
    }
    coeff_a(qf_output_d + i + i) = qf_ws_d(XX + i);
  }

  // Decrement level & exit
  //
  qf_ii_d--;
  
  // exit gracefully
  //
  return true;
}

// method: dst_real_cc
//
// arguments:
//  VectorFloat& coeff: (output) dst coefficient array
//  const VectorFloat& data: (input) data array
//
// return: a logical_1 value indicating status
//
// this method implements a discrete sine transform
//
// this code very closely follows the algorithm described in:
//    H.Guo et.al,"The Quick Discrete Fourier Transform", IEEE Transactions
//    on Acoustics Speech and Signal Processing 1994, pp. 445-448, 1994.
//
// this implementation is adapted from the code written by G.A.Sitton at
// Rice University
//
boolean FourierTransform::dst_real_cc(VectorFloat& coeff_a,
				      const VectorFloat& data_a) {
  // declare local variables
  //
  float  t0, t1, t2;
  long XX;
  long xx ;
  long i, j, k, n, nq, n2, n2m1, temp;
  long input_temp, output_temp;
  
  // Increment level & start	
  //
  qf_ii_d++;
  n = (long)qf_ws_d(qf_nc_d + qf_ii_d);
  xx = (long)qf_ws_d(qf_ii_d);
  
  // Last recursion level ?	
  //
  if (n == N2) {

    // Form odd function; eqn (24)
    //
    qf_ws_d(xx + 0) = data_a(qf_input_d + 0) - data_a(qf_input_d + 6);
    qf_ws_d(xx + 1) = data_a(qf_input_d + 1) - data_a(qf_input_d + 5);
    qf_ws_d(xx + 2) = data_a(qf_input_d + 2) - data_a(qf_input_d + 4);
    
    // Save odd DST coefficients
    //
    t0 = (qf_ws_d(xx + 0) + qf_ws_d(xx + 2)) * sqrt2;
    
    coeff_a(qf_output_d + 1) = t0 + qf_ws_d(xx + 1);
    coeff_a(qf_output_d + 3) = qf_ws_d(xx + 0) - qf_ws_d(xx + 2) ;
    coeff_a(qf_output_d + 5) = t0 - qf_ws_d(xx + 1);
    
    // Form twiddled even function; eqn (24)
    //
    qf_ws_d(xx + 0) = (data_a(qf_input_d + 0) + data_a(qf_input_d + 6)) * sec1;
    qf_ws_d(xx + 1) = (data_a(qf_input_d + 1) + data_a(qf_input_d + 5)) * sec2;
    qf_ws_d(xx + 2) = (data_a(qf_input_d + 2) + data_a(qf_input_d + 4)) * sec3;
    
    // Save even DST coefficients
    //
    t0 = (qf_ws_d(xx + 0) + qf_ws_d(xx + 2)) * sqrt2;
    t1 = t0 + qf_ws_d(xx + 1);
    t2 = qf_ws_d(xx + 0) - qf_ws_d(xx + 2);
    coeff_a(qf_output_d + 0) = t1 + data_a(qf_input_d + 3);
    coeff_a(qf_output_d + 2) = t1 + t2 - data_a(qf_input_d + 3);
    t0 = t0 - qf_ws_d(xx + 1);
    coeff_a(qf_output_d + 4) = t0 + t2  + data_a(qf_input_d + 3);
    coeff_a(qf_output_d + 6) = t0 - data_a(qf_input_d + 3);
  }
  else {
    temp = n;
    n2 = temp >> 1;
    n2m1 = n2 - 1;
    nq = qf_ws_d(qf_ic_d + qf_ii_d);
    
    if ((k = qf_ws_d(qf_mc_d + qf_ii_d)) > n2m1)
      k = n2m1;
    
    // Form even-odd functions
    //
    if (qf_nn_d <= n2) {
      temp = qf_nn_d - 1;
      for (i = 0, j = nq; i < temp; ++i, j += nq) {
	qf_ws_d(xx + i) = data_a(qf_input_d + i);
	qf_ws_d(xx + i + n2m1) = qf_ws_d(xx + i) * qf_ws_d(qf_sc_d + j);
      }
    }
    else  {
      temp = n - qf_nn_d;
      for (i = 0, j = nq; i < temp; ++i, j += nq) {
	qf_ws_d(xx + i) = data_a(qf_input_d + i);
	qf_ws_d(xx + i + n2m1) = qf_ws_d(xx + i) * qf_ws_d(qf_sc_d + j);
      }
      
      for (; i < n2m1; ++i, j += nq) {
	qf_ws_d(xx + i) = data_a(qf_input_d + i) - data_a(qf_input_d + (n - 2) - i);
	qf_ws_d(xx + i + n2m1) = (data_a(qf_input_d + i) +
			       data_a(qf_input_d + (n - 2) - i)) * qf_ws_d(qf_sc_d + j);
      }
    }
    
    // Do recursive DSTs; similar to eqn (26)
    //
    XX = xx + (n - 2);
    input_temp = qf_input_d;
    output_temp = qf_output_d;
    qf_input_d = xx;
    qf_output_d = XX;
    dst_real_cc(qf_ws_d, qf_ws_d);

    qf_input_d = xx + n2m1;
    qf_output_d = XX + n2m1;
    dst_real_cc(qf_ws_d, qf_ws_d);

    qf_input_d = input_temp;
    qf_output_d = output_temp;
    
    // Save DST coefficients
    //
    coeff_a(qf_output_d + 1) = qf_ws_d(XX + 0);
    
    if (n2 >= qf_nn_d) {
      t0 = coeff_a(qf_output_d + 0) = qf_ws_d(XX + n2m1);
      
      for (i = 1; i < k; i += 2) {
	t1 = qf_ws_d(XX + i + n2m1);
	j = i + i;
	coeff_a(qf_output_d + j) = t1 + t0;
	coeff_a(qf_output_d + j + 1) = qf_ws_d(XX + i);
	t0 = qf_ws_d(XX + i + n2m1 + 1);
	coeff_a(qf_output_d + j + 2) = t0 + t1;
	coeff_a(qf_output_d + j + 3) = qf_ws_d(XX + i + 1);
      }
      coeff_a(qf_output_d + i + i) = t0;
    }
    else {
      t1 = qf_ws_d(XX + n2m1);
      t0 = data_a(qf_input_d + n2m1);
      coeff_a(qf_output_d + 0) = t1 + t0;
      
      for (i = 1; i < k; i += 2) {
	j = i + i;
	t2 = qf_ws_d(XX + i + n2m1);
	coeff_a(qf_output_d + j) = t2 + t1 - t0;
	coeff_a(qf_output_d + j + 1) = qf_ws_d(XX + i);
	t1 = qf_ws_d(XX + i + n2m1 + 1);
	coeff_a(qf_output_d + j + 2) = t1 + t2 + t0;
	coeff_a(qf_output_d + j + 3) = qf_ws_d(XX + i + 1);
      }
      
      coeff_a(qf_output_d + i + i) = t1 - t0;
    }
  }
  
  // Decrement level & exit
  //
  qf_ii_d--;
  
  // exit gracefully
  //
  return true;
}