ft_ditf_0.cc

这是处理语音信号的程序

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

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

// method: ditf_real_cc
//
// This method is based on the DITF algorithm by Ali Saidi and implements the
// butterfly version of the algorithm rather than the recursive version
// for reasons of efficiency.
//
//   Ali Saidi, "Decimation-In-Time-Frequency FFT Algorithm", Proceedings of
//   ICASSP 1994, pp. 453-456, SanFrancisco, CA, April 1994.
//
// 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
//
// this method performs the real ditf algorithm
//
logical_1 Fourier_transform::ditf_real_cc(float_8* output_a,
					  float_8* input_a) {
  
  // declare local variables
  //
  int_4 transition_stage = (int_4)0;
  int_4 trans_butterfly_size = (int_4)0;
  int_4 m, m1, p, l, i1, i2, i, ii, jj, j, ib, jb, kb, k, kk, temp;
  int_4 n1, n2, n3, n4, n5, pn2;
  float_8 real_twid, imag_twid, tempr, tempi;
    
  // initialize the lookup tables
  //
  if (ditf_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"ditf_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];
    ditf_temp_d[k] = 0;
  }
  
  // check the order
  //
  if (is_power_cc((int_4&)m, (int_4)2) == ISIP_FALSE) {
    error_handler_cc((char_1*)"ditf_real_cc",
                     (char_1*)"order must be a power of two");
    return ISIP_FALSE;
  }
  
  // initialize indices for the butterflies
  //
  n1 = N_d;   
  n2 = N_d >> (int_4)1;
  
  // only do the DIT portion for orders greater than 4
  //
  if ( (int_4)4 < N_d ) {
    
    // calculate the location of the transition stage
    //
    temp = m;
    transition_stage = temp >> (int_4)1;
    
    // do the first decimation-in-time twiddleless butterfly
    //
    for (j = 0; j < n2; ++j) {
      l = j + n2;
      tempr = output_a[j];
      tempi = ditf_temp_d[j];
      output_a[j] = tempr + output_a[l];
      ditf_temp_d[j] = tempi + ditf_temp_d[l];
      output_a[l] = tempr - output_a[l];
      ditf_temp_d[l] = tempi - ditf_temp_d[l];
    }
    
    // initialize a counter for the internal dit bit-reversals
    //
    n3 = 1;
	
    // do the remaining dit butterflies
    //
    for (i = 1; i < transition_stage; ++i) {
      n1 = n2;
      n2 = n2 >> (int_4)1;
      
      // procedure for bit reversal
      //
      jb = 0;
      n3 = n3 << (int_4)1;
      temp = n3;
      n5 = temp >> (int_4)1;
      n4 = n3 - 1;
      
      for (ib = 0; ib < n3; ++ib) {
	ditf_indices_d[ib] = ib;
      }
      for (ib = 0; ib < n4 ; ++ib) {
	
	if (ib < jb) {
	  ditf_indices_d[jb] = ib;
	  ditf_indices_d[ib] = jb;
	}
	kb = n5;
	while (kb <= jb) {
	  jb = jb - kb;
	  kb = kb >> (int_4)1;
	}
	jb = jb + kb;
      }
	    
      // carry out the dit butterfly
      //
      p = 0;
      for (j = 0; j < n3 ; ++j) {
	if ( j != 0 ) {
	  real_twid = ditf_wr_d[n2*ditf_indices_d[j]];
	  imag_twid = ditf_wi_d[n2*ditf_indices_d[j]];
	  
	  for (k = p; k < p+n2; ++k) {
	    l = k + n2;
	    tempr = (real_twid * output_a[l]) +
	      (imag_twid * ditf_temp_d[l]);
	    tempi = (real_twid * ditf_temp_d[l]) -
	      (imag_twid * output_a[l]);
	    output_a[l] = output_a[k] - tempr;
	    ditf_temp_d[l] = ditf_temp_d[k] - tempi;
	    output_a[k] = output_a[k] + tempr;
	    ditf_temp_d[k] = ditf_temp_d[k] + tempi;
	  }
	}

	else {
	  pn2 = p + n2;
	  for (k = p; k < pn2; ++k) {
	    l = k + n2;
	    tempr = output_a[l];
	    tempi = ditf_temp_d[l];
	    output_a[l] = output_a[k] - tempr;
	    ditf_temp_d[l] = ditf_temp_d[k] - tempi;
	    output_a[k] = output_a[k] + tempr;
	    ditf_temp_d[k] = ditf_temp_d[k] + tempi;
	  }
	}
	p = p + n1;
      }
    }
    
    // find the transition stage factors
    //
	
    // procedure for bit reversal
    //
    n3 = n3 << 1;
    int_4 jb = 0;
    n4 = n3 - 1;
    temp = n3;
    n5 = temp >> (int_4)1;
    for (ib = 0; ib < n3; ++ib) {
      ditf_indices_d[ib] = ib;
    }
    for (ib = 0; ib < n4 ; ++ib){
	    
      if (ib < jb) {
	ditf_indices_d[jb] = ib;
	ditf_indices_d[ib] = jb;
      }
      
      kb = n5;
      while (kb <= jb) {
	jb = jb - kb;
	kb = kb >> (int_4)1;
      }
      jb = jb + kb;
    }
    
	
    n5 = 0;
    trans_butterfly_size = N_d >> transition_stage;
    for (jj = 0; jj < n3; ++jj) {
      for (ii = 0; ii < trans_butterfly_size; ++ii) {
	ditf_trans_factor_indices_d[n5] = ii * ditf_indices_d[jj];
	++n5;
      }
    }
	
    // multiply by the transition stage factors.
    //
    for (i = 0; i < N_d; ++i) {
      real_twid = ditf_wr_d[ditf_trans_factor_indices_d[i]];
      imag_twid = ditf_wi_d[ditf_trans_factor_indices_d[i]];
      tempr = output_a[i];
      tempi = ditf_temp_d[i];
      output_a[i] = (real_twid * tempr) + (imag_twid * tempi);
      ditf_temp_d[i] = (real_twid * tempi) - (imag_twid * tempr);
    }
  }
  
  // reset the stage counter
  //
  n2 = N_d >> transition_stage;
  
  // carry out the decimation-in-frequency portion
  //
  m1 = m - 1;
  for (i = transition_stage; i < m1; ++i) {
    n1 = n2;
    n2 = n2 >> (int_4)1;
    i1 = 0;
    i2 = (int_4)((N_d)/n1);
    for (j = 0; j < n2; ++j) {
      if ( j != 0 ) {
	real_twid = ditf_wr_d[i1];
	imag_twid = ditf_wi_d[i1];
	i1 = i1 + i2;
	
	for (k = j; k < N_d; k += n1) {
	  l = k + n2;
	  tempr = output_a[k] - output_a[l];
	  output_a[k] = output_a[k] + output_a[l];
	  tempi = ditf_temp_d[k] - ditf_temp_d[l];
	  ditf_temp_d[k] = ditf_temp_d[k] + ditf_temp_d[l];
	  output_a[l] = (real_twid * tempr) + (imag_twid * tempi);
	  ditf_temp_d[l] =  (real_twid * tempi) - (imag_twid * tempr);
	}
      }
      else {
	i1 = i1 + i2;
	for (k = j; k < N_d; k += n1) {
	  l = k + n2;
	  tempr = output_a[k] - output_a[l];
	  output_a[k] = output_a[k] + output_a[l];
	  tempi = ditf_temp_d[k] - ditf_temp_d[l];
	  ditf_temp_d[k] = ditf_temp_d[k] + ditf_temp_d[l];
	  output_a[l] = tempr;
	  ditf_temp_d[l] = tempi;
	}
      }
    }
  }
    
  // Calculate last twiddleless stage of the DIF portion
  //
  for (j = 0; j < N_d; j += 2) {
    l = j + 1;
    tempr = output_a[j];
    tempi = ditf_temp_d[j];
    output_a[j] = tempr + output_a[l];
    ditf_temp_d[j] = tempi + ditf_temp_d[l];
    output_a[l] = tempr - output_a[l];
    ditf_temp_d[l] = tempi - ditf_temp_d[l];
  }
  
  // procedure for bit reversal
  //
  temp = 0;
  j = 0;
  n1 = N_d - 1;
  n2 = N_d >> (int_4)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 = ditf_temp_d[j];
      ditf_temp_d[j] = ditf_temp_d[i];
      ditf_temp_d[i] = tempr;
    }
    k = n2;
    while (k <= j) {
      j = j - k;
      k = k >> (int_4)1;
    }
    j = j + k;
  }
  
  // interlace the output and store in the output_a array
  //
  for (kk = N_d - 1; kk >= 0; --kk) {
    output_a[2*kk+1] = ditf_temp_d[kk];
    output_a[2*kk] = output_a[kk];    
  }
  
  // exit gracefully
  //
  return ISIP_TRUE;
}