ft_fh_1.cc

这是处理语音信号的程序

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

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

#include "fourier_transform_constants.h"

// definition required to initialize some functions and arrays specific
// to the fht provided by ron mayer
//
#define FT_FHT_GOOD_TRIG

// method: fht_real_cc
//
// this code is based on the algorithm as described in the work by
// ron mayer whose code is available as a free-ware. the algorithm is
// described in the paper:
//   Hsieh S. Hou,"The Fast Hartley Transform Algorithm",
//   IEEE Transactions on Computers, pp. 147-153, February 1987.
//
// 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 2*N
//			input data memory should be allocated by the
//			calling program.
//
// return: a logical_1 value indicating status
//
// implements a complex Fast Hartley transform
//
logical_1 Fourier_transform::fh_complex_cc(float_8* output_a,
						    float_8* input_a) {
  // declare local variables
  //
  float_8 a, b, c, d;
  int_4 i, j, k, m, N_d_2;
  float_8 q,r,s,t;

  // initialize lookup tables and temporary memory
  //
  if (fh_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"fh_complex_cc",
		     (char_1*)"error generating the lookup table");
    return ISIP_FALSE;
  }
  USE_STATIC;
  
  // check if inout data is a power of 2
  //
  if (is_power_cc((int_4&)m, (int_4)2) == ISIP_FALSE) {
    error_handler_cc((char_1*)"fh_complex_cc",
		     (char_1*)"order must be a power of 2");
    return ISIP_FALSE;
  }
  
  // copy input data to  temporary memory and zero the output so that
  // the input data is retained after calculations.
  //
  for (k = 0; k < N_d; ++k) {
    fh_temp_output_d[k] = input_a[k*2];
    fh_temp_input_d[k] = input_a[k*2 + 1];
  }

  // prepare data for FHT computations
  //
  N_d_2 = N_d >> 1;
  for (i = 1,j = N_d - 1, k = N_d_2; i < k; ++i, --j) {
    a = fh_temp_output_d[i];
    b = fh_temp_output_d[j];
    q = a + b;
    r = a - b;
    c = fh_temp_input_d[i];
    d = fh_temp_input_d[j];
    s = c + d;
    t = c - d;
    fh_temp_output_d[i] = (q+t) * 0.5;
    fh_temp_output_d[j] = (q-t) * 0.5;
    fh_temp_input_d[i] = (s-r) * 0.5;
    fh_temp_input_d[j] = (s+r) * 0.5;
  }
  
  // compute hartley transform of the real and imaginary parts separately
  //
  fh_real_cc((float_8*)output_a,(float_8*)fh_temp_output_d);
  fh_real_cc((float_8*)fh_temp_d,(float_8*)fh_temp_input_d);
  
  // interlace output and put into output_d array
  //
  for (k = N_d - 1; k >= 0; --k) {
    output_a[k*2] = output_a[k];
    output_a[k*2+1] = fh_temp_d[k];
  }
  
  // zero the imaginary portion of the N_d/2 coefficient
  //
  output_a[N_d+1] = 0.0;
  
  // exit gracefully
  //
  return ISIP_TRUE;
}