ft_19.cc

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

// file: $isip/class/algo/FourierTransform/ft_19.cc
// version: $Id: ft_19.cc,v 1.2 2002/05/31 21:57:29 picone Exp $
//

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

// method: triInit
//
// arguments:
//  long order: (input) new order
//
// return: a boolean value indicating status
//
// creates lookup tables for sin and cosine terms. note that space is
// allocated in this method only when the order changes.
//
boolean FourierTransform::triInit(long order_a) {

  return Error::handle(name(), L"triInit", ERR, __FILE__, __LINE__);
}

// method: triReal
//
// 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 discrete fourier transform using a
// table lookup. output data vector length of output vector will
// always be 2*order_d memory is allocated internally, not by the
// calling program. 
//
boolean FourierTransform::triReal(VectorFloat& output_a,
				 const VectorFloat& input_a) {

  // assign the order
  //  L: input data length, N: DFT points, A: N*L matrix
  //
  long L = input_a.length();
  long N = (olen_d == DEF_LENGTH) ? L : (long)olen_d;

  order_d = (direction_d == FORWARD) ? N : L;

  // declare local variables
  //
  long n, m, index;
  long order = (long)order_d, sym_pos, two_order = order * 2;
  float sumr, sumi;
  double scale_factor = Integral::TWO_PI / order, t;
  VectorFloat wd_cos(two_order), wd_sin(two_order);
  
  // allocate memory for output vector
  //
  output_a.setLength(order * 2);

  // initialize the look-up table
  //
  wd_cos(0) = 1;
  wd_cos(1) = cos(Integral::TWO_PI / order);
  wd_sin(0) = 0;
  wd_sin(1) = sin(Integral::TWO_PI / order);
  for (m = 2; m < two_order; m++) {
    t = scale_factor * m;
    wd_cos(m) = cos(t);
    wd_sin(m) = sin(t);
    wd_cos(two_order - m) = wd_cos(m);
    wd_sin(two_order - m) = -wd_sin(m);
  }
  
  // compute the first coefficient when k = 0
  //
  for (n = 0, sumr = 0; n < L; n++) {
    sumr += input_a(n);
  }
  output_a(0) = sumr;
  output_a(1) = 0;
  
  // loop over all frequency samples
  //
  long half_order = order / 2;
  for (long k = 1; k <= half_order; k++) {
    
    // declare local accumulators: n = 0
    //
    sumr = input_a(0);
    sumi = 0;

    // loop over the data
    //
    for (n = 1, m = k * n; n < L; n++, m += k) {
      
      // compute the real an imaginary parts
      //
      index = m % two_order;
      sumr += wd_cos(index) * input_a(n);
      sumi -= wd_sin(index) * input_a(n);
    }

    // put the output data in an interlaced format
    //
    output_a(2 * k) = sumr;
    output_a(2 * k + 1) = sumi;
    sym_pos = order - k;
    output_a(2 * sym_pos) = sumr;
    output_a(2 * sym_pos + 1) = -sumi;
  }
  
  // exit gracefully
  //
  return true;
}

// method: triComplex
//
// 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 discrete fourier transform
//
boolean FourierTransform::triComplex(VectorFloat& output_a,
				     const VectorFloat& input_a) {
  // assign the order
  //  L: input data length, N: DFT points, A: N*L matrix
  //
  long L = input_a.length() / 2;
  long N = (olen_d == DEF_LENGTH) ? L : (long)olen_d;

  order_d = (direction_d == FORWARD) ? N : L;

  // declare local variables
  //
  long n, m, index;
  long order = (long)order_d, two_order = order * 2;
  float sumr, sumi;
  double scale_factor = Integral::TWO_PI / order, t;
  VectorFloat wd_cos(two_order), wd_sin(two_order);
  
  // allocate memory for output vector
  //
  output_a.setLength(order * 2);

  // initialize the look-up table
  //
  wd_cos(0) = 1;
  wd_cos(1) = cos(Integral::TWO_PI / order);
  wd_sin(0) = 0;
  wd_sin(1) = sin(Integral::TWO_PI / order);
  for (m = 2; m < two_order; m++) {
    t = scale_factor * m;
    wd_cos(m) = cos(t);
    wd_sin(m) = sin(t);
    wd_cos(two_order - m) = wd_cos(m);
    wd_sin(two_order - m) = -wd_sin(m);
  }
  
  // loop over all frequency samples
  //
  for (long k = 0; k < order; k++) {
    
    // declare local accumulators: n = 0
    //
    sumr = 0;
    sumi = 0;

    // loop over the data
    //
    for (n = 0, m = k * n; n < L; n++, m += k) {
      
      // compute the real an imaginary parts
      //
      index = m % two_order;
      sumr += ((wd_cos(index) * input_a(2 * n)) +
	       (wd_sin(index) * input_a(2 * n + 1)));
      sumi += ((wd_cos(index) * input_a(2 * n + 1)) -
	       (wd_sin(index) * input_a(2 * n)));
    }

    // put the output data in an interlaced format
    //
    output_a(2 * k) = sumr;
    output_a(2 * k + 1) = sumi;
  }
  
  // exit gracefully
  //
  return true;
}