ft_rad2_0.cc

这是处理语音信号的程序

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

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

// method: rad2_real_cc
//
// this method performs a radix-2 cooley-tukey decimation-in-frequency
// algorithm.  the code follows the exact structure of that in the book:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 731-732, 1992.
//
// 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 discrete fourier transform
//
logical_1 Fourier_transform::rad2_real_cc(float_8* output_a,
						   float_8* input_a){

  // declare local variables
  //
  int_4 i,j,k,l,m,i1,i2;
  int_4 n1,n2;
  float_8 real_twid,imag_twid;
  float_8 tempr,tempi; 
  
  // initialize the lookup tables
  //
  if (rad2_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"rad2_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] = input_a[k];
    rad2_temp_d[k] = 0;
  }
  
  // check the order
  //
  if (is_power_cc((int_4&)m, (int_4)2) == ISIP_FALSE) {
    error_handler_cc((char_1*)"rad2_real_cc",
		     (char_1*)"order must be a power of two");
    return ISIP_FALSE;
  }

  // actual dif radix-2 computation
  //
  n2 = N_d;

  // looping  through the first  stage
  //
  n1 = n2;
  n2 = n2 >> 1;
  i1 = 0;
  i2 = (int_4)((N_d)/n1);
  
  for (j = 0; j < n2; ++j) {
    
    // loop through each butterfly
    //
    real_twid = rad2_wr_d[i1];
    imag_twid = rad2_wi_d[i1];
    i1 = i1 + i2;
    
    for (k = j; k < N_d; k += n1) {
	
      // perform the butterfly computation
      //
      l = k + n2;
      tempr = output_a[k] - output_a[l];
      output_a[k] = output_a[k] + output_a[l];
      output_a[l] = (real_twid * tempr);
      rad2_temp_d[l] = -(imag_twid * tempr);
    }
  }
  
  for (i = 1; i < m; ++i) {
    
    // looping  through each of the remaining stages
    //
    n1 = n2;
    n2 = n2 >> 1;
    i1 = 0;
    i2 = (int_4)((N_d)/n1);
    
    for (j = 0; j < n2; ++j) {

      // loop through each butterfly
      //
      real_twid = rad2_wr_d[i1];
      imag_twid = rad2_wi_d[i1];
      i1 = i1 + i2;
      
      for (k = j; k < N_d; k += n1) {

	// perform the butterfly computation; computation within square braces
	// in eqn (9.3.34) multiplied by the twiddle factor
	//
	l = k + n2;
        tempr = output_a[k] - output_a[l];
        output_a[k] = output_a[k] + output_a[l];
        tempi = rad2_temp_d[k] - rad2_temp_d[l];
        rad2_temp_d[k] = rad2_temp_d[k] + rad2_temp_d[l];
        output_a[l] = (real_twid * tempr) + (imag_twid * tempi);
        rad2_temp_d[l] =  (real_twid * tempi) - (imag_twid * tempr);
      }
    }
  }
  
  // procedure for bit reversal
  //
  j = 0;
  n1 = N_d - 1;
  n2 = N_d >> 1;
  
  for (i = 0; i < n1 ; ++i) {
    if (i < j) {
      tempr = output_a[j];
      output_a[j] = output_a[i];
      output_a[i] = tempr;
      tempr = rad2_temp_d[j];
      rad2_temp_d[j] = rad2_temp_d[i];
      rad2_temp_d[i] = tempr;
    }
    k = n2;
    while (k <= j) {
      j -= k;
      k = k >> 1;
    }
    j += k;
  }
  
  // interlace the output and store in the output array
  //
  for (k = N_d - 1; k >= 0; --k) {
    output_a[k*2+1] = rad2_temp_d[k];
    output_a[k*2] = output_a[k];    
  }
 
  // exit gracefully
  //
  return ISIP_TRUE;
}