ft_sr_1.cc

这是处理语音信号的程序

// file: $PDSP/class/fourier_transform/v3.0/ft_sr_1.cc
//
 
// isip include files
//
#include "fourier_transform.h"
#include "fourier_transform_constants.h"

// method: sr_complex_cc
//
// this code is based on the code by G.A.Sitton of the Rice University.
// the theory behind the implementation closely follows the description of the
// split-radix algorithm in:
//
// J.G.Proakis, D.G.Manolakis, "Digital Signal Processing - Principles,
// Algorithms, and Applications" 2nd Edition, Macmillan Publishing Company,
// New York, pp. 727-730, 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 2*N
//                      memory is allocated by the calling
//                      program.
//
// return : logical_1 indicating status of computation
//
// this method computes the fourier transform using the split-radix
// algorithm 
//
logical_1 Fourier_transform::sr_complex_cc(float_8* output_a,
						    float_8* input_a) {

  // declare local variables
  //
  int_4	i, j, k, m, q, ii, i0, i1, i2, i3, is, id, temp;
  float_8 cc1, ss1, cc3, ss3, r1, r2, r3, s1, s2, t;
  int_4 n1, n2, n4;
  
  // initialize lookup tables and temporary memory
  //
  if (sr_init_cc(N_d) == ISIP_FALSE) {
    error_handler_cc((char_1*)"sr_real_cc",
		     (char_1*)"error generating the lookup table");
    return ISIP_FALSE;
  }

  // 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*)"sr_real_cc",
		     (char_1*)"order must be a power of 2");
    return ISIP_FALSE;
  }
  
  // copy input real data to temporary memory and input imaginary data
  // into output so that the input data is retained
  // after calculations.
  //
  for (k = 0; k < N_d; ++k) {
    sr_temp_d[k] = input_a[k*2];
    output_a[k] = input_a[k*2+1];
  }
  
  // bit reversal routine
  //
  n1 = N_d - 1;
  n2 = N_d >> 1;
  for (i = j = 0; i < n1; ++i, j += k) {
    if (i < j) {
      t = sr_temp_d[j];
      sr_temp_d[j] = sr_temp_d[i];
      sr_temp_d[i] = t;
      t = output_a[j];
      output_a[j] = output_a[i];
      output_a[i] = t;
    }
    k = n2;
    while (k <= j) {
      j -= k;
      k = k >> 1;
    } 
  }
  
  // Length two butterflies	
  //
  is = 0;
  id = 1;
  n1 = N_d - 1;
  
  do {
    id = id << 2;      
    for (i0 = is; i0 < N_d; i0 += id) {
      i1 = i0 + 1;
      
      r1 = sr_temp_d[i0];
      sr_temp_d[i0] = r1 + sr_temp_d[i1];
      sr_temp_d[i1] = r1 - sr_temp_d[i1];
	
      r1 = output_a[i0];
      output_a[i0] = r1 + output_a[i1];
      output_a[i1] = r1 - output_a[i1]; 
      }
    temp = id - 1;
    is = temp << 1;
  } while (is < n1);
  
  // L shaped butterflies	
  //
  n2 = 2;
  for (k = 2; k <= m; ++k) {
    temp = n2;
    n4 = temp >> 1;
    n2 = n2 << 1;
    is = 0;
    temp = n2;
    id = temp << 1;

    // Unity twiddle factor loops	
    //
    do {
      for (i0 = is; i0 < n1; i0 += id) {
	i1 = i0 + n4;
	i2 = i1 + n4;
	i3 = i2 + n4;	
	r3 = sr_temp_d[i2] + sr_temp_d[i3];
	r2 = sr_temp_d[i2] - sr_temp_d[i3];
	r1 = output_a[i2] + output_a[i3];
	s2 = output_a[i2] - output_a[i3];
	
	sr_temp_d[i2] = sr_temp_d[i0] - r3;
	sr_temp_d[i0] = sr_temp_d[i0] + r3;
	sr_temp_d[i3] = sr_temp_d[i1] - s2;
	sr_temp_d[i1] = s2 + sr_temp_d[i1];
	
	output_a[i2] = output_a[i0] - r1;
	output_a[i0] = output_a[i0] + r1;
	output_a[i3] = output_a[i1] + r2;
	output_a[i1] = output_a[i1] - r2;
	
      } 
      is = (id << 1) - n2;
      id = id << 2; 
      
    } while (is < n1);
    
    q = ii = (int_4)(N_d / (float_8)n2);
    
    // Non-trivial twiddle factor loops	
    //
    for (j = 1; j < n4; ++j, q += ii) {
      cc1 = sr_wr_d[q];
      ss1 = sr_wi_d[q];
      i3 = q * 3;
      cc3 = sr_wr_d[i3];
      ss3 = sr_wi_d[i3];
      
      is = j;
      id = n2 << 1;
      
      do {
	for (i0 = is; i0 < n1; i0 += id) {
	  i1 = i0 + n4;
	  i2 = i1 + n4;
	  i3 = i2 + n4;
	     
	  r1 = (sr_temp_d[i2] * cc1) + (output_a[i2] * ss1);
	  s1 = (output_a[i2] * cc1) - (sr_temp_d[i2] * ss1);
	  r2 = (sr_temp_d[i3] * cc3) + (output_a[i3] * ss3);
	  s2 = (output_a[i3] * cc3) - (sr_temp_d[i3] * ss3);
	  
	  r3 = r1 + r2;
	  r2 = r1 - r2;
	  r1 = s1 + s2;
	  s2 = s1 - s2;
	  
	  sr_temp_d[i2] = sr_temp_d[i0] - r3;
	  sr_temp_d[i0] = r3 + sr_temp_d[i0];
	  sr_temp_d[i3] = sr_temp_d[i1] - s2;
	  sr_temp_d[i1] = s2 + sr_temp_d[i1];
	  
	  output_a[i2] = output_a[i0] - r1;
	  output_a[i0] = r1 + output_a[i0];
	  output_a[i3] = output_a[i1] + r2;
	  output_a[i1] = output_a[i1] - r2;
	}
	is = (id << 1) - n2 + j;
	id = id << 2;
      } while (is < n1);
    }
  }
  
  // interlace the output and store in the output_a array
  //
  for (k = N_d - 1; k >= 0; --k) {
    output_a[k*2+1] = output_a[k];
    output_a[k*2] = sr_temp_d[k];    
  }

  // exit gracefully
  //
  return ISIP_TRUE;
}