ft_sr_0.cc

这是处理语音信号的程序

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

// method: sr_real_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 N
//                      memory is allocated by the calling
//                      program.
//
// return: a logical_1 value indicating status
//
// routine used to compute the real split-radix fft
//
logical_1 Fourier_transform::sr_real_cc(float_8* output_a,
						 float_8* input_a) {

  // declare local variables
  //
  float_8   t1, t2, t3, t4, t5, t6, t;
  float_8   cc1, ss1, cc3, ss3;
  int_4     i, j, k, l, m, n_d_2, inc;
  int_4     n1, n2, n4, n8;
  int_4     is, id, i1, i2, i3, i4, i5, i6, i7, i8;
  
  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;
  }
  
  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 data to temporary memory so that the input data is retained
  // after calculations.
  //
  for (k = 0; k < N_d; ++k) {
    sr_temp_d[k] = input_a[k];
    output_a[k] = 0;
  }

  // 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;
    }

    k = n2;
    while (k <= j) {
      j -= k;
      k = k >> 1;
    } 
  }
  
  //  First (twiddleless) L Butterflies; completes stage A in fig. 9.23   
  //
  is = 0;
  id = 4;
  n1 = N_d - 1;
  
  do {
    for (i = is; i < n1; i += id) {
      t1 = sr_temp_d[i];
      i1 = i + 1;
      t2 = sr_temp_d[i1];
      sr_temp_d[i] = t1 + t2;
      sr_temp_d[i1] = t1 - t2;
    }
    is = (id << 1) - 2;
    id = id << 2;
  } while (is < n1); 

  //  Other (twiddled) L Butterflies  
  //
  n2 = (int_4)2;
  for (k = 2; k <= m; ++k) {
    n4 = n2 >> 1;
    n8 = n4 >> 1;
    n2 = n2 << 1;
    is = 0;
    id = n2 << 1; 
    do {
      for (i = is; i < N_d; i += id) {
        i2 = i + n4;
        i3 = i2 + n4;
        i4 = i3 + n4;
        t3 = sr_temp_d[i3];
        t4 = sr_temp_d[i4];
        t1 = sr_temp_d[i];
        t2 = t4 + t3;
        sr_temp_d[i4] = t4 - t3;
        sr_temp_d[i3] = t1 - t2;
        sr_temp_d[i] = t1 + t2;
 
        if (n4 != 1) {
          i1 = i + n8;
          i2 = n8 + i2;
          i3 = n8 + i3;
          i4 = n8 + i4;
          t3 = sr_temp_d[i3];
          t4 = sr_temp_d[i4];
          t1 = (t3 + t4) * ISIP_INV_SQRT2;
          t2 = (t3 - t4) * ISIP_INV_SQRT2;
          t3 = sr_temp_d[i1];
          t4 = sr_temp_d[i2];
          sr_temp_d[i4] = t4 - t1;
          sr_temp_d[i3] = -(t4 + t1);
          sr_temp_d[i2] = t3 - t2;
          sr_temp_d[i1] = t3 + t2;
        }
      }
      is = (id << 1) - n2;
      id = id << 2;
    } while (is < n1);
    
    inc = N_d/n2;
    l = inc;
    for (j = 1; j < n8; ++j) {
      cc1 = sr_wr_d[l];
      ss1 = sr_wi_d[l];

      // i1 gives the twiddle factors which are multiples of 3 and are needed
      // at the B stage in fig. 9.23
      i1 = l * 3;
      cc3 = sr_wr_d[i1];
      ss3 = sr_wi_d[i1];
      l = (j + 1) * inc;
      is = 0;
      id = n2 << 1;
      do {
	for (i = is; i < N_d; i += id) {
	  i1 = i + j;
	  i2 = i1 +  n4;
	  i3 = i2 +  n4;
	  i4 = i3 +  n4;
	  i5 = i - j + n4;
	  i6 = i5 +  n4;
	  i7 = i6 +  n4;
	  i8 = i7 +  n4;
	  t3 = sr_temp_d[i3];
	  t4 = sr_temp_d[i7];
	  t1 = (t3*cc1) + (t4*ss1);
	  t2 = (t4*cc1) - (t3*ss1);
	  t5 = sr_temp_d[i4];
	  t6 = sr_temp_d[i8];
	  t3 = (t5*cc3) + (t6*ss3);
	  t4 = (t6*cc3) - (t5*ss3);
	  t5 = t1 + t3;
	  t3 = t1 - t3;
	  t6 = t2 + t4;
	  t4 = t2 - t4;
	  t1 = sr_temp_d[i6];
	  sr_temp_d[i8] = t1 + t6;
	  sr_temp_d[i3] = t6 - t1;
	  t2 = sr_temp_d[i2];
	  sr_temp_d[i4] = t2 - t3;
	  sr_temp_d[i7] = -(t2 + t3);
	  t1 = sr_temp_d[i1];
	  sr_temp_d[i1] = t1 + t5;
	  sr_temp_d[i6] = t1 - t5;
	  t5 = sr_temp_d[i5];
	  sr_temp_d[i2] = t5 + t4;
	  sr_temp_d[i5] = t5 - t4;
	}
	is = (id << 1) - n2;
	id = id << 2;
      } while (is < n1);
    }
  }
    
  // interlace the output and store in the output_a array
  //
  n_d_2 = N_d/2;
  for (k = 1; k <= n_d_2; ++k) {
    output_a[k*2] = output_a[2*N_d - k*2] = sr_temp_d[k];
    output_a[2*N_d - k*2 + 1 ] = -sr_temp_d[N_d - k];
    output_a[k*2+1] = sr_temp_d[N_d - k];
  }
  
  // manually set values to special points in the output data
  // zero imaginary part of dc, and N/2 point
  //
  output_a[0] = sr_temp_d[0];
  output_a[1] = (float_8)0.0;
  output_a[N_d+1] = (float_8)0.0;
  
  // exit gracefully
  //
  return ISIP_TRUE;
}