ft_fh_0.cc

这是处理语音信号的程序

// file: $PDSP/class/fourier_transform/v3.0/ft_fh_0.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 mafyer
//
#define FT_FHT_GOOD_TRIG
#define ISIP_SQRT2 (float_8)1.41421356237

// 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 N
//			input data memory should be allocated by the
//			calling program.
//
// return: a logical_1 value indicating status
//
// implements a real fast hartley transform
//
logical_1 Fourier_transform::fh_real_cc(float_8* output_a,
						 float_8* input_a){
  
  // declare local variables
  //
  float_8 a, b;
  float_8 c1, s1, s2, c2, s3, c3, s4, c4;
  float_8 f0, g0, f1, g1, f2, g2, f3, g3;
  int_4 i, j, kk, k, k1, k2, k3, k4, kx, m, temp;
  float_8* fi, *fn, *gi;
  int_4 N_d_2 = N_d >> (int_4)1;

  // initializes the variables in the lookup table
  //
  TRIG_VARS;
  USE_STATIC;
  
  // initialize the lookup tables
  //
  if (fh_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"fh_real_cc",
		     (char_1*)"error initializing the algorithm");
    return ISIP_FALSE;
  }
  
  // copy input data to temporary memory so that the input data is retained
  // after calculations.
  //
  for (k = 0; k < N_d; ++k) {
    output_a[k] = (float_8)0;
    fh_temp_d[k] = input_a[k];
  }
  
  // check the order
  //
  if (is_power_cc((int_4&)m, (int_4)2) == ISIP_FALSE) {
    error_handler_cc((char_1*)"fh_real_cc",
		     (char_1*)"order must be a power of two");
    return ISIP_FALSE;
  }

  // perform bit reversal
  //
  for (k1 = 1, k2 = 0; k1 < N_d; ++k1) {
    for (k = N_d >> 1; (!((k2 ^= k) & k)); k >>= 1) {};
    if (k1 > k2) {
      a = fh_temp_d[k1];
      fh_temp_d[k1] = fh_temp_d[k2];
      fh_temp_d[k2] = a;
    }
  }

  // get the exponent of two
  //
  for (k = 0; (1 << k) < N_d; ++k);
  k  &= (int_4)1;

  // twiddleless radix-2 computations when order is an even power of 2
  // completes the stage before the g's in fig.2
  //
  if (k == 0) {
    for (fi = fh_temp_d, fn = fh_temp_d + N_d ;fi < fn; fi += 4) {
      f1 = fi[0] - fi[1];
      f0 = fi[0] + fi[1];
      f3 = fi[2] - fi[3];
      f2 = fi[2] + fi[3];
      fi[2] = f0 - f2;	
      fi[0] = f0 + f2;
      fi[3] = f1 - f3;	
      fi[1] = f1 + f3;
    }
  }

  // twiddleless radix-2 computations when order is an odd power of 2
  // completes the stages till g's in fig. 1
  //
  else {
    
    for (fi = fh_temp_d, fn = fh_temp_d + N_d, gi = fi + 1;
	 fi < fn; fi += 8, gi += 8) {
      c1 = fi[0] - gi[0];
      s1 = fi[0] + gi[0];
      c2 = fi[2] - gi[2];
      s2 = fi[2] + gi[2];
      c3 = fi[4] - gi[4];
      s3 = fi[4] + gi[4];
      c4 = fi[6] - gi[6];
      s4 = fi[6] + gi[6];
      f1 = s1 - s2;	
      f0 = s1 + s2;
      g1 = c1 - c2;	
      g0 = c1 + c2;
      f3 = s3 - s4;	
      f2 = s3 + s4;
      g3 = c4 * ISIP_SQRT2;		
      g2 = c3 * ISIP_SQRT2;
      fi[4] = f0 - f2;
      fi[0] = f0 + f2;
      fi[6] = f1 - f3;
      fi[2] = f1 + f3;
      gi[4] = g0 - g2;
      gi[0] = g0 + g2;
      gi[6] = g1 - g3;
      gi[2] = g1 + g3;
    }
  }

  // stop computations here if order is less than 16
  //
  if (N_d < 16) { 
        
    if (data_type_d == FT_REAL) {
      
      //  interlace the output and put into output_a
      //
      output_a[0] = fh_temp_d[0];
      output_a[1] = 0.0;

      for ( kk = 1; kk <= N_d_2; ++kk ){
	output_a[2*kk] = output_a[2*N_d - 2*kk ] = fh_temp_d[kk];
	output_a[2*N_d - 2*kk + 1 ] = fh_temp_d[N_d - kk];
	output_a[2*kk+1] = -fh_temp_d[N_d - kk];
      }
      output_a[N_d+1] = 0.0;
    }
    
    if (data_type_d == FT_COMPLEX) {

      // get the output data
      //
      for (kk = 0; kk < N_d; ++kk) {
	output_a[kk] = fh_temp_d[kk];
      }
    }
    
    // exit gracefully
    //
    return ISIP_TRUE;
  }

  // if order is greater than 16, continue twiddled computations
  //

  // does the computations for the even indices
  //
  do {
    k  += (int_4)2;
    k1  = 1 << k;
    temp = k1;
    k2  = temp << (int_4)1;
    temp = k2;
    k4  = temp << (int_4)1;
    k3  = (int_4)k1 + (int_4)k2;
    temp = k1;
    kx  = temp >> 1;
    fi  = fh_temp_d;
    gi  = (float_8*)((float_8*)fi + (int_4)kx);
    fn  = (float_8*)((float_8*)fh_temp_d + N_d);
    do {
      f1 = fi[0 ] - fi[k1];
      f0 = fi[0 ] + fi[k1];
      f3 = fi[k2] - fi[k3];
      f2 = fi[k2] + fi[k3];
      fi[k2]  = f0 - f2;
      fi[0 ]  = f0 + f2;
      fi[k3]  = f1 - f3;
      fi[k1]  = f1 + f3;
      g1      = gi[0 ] - gi[k1];
      g0      = gi[0 ] + gi[k1];
      g3      = gi[k3] * ISIP_SQRT2;
      g2      = gi[k2] * ISIP_SQRT2;
      gi[k2]  = g0 - g2;
      gi[0 ]  = g0 + g2;
      gi[k3]  = g1 - g3;			
      gi[k1]  = g1 + g3;
      gi = (float_8*)((float_8*)gi + (int_4)k4);
      fi = (float_8*)((float_8*)fi + (int_4)k4);
    } while ((float_8*)fi < (float_8*)fn);
    
    TRIG_INIT(k, c1, s1);

    // does computations for odd indices, this includes multiplication by
    // the twiddles. note that these computations use some sort look ahead
    // and avoids redundancy of computations.
    //
    for (i = 1; (int_4)i < (int_4)kx; ++i) {

      // uses the euler form of approximation to compute the cosines
      //
      TRIG_NEXT(k, c1, s1);
      c2 = (c1 * c1) - (s1 * s1);
      s2 = (c1 * s1) * 2;
      fn = (float_8*)((float_8*)fh_temp_d + N_d);
      fi = (float_8*)((float_8*)fh_temp_d + (int_4)i);
      gi = (float_8*)((float_8*)fh_temp_d +  (int_4)k1 - (int_4)i);
      do {
	b = s2*fi[k1] - c2*gi[k1];
	a = c2*fi[k1] + s2*gi[k1];
	f1 = fi[0 ] - a;
	f0 = fi[0 ] + a;
	g1 = gi[0 ] - b;
	g0 = gi[0 ] + b;
	b = s2*fi[k3] - c2*gi[k3];
	a = c2*fi[k3] + s2*gi[k3];
	f3 = fi[k2] - a;
	f2 = fi[k2] + a;
	g3 = gi[k2] - b;
	g2 = gi[k2] + b;
	b = s1*f2  - c1*g3;
	a = c1*f2  + s1*g3;
	fi[k2]  = f0 - a;
	fi[0 ]  = f0 + a;
	gi[k3]  = g1 - b;
	gi[k1]  = g1 + b;
	b = c1*g2 - s1*f3;
	a = s1*g2 + c1*f3;
	gi[k2]  = g0 - a;
	gi[0 ]  = g0 + a;
	fi[k3]  = f1 - b;
	fi[k1]  = f1 + b;
	gi = (float_8*)((float_8*)gi + (int_4)k4);
	fi = (float_8*)((float_8*)fi + (int_4)k4);
      } while ((float_8*)fi < (float_8*)fn);
    }
    TRIG_RESET(k, c1, s1);
  } while ((int_4)k4 < N_d);
  
  a = b = (float_8)0;  

  if (data_type_d == FT_COMPLEX) {
    
    // get the output data
    //
    for (kk = 0; kk < N_d; ++kk) {
      output_a[kk] = fh_temp_d[kk];
    }
  }
  
  if (data_type_d == FT_REAL) {
    
    // get fourier coefficients from hartley coefficients
    //
    for (i = 1, j = (N_d - 1), k = N_d_2; (int_4)i < (int_4)k; ++i, --j) {
      a = fh_temp_d[i];
      b = fh_temp_d[j];
      fh_temp_d[j] = (a - b) * (float_8)0.5;
      fh_temp_d[i] = (a + b) * (float_8)0.5;
    }
    
    //  interlace the output and put into output_a
    //
    output_a[0] = fh_temp_d[0];
    output_a[1] = (float_8)0;

    for (kk = 1; kk <= N_d_2; ++kk) {
      output_a[kk*2] = output_a[2*N_d - kk*2 ] = fh_temp_d[kk];
      output_a[kk*2+1] = -fh_temp_d[N_d - kk];
      output_a[2*N_d - kk*2 + 1 ] = fh_temp_d[N_d - kk];
    }
    output_a[N_d+1] = 0.0;
  }
  
  // exit gracefully
  //
  return ISIP_TRUE;
}