extf_func_4.cc

这是处理语音信号的程序

// file: extract_feature/extf_func_4.cc
//
// isip include files
//
#include "extract_feature.h"
#include "extract_feature_constants.h"

#include <Signal.h>
#include <signal_constants.h>

#include <fourier_transform.h>
#include <fourier_transform_constants.h>

//method: func_aud_weight_cc
//
// arguments:
//  float_8* cb: (output) array of weighting coefficients to simulate 
//                      critical band spectral resolution
//  float_8* eql:(output) array of weighting coefficients to simulate 
//                      equal loudness sensitivity of hearing
//  int_4* ibegen:(output) three-dimensional array which indicates where
//                        to begin and where to end integration of
//			  speech spectrum and where the given
//                        weighting function starts in array cb 
//
//  int_4* npoint: (input) number of points in the fft spectrum
//  int_4*  nfilt_a: (output) number of samples of auditory spectrum equally
//

// This method computes the auditory weighting functions
//
logical_1 func_aud_weight_cc(int_4 npoint_a,int_4& nfilt_a,float_8* cb_a,
		      float_8* eql_a, int_4* ibegen_a) {
  
  // Local variables 
  //
  static float_8 freq, zdel, rsss;
  static int_4 i, j;
  static float_8 x, z, f0, z0, f2samp, fh, fl, fnqbar;
  static int_4 icount;
  static float_8 fsq;
  static int_4 sf=8000;
  static float_8 c_b28 = 10.;

  // Parameter adjustments --- no int_4er needed
  //
  //  --(cb_a - 1);
  //  ibegen_a -= 24;
  
  float_8 r_1 = sf / (float_8)1200.;
  fnqbar = log((float_8)sf / (float_8)1200. + sqrt(r_1 * r_1 + (float_8)1.)) 
    * (float_8) 6.;
  
  // compute number of filters for less than 1 Bark spacing 
  //

  nfilt_a = (int_4 ) fnqbar + 2;
  
  // frequency -> fft spectral sample conversion 
  //

  f2samp = (float_8) (npoint_a - 1) / (sf / (float_8)2.);
  
  // compute filter step in Bark 
  //
  zdel = fnqbar / (float_8) (nfilt_a - 1);
  
  
  // loop over all weighting functions 
  //

  icount = 1;
  int_4 i_1 = nfilt_a - 1;
  for (j = 2; j <= i_1; ++j) {
    (ibegen_a - 24)[j + 69] = icount;

    // get center frequency of the j-th filter in Bark 
    //

    z0 = zdel * (float_8) (j - 1);

    //get center frequency of the j-th filter in Hz 
    //

    f0 = (exp((float_8)z0 / 6.) - exp(-(float_8)z0 / 6)) 
      * (float_8)600. / (float_8) 2.;

    // get low-cut frequency of j-th filter in Hz 
    //

    fl = (exp((z0 - (float_8)2.5) / 6) - exp(-(float_8)(z0 - (float_8)2.5) 
					    / 6)) * (float_8)600. / (float_8)2.;
    r_1 = fl * f2samp;
    (ibegen_a - 24)[j + 23] = (int_4)(((float_8)r_1)+0.5) + 1;
    if ((ibegen_a - 24)[j + 23] < 1) {
      (ibegen_a - 24)[j + 23] = 1;
    }

    // get high-cut frequency of j-th filter in Hz 
    //

    fh = (exp((z0 + (float_8)1.3) / 6) - exp(-(float_8)(z0 + (float_8)1.3) 
					    / 6)) * (float_8)600. / (float_8)2.;
    r_1 = fh * f2samp;
    (ibegen_a - 24)[j + 46] = (int_4)(((float_8)r_1)+0.5) + 1;
    if ((ibegen_a - 24)[j + 46] > npoint_a) {
      (ibegen_a - 24)[j + 46] = npoint_a;
    }

    // do-loop over the power spectrum 
    //

    int_4 i_2 = (ibegen_a - 24)[j + 46];
    for (i = (ibegen_a - 24)[j + 23]; i <= i_2; ++i) {

      //     get frequency of j-th spectral point in Hz 
      //

      freq = (float_8) (i - 1) / f2samp;
      
      // get frequency of j-th spectral point in Bark 
      //

      x = freq / (float_8)600.;
      r_1 = x;
      z = log(x + sqrt(r_1 * r_1 + (float_8)1.)) * (float_8)6.;
      
      // normalize by center frequency in barks: 
      //

      z -= z0;
      
      // compute weighting 
      //
      
      if (z <= (float_8)-.5) {
	float_8 d_1 = (float_8) (z + (float_8).5);
	(cb_a - 1)[icount] = pow(c_b28, d_1);
      } else if (z >= (float_8).5) {
	float_8 d_1 = (float_8) ((z - (float_8).5) * (float_8)-2.5);
	(cb_a - 1)[icount] = pow(c_b28, d_1);
      } else {
	(cb_a - 1)[icount] = (float_8)1.;
      }

      //     calculate the LOG 40 db equal-loudness curve 
      //     at center frequency 
      //
      r_1 = f0;
      fsq = r_1 * r_1;
      r_1 = fsq;
      float_8 r_2 = fsq + (float_8)1.6e5;
      rsss = r_1 * r_1 * (fsq + (float_8)1.44e6) / (r_2 * r_2 * (fsq + (float_8)9.61e6));

      //     take log and put the equal-loundness curve into eql_a array 
      //
      eql_a[j] = log((float_8)rsss);
      ++icount;
    }
  }

  //return without error
  //
  return ISIP_TRUE;
}