ft_qf_0.cc

这是处理语音信号的程序

// file: $PDSP/class/fourier_transform/v3.0/ft_qf_0.cc
//

// system include files
//
#include <math.h>

// isip include files
//
#include "fourier_transform.h"
#include "ft_qf.h"

#include "fourier_transform_constants.h"

// method: qf_real_cc
//
// 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.
//
// arguments:
//  float_8* output_a: (output) output data array 
//                      length of output array will always be 2*N
//                      memory is allocated internally, not by the calling
//                      program.
//  float_8* input_a: (input) input data array
//                      length of input array is N
//                      memory is allocated by the calling
//                      program.
//
// return: logical_1 indicating status
//
// implements a real Quick Fourier Transform
//
logical_1 Fourier_transform::qf_real_cc(float_8* output_a,
						 float_8* input_a) {
  
  // declare local variables
  //
  float_8* y;
  float_8* xx;
  int_4	i, k, m, n2, temp;

  // initialize the lookup tables
  //
  if (qf_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"qf_real_cc",
		     (char_1*)"error initializing the algorithm");
    return ISIP_FALSE;
  }
  
  // check the order
  //
  if (is_power_cc((int_4&)m, (int_4)2) == ISIP_FALSE) {
    error_handler_cc((char_1*)"qf_real_cc",
		     (char_1*)"order must be a power of two");
    return ISIP_FALSE;
  }

  // copy input to temporary memory so that the input data is retained
  // after calculations.
  //
  for (k = 0; k < N_d; ++k) {
    qf_real_temp_d[k] = input_a[k];
  }

  // Set up constants	
  //
  if (N_d < MM)
    nn = MM;
  else
    nn = N_d;
  
  mm = N_d;
  n2 = N_d >> 1;
  ii = (int_4)0;
  xx = ws[(int_4)ii];
  
  for (i = 1; i < mc[0]; ++i)
    mc[i] = N_d;

  // form even function; eqn (15)
  //
  if (nn <= n2) {
    y = qf_real_temp_d;
  }
  
  else {
    xx[0] = qf_real_temp_d[0];
    temp = -nn + N_d;
    for (i = 1; i <= temp; i += 2){
      xx[i] = qf_real_temp_d[i];
      xx[i + 1] = qf_real_temp_d[i + 1];
    }
    
    xx[i] = qf_real_temp_d[i] + qf_real_temp_d[N_d - i];
    ++i;
    
    for (;i < n2; i += 2) {
      xx[i] = qf_real_temp_d[i] + qf_real_temp_d[N_d - i];
      xx[i+1] = qf_real_temp_d[i + 1] +  qf_real_temp_d[N_d - i - 1];
    }
    
    xx[i] = qf_real_temp_d[i];
    ++i;
    
    for (;i < MM; ++i)
      xx[i] = 0.0;
    
    y = xx;
  }
  
  // Do a real DCT
  //
  qf_real_dct_cc((float_8*)qf_real_coeff_d, (float_8*)y);
  
  //    Form negative odd function; eqn (16)
  //
  if (nn <= n2) {
    xx[1] = -qf_real_temp_d[1];
    
    for (i = 2; i < nn; i += 2) {
      xx[i] = -qf_real_temp_d[i];
      xx[i+1] = -qf_real_temp_d[i + 1];
    }
  }
  else {
    temp = -nn + N_d;
    for (i = 1; i <= temp; i += 2) {
      xx[i] = -qf_real_temp_d[i];
      xx[i+1] = -qf_real_temp_d[i + 1];
    }
    
    xx[i] = qf_real_temp_d[-i + N_d] - qf_real_temp_d[i];
    ++i;
    
      for (;i < n2; i += 2) {
	xx[i] = qf_real_temp_d[-i + N_d] - qf_real_temp_d[i];
	xx[i+1] = qf_real_temp_d[-(i+1) + N_d] - qf_real_temp_d[i+1];
      }
  }
  
  for (; i < MM; ++i)
    xx[i] = 0.0;
  
  // Do a real DST
  //
  qf_real_dst_cc((float_8*)&qf_imag_coeff_d[1], (float_8*)&xx[1]);
  qf_imag_coeff_d[0] = qf_imag_coeff_d[n2] = (float_8)0.0;
  
  // create symmetric part of the spectrum
  //
  for (k = 1; k <= N_d/2; ++k) {
    qf_imag_coeff_d[-k + N_d ] = -qf_imag_coeff_d[k];
    qf_real_coeff_d[- k + N_d] = qf_real_coeff_d[k];
  }
  qf_imag_coeff_d[N_d/2] = 0.0;
  
  // interlace output into output_a
  //
  for (k = N_d - 1; k >= 0; --k) {
    output_a[k*2] = qf_real_coeff_d[k];
    output_a[k*2+1] = qf_imag_coeff_d[k];
  }

  // deallocate memory
  //
  qf_close_cc();
  
  // exit gracefully
  //
  return ISIP_TRUE;
}