mfcc-core.c

about sound recognition.i want to downlod

float CalcLogRawE(float *wave, int framesize)
{		   
  int i;
  double raw_E = 0.0;
  float energy;

  for(i = 1; i <= framesize; i++)
    raw_E += wave[i] * wave[i];
  energy = (float)log(raw_E);

  return(energy);
}

/** 
 * Apply pre-emphasis filter.
 * 
 * @param wave [i/o] waveform data in the current frame
 * @param para [in] configuration parameters
 */
void PreEmphasise (float *wave, Value para)
{
  int i;
   
  for(i = para.framesize; i >= 2; i--)
    wave[i] -= wave[i - 1] * para.preEmph;
  wave[1] *= 1.0 - para.preEmph;  
}

/** 
 * Apply hamming window.
 * 
 * @param wave [i/o] waveform data in the current frame
 * @param framesize [in] frame size
 */
void Hamming(float *wave, int framesize)
{
  int i;
#ifdef MFCC_SINCOS_TABLE
  for(i = 1; i <= framesize; i++)
    wave[i] *= costbl_hamming[i-1];
#else
  float a;
  a = 2 * PI / (framesize - 1);
  for(i = 1; i <= framesize; i++)
    wave[i] *= 0.54 - 0.46 * cos(a * (i - 1));
#endif
}

/** 
 * Apply FFT
 * 
 * @param xRe [i/o] real part of waveform
 * @param xIm [i/o] imaginal part of waveform
 * @param p [in] 2^p = FFT point
 */
void FFT(float *xRe, float *xIm, int p)
{
  int i, ip, j, k, m, me, me1, n, nv2;
  double uRe, uIm, vRe, vIm, wRe, wIm, tRe, tIm;
  
  n = 1<<p;
  nv2 = n / 2;
  
  j = 0;
  for(i = 0; i < n-1; i++){
    if(j > i){
      tRe = xRe[j];      tIm = xIm[j];
      xRe[j] = xRe[i];   xIm[j] = xIm[i];
      xRe[i] = tRe;      xIm[i] = tIm;
    }
    k = nv2;
    while(j >= k){
      j -= k;      k /= 2;
    }
    j += k;
  }

  for(m = 1; m <= p; m++){
    me = 1<<m;                me1 = me / 2;
    uRe = 1.0;                uIm = 0.0;
#ifdef MFCC_SINCOS_TABLE
    wRe = costbl_fft[m-1];    wIm = sintbl_fft[m-1];
#else
    wRe = cos(PI / me1);      wIm = -sin(PI / me1);
#endif
    for(j = 0; j < me1; j++){
      for(i = j; i < n; i += me){
	ip = i + me1;
	tRe = xRe[ip] * uRe - xIm[ip] * uIm;
	tIm = xRe[ip] * uIm + xIm[ip] * uRe;
	xRe[ip] = xRe[i] - tRe;   xIm[ip] = xIm[i] - tIm;
	xRe[i] += tRe;            xIm[i] += tIm;
      }
      vRe = uRe * wRe - uIm * wIm;   vIm = uRe * wIm + uIm * wRe;
      uRe = vRe;                     uIm = vIm;
    }
  }
}


/** 
 * Convert wave -> (spectral subtraction) -> mel-frequency filterbank
 * 
 * @param wave [in] waveform data in the current frame
 * @param fbank [out] the resulting mel-frequency filterbank
 * @param fb [in] filterbank information
 * @param para [in] configuration parameters
 * @param ssbuf [in] noise spectrum, or NULL if not apply subtraction
 */
void MakeFBank(float *wave, double *fbank, FBankInfo fb, Value para, float *ssbuf)
{
  int k, bin, i;
  double Re, Im, A, P, NP, H, temp;

  for(k = 1; k <= para.framesize; k++){
    fb.Re[k - 1] = wave[k];  fb.Im[k - 1] = 0.0;  /* copy to workspace */
  }
  for(k = para.framesize + 1; k <= fb.fftN; k++){
    fb.Re[k - 1] = 0.0;      fb.Im[k - 1] = 0.0;  /* pad with zeroes */
  }
  
  /* Take FFT */
  FFT(fb.Re, fb.Im, fb.n);

  if (ssbuf != NULL) {
    /* Spectral Subtraction */
    for(k = 1; k <= fb.fftN; k++){
      Re = fb.Re[k - 1];  Im = fb.Im[k - 1];
      P = sqrt(Re * Re + Im * Im);
      NP = ssbuf[k - 1];
      if((P * P -  para.ss_alpha * NP * NP) < 0){
	H = para.ss_floor;
      }else{
	H = sqrt(P * P - para.ss_alpha * NP * NP) / P;
      }
      fb.Re[k - 1] = H * Re;
      fb.Im[k - 1] = H * Im;
    }
  }

  /* Fill filterbank channels */ 
  for(i = 1; i <= para.fbank_num; i++)
    fbank[i] = 0.0;
  
  for(k = fb.klo; k <= fb.khi; k++){
    Re = fb.Re[k-1]; Im = fb.Im[k-1];
    A = sqrt(Re * Re + Im * Im);
    bin = fb.loChan[k];
    Re = fb.loWt[k] * A;
    if(bin > 0) fbank[bin] += Re;
    if(bin < para.fbank_num) fbank[bin + 1] += A - Re;
  }

  /* Take logs */
  for(bin = 1; bin <= para.fbank_num; bin++){ 
    temp = fbank[bin];
    if(temp < 1.0) temp = 1.0;
    fbank[bin] = log(temp);  
  }
}

/** 
 * Calculate 0'th cepstral coefficient.
 * 
 * @param fbank [in] filterbank
 * @param para [in] configuration parameters
 * 
 * @return 
 */
float CalcC0(double *fbank, Value para)
{
  int i; 
  float S;
  
  S = 0.0;
  for(i = 1; i <= para.fbank_num; i++)
    S += fbank[i];
  return S * sqrt2var;
}

/** 
 * Apply DCT to filterbank to make MFCC.
 * 
 * @param fbank [in] filterbank
 * @param mfcc [out] output MFCC vector
 * @param para [in] configuration parameters
 */
void MakeMFCC(double *fbank, float *mfcc, Value para)
{
#ifdef MFCC_SINCOS_TABLE
  int i, j, k;
  k = 0;
  /* Take DCT */
  for(i = 0; i < para.mfcc_dim; i++){
    mfcc[i] = 0.0;
    for(j = 1; j <= para.fbank_num; j++)
      mfcc[i] += fbank[j] * costbl_makemfcc[k++];
    mfcc[i] *= sqrt2var;
  }
#else
  int i, j;
  float B, C;
  
  B = PI / para.fbank_num;
  /* Take DCT */
  for(i = 1; i <= para.mfcc_dim; i++){
    mfcc[i - 1] = 0.0;
    C = i * B;
    for(j = 1; j <= para.fbank_num; j++)
      mfcc[i - 1] += fbank[j] * cos(C * (j - 0.5));
    mfcc[i - 1] *= sqrt2var;     
  }
#endif
}

/** 
 * Re-scale cepstral coefficients.
 * 
 * @param mfcc [i/o] a MFCC vector
 * @param para [in] configuration parameters
 */
void WeightCepstrum (float *mfcc, Value para)
{
#ifdef MFCC_SINCOS_TABLE
  int i;
  for(i=0;i<para.mfcc_dim;i++) {
    mfcc[i] *= sintbl_wcep[i];
  }
#else
  int i;
  float a, b, *cepWin;
  
  if((cepWin = (float *)mymalloc(para.mfcc_dim * sizeof(float))) == NULL){
    j_error("WeightCepstrum: failed to malloc\n");
  }
  a = PI / para.lifter;
  b = para.lifter / 2.0;
  
  for(i = 0; i < para.mfcc_dim; i++){
    cepWin[i] = 1.0 + b * sin((i + 1) * a);
    mfcc[i] *= cepWin[i];
  }
  
  free(cepWin);
#endif
}

/************************************************************************/
/************************************************************************/
/************************************************************************/
/************************************************************************/
/************************************************************************/

static double *fbank;   ///< Local buffer to hold filterbank
static FBankInfo fb;	///< Local buffer to hold filterbank information

/** 
 * Initialize calculation functions and work areas.
 * 
 * @param para [in] configuration parameters
 * @param bf [out] returns pointer to newly allocated window buffer
 * @param bflen [out] length of @a bf
 */
void
WMP_calc_init(Value para, float **bf, int *bflen)
{
  /* Get filterbank information and initialize tables */
  fb = InitFBank(para);
  
  if((fbank = (double *)mymalloc((para.fbank_num+1)*sizeof(double))) == NULL){
    j_error("Error: WMP_calc_init failed\n");
  }

  if((*bf = (float *)mymalloc(fb.fftN * sizeof(float))) == NULL){
    j_error("Error: WMP_calc_init failed\n");
  }

  *bflen = fb.fftN;

}

/** 
 * Calculate MFCC and log energy for one frame.  Perform spectral subtraction
 * if @a ssbuf is specified.
 * 
 * @param mfcc [out] buffer to hold the resulting MFCC vector
 * @param bf [i/o] work area for FFT
 * @param para [in] configuration parameters
 * @param ssbuf [in] noise spectrum, or NULL if not using spectral subtraction
 */
void
WMP_calc(float *mfcc, float *bf, Value para, float *ssbuf)
{
  float energy = 0.0;
  float c0 = 0.0;
  int p;

  if (para.zmeanframe) {
    ZMeanFrame(bf, para.framesize);
  }

  if (para.energy && para.raw_e) {
    /* calculate log raw energy */
    energy = CalcLogRawE(bf, para.framesize);
  }
  /* pre-emphasize */
  PreEmphasise(bf, para);
  /* hamming window */
  Hamming(bf, para.framesize);
  if (para.energy && ! para.raw_e) {
    /* calculate log energy */
    energy = CalcLogRawE(bf, para.framesize);
  }
  /* filterbank */
  MakeFBank(bf, fbank, fb, para, ssbuf);
  /* 0'th cepstral parameter */
  if (para.c0) c0 = CalcC0(fbank, para);
  /* MFCC */
  MakeMFCC(fbank, mfcc, para);
  /* weight cepstrum */
  WeightCepstrum(mfcc, para);
  /* set energy to mfcc */
  p = para.mfcc_dim;
  if (para.c0) mfcc[p++] = c0;
  if (para.energy) mfcc[p++] = energy;
}

/** 
 * Free work area for MFCC computation
 * 
 * @param bf [in] window buffer previously allocated by WMP_calc_init()
 */
void
WMP_calc_fin(float *bf)
{
  FreeFBank(fb);
  free(fbank);
  free(bf);
}