ft_dct_0.cc

这是处理语音信号的程序

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

// system include files
//
#include <math.h>

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

#include "fourier_transform_constants.h"

// method: qf_real_dct_cc
//
// This code very closely follows the algorithm described in:
//    H.Guo et.al,"The Quick Discrete Fourier Transform", IEEE Transactions
//    on Acoustics Speech and Signal Processing 1994, pp. 445-448, 1994
//
// This implementation is  adapted from the code written by G.A.Sitton at
// Rice University.
//
// arguments:
//  float_8* coeff_a : (output) output dct coefficient array
//  float_8* data_a  : (input) input data array
//
// return: a logical_1 value indicating status
//
// implements a disrete cosine transform
//
logical_1 Fourier_transform::qf_real_dct_cc(float_8* coeff_a,
						     float_8* data_a){

  // declare local variables
  //
  float_8 t0, t1, t2, t3, t4;
  float_8* XX;
  float_8* xx;
  int_4	i, j, k, n, nq, n2, n2p1, temp;
  
  // Increment level & start	
  //
  n = nc[++ii];
  xx = ws[ii];

  // Last recursion level ?	
  //
  if (n == N2) {
    
    // Form even function
    //
    xx[0] = (float_8)(data_a[0] + data_a[8]);
    xx[1] = (float_8)(data_a[1] + data_a[7]);
    xx[2] = (float_8)(data_a[2] + data_a[6]);
    xx[3] = (float_8)(data_a[3] + data_a[5]);
  
    // Save even DCT coefficients
    //
    t0 = xx[0] - data_a[4];
    t3 = xx[0] + data_a[4];
    t2 = xx[1] + xx[3];
    t1 = (xx[1] - xx[3]) * sqrt2;
    t4 = t3 + xx[2];
    
    coeff_a[0] = t4 + t2;
    coeff_a[2] = t0 + t1;
    coeff_a[4] = t3 - xx[2];
    coeff_a[6] = t0 - t1;
    coeff_a[8] = t4 - t2;
    
    // Form twiddled odd function
    //
    xx[0] = (float_8)(data_a[0] - data_a[8]) * sec0;
    xx[1] = (float_8)(data_a[1] - data_a[7]) * sec1;
    xx[2] = (float_8)(data_a[2] - data_a[6]) * sec2;
    xx[3] = (float_8)(data_a[3] - data_a[5]) * sec3;
    
    // Save odd DCT coefficients
    //
    t1 = (xx[1] - xx[3]) * sqrt2;
    t2 = xx[1] + xx[3];
    t4 = xx[0] + xx[2];
    
    coeff_a[1] = t4 + (t2 + (t0 = (xx[0] + t1)));
    coeff_a[3] = t0 + (t3 = (xx[0] - xx[2]));
    coeff_a[5] = t3 + (t1 = (xx[0] - t1));
    coeff_a[7] = (t1 + t4) - t2; 
  }
  
  else {
    temp = n;
    n2 = temp >> (int_4)1;
    n2p1 = n2 + 1;
    nq = ic[ii];  
      
    if ((k = mc[ii]) > n2)
      k = n2;
    
    // Form even-odd functions
    //
    if (nn <= n2p1)
      for (i = j = 0; i < nn; ++i, j += nq)
	xx[i + n2p1] = (xx[i] = data_a[i]) * sc[j];
    else {
      temp = n - nn;
      for (i = j = 0; i <= temp; ++i, j += nq)
	xx[i + n2p1] = (xx[i] = (float_8)data_a[i]) * sc[j];
      
      for (; i < n2; ++i, j += nq) {
	xx[i] = (float_8)(data_a[i] + data_a[n - i]);
	xx[i + n2p1] = (float_8)(data_a[i] - data_a[n - i]) * sc[j];
      }
      
      xx[i] = data_a[i];
      xx[i + n2p1] = 0.0;
    }
    
    // Do recursive DCTs
    //
    qf_real_dct_cc((float_8*)(XX = &xx[n + 2]), (float_8*)xx);
    qf_real_dct_cc((float_8*)&XX[n2p1], (float_8*)&xx[n2p1]);
    
    // Save DCT coefficients
    //
    t0 = XX[n2p1];
    
    for (i = 0; i < k; i += 2) {
      coeff_a[j = i + i] = XX[i];
      coeff_a[j + 1] = t0 + (t1 = XX[(i + n2p1) + 1]);
      coeff_a[j + 2] = XX[i + 1];
      coeff_a[j + 3] = t1 + (t0 = XX[(i + n2p1) + 2]);
    }
    coeff_a[i + i] = XX[i];
  }
  
  // Decrement level & exit
  //
  --ii;
  
  // exit gracefully
  //
  return ISIP_TRUE;

}