ft_rad4_1.cc

这是处理语音信号的程序

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

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

// method: rad4_complex_cc
//
// This method performs a standard radix-4 decimation-in-frequency algorithm
// The code follows the exact structure of that in the book:
//
//   C.S. Burrus and T.W. Parks, "DFT/FFT and Convolution Algorithms - Theory
//   and Implementation", 2nd Edition, John Wiley and Sons, New York, NY,1985
// 
// 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
//
logical_1 Fourier_transform::rad4_complex_cc(float_8* output_a,
						      float_8* input_a) {

  // declare local variables
  //
  int_4 i, j, k, m, i1, i2, i3;
  int_4 n1, n2, N_d_1;
  float_8 r1, r2, r3, r4, s1, s2, s3, s4, co1, co2, co3, si1, si2, si3;
  int_4 a, b, c;
  
  // check the order
  //
  if (is_power_cc((int_4&)m, (int_4)4) == ISIP_FALSE) {
    error_handler_cc((char_1*)"rad4_complex_cc",
		     (char_1*)"order must be a power of 4");
    return ISIP_FALSE;
  }
  n2 = N_d;

  // create lookup table ( we use the same lookup table generator as
  // we used for rad2 computations )
  //
  if (rad4_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"rad4_init_cc",
		     (char_1*)"error initializing the algorithm");
    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) {
    output_a[k] = input_a[k*2];
    rad4_temp_d[k] = input_a[k*2 + 1];
  }
  
  // main computation loop; looping through all the stages; there are
  // log4(N_d) stages.
  //
  for (k = 0; k < m; ++k) {
    n1 = n2;
    n2  = n2 >> (int_4)2;
    a = 0;

    // go through each butterfly in a stage; there are N_d/4 butterflies
    // in each stage.
    //
    for (j = 0; j < n2; ++j) {
      b = a + a;
      c = a + b;
      co1 = rad4_wr_d[(int_4)a];
      co2 = rad4_wr_d[(int_4)b];
      co3 = rad4_wr_d[(int_4)c];
      si1 = rad4_wi_d[(int_4)a];
      si2 = rad4_wi_d[(int_4)b];
      si3 = rad4_wi_d[(int_4)c];
      a = (j + 1)*(int_4)pow(4,k);
      
      // compute each radix-4 butterfly;
      //
      for (i = j; i < N_d; i += n1) {
	i1 = i + n2;
	i2 = i1 + n2;
	i3 = i2 + n2;
	r1 = output_a[i] + output_a[i2];
	r3 = output_a[i] - output_a[i2];
	s1 = rad4_temp_d[i] + rad4_temp_d[i2];
        s3 = rad4_temp_d[i] - rad4_temp_d[i2];
        r2 = output_a[i1] + output_a[i3];
        r4 = output_a[i1] - output_a[i3];
        s2 = rad4_temp_d[i1] + rad4_temp_d[i3];
        s4 = rad4_temp_d[i1] - rad4_temp_d[i3];
        output_a[i] = r1 + r2;
        r2 = r1 - r2;      
        r1 = r3 - s4;
        r3 = r3 + s4;     
        rad4_temp_d[i] = s1 + s2;
        s2 = s1 - s2;      
        s1 = s3 + r4;
        s3 = s3 - r4;
	output_a[i1] = (r3 * co1) + (s3 * si1);
	rad4_temp_d[i1] = (s3 * co1) - (r3 * si1);
	output_a[i2] = (r2 * co2) + (s2 * si2);
	rad4_temp_d[i2] = (s2 * co2) - (r2 * si2);
	output_a[i3] = (r1 * co3) + (s1 * si3);
	rad4_temp_d[i3] = (s1 * co3) - (r1 * si3);
      }
    }
  }
  
  // procedure for bit reversal
  //
  float_8 temp = 0;
  j = 0;
  N_d_1 = N_d - 1;
  for (i = 0; i < N_d_1; i++) {
    if (i < j) {
      temp = output_a[j];
      output_a[j] = output_a[i];
      output_a[i] = temp;
      temp = rad4_temp_d[j];
      rad4_temp_d[j] = rad4_temp_d[i];
      rad4_temp_d[i] = temp;
    }
    
    k = N_d >> (int_4)2;
    while (k*3 <= j){
      j -= k*3;
      k  = k >> (int_4)2;
    }
    j += k;
  }

  N_d_1 = N_d - 1;
  for (k = N_d_1; k >= 0; --k) {
    output_a[k*2+1] = rad4_temp_d[k];
    output_a[k*2] = output_a[k];    
  }

  // exit gracefully
  //
  return ISIP_TRUE;
}