mfcc-core.c

about sound recognition.i want to downlod

/**
 * @file   mfcc-core.c
 * @author Akinobu Lee
 * @date   Mon Aug  7 11:55:45 2006
 * 
 * <JA>
 * @brief  MFCC 泼魔翁の纷换
 *
 * ここでは·岭をかけて艰り叫された不兰侨妨デ〖タから MFCC 泼魔翁を
 * 换叫するコア簇眶が羌められていますˉ
 * </JA>
 * 
 * <EN>
 * @brief  Compute MFCC parameter vectors
 *
 * These are core functions to compute MFCC vectors from windowed speech data.
 * </EN>
 * 
 * $Revision: 1.3 $
 * 
 */
/*
 * Copyright (c) 1991-2006 Kawahara Lab., Kyoto University
 * Copyright (c) 2000-2005 Shikano Lab., Nara Institute of Science and Technology
 * Copyright (c) 2005-2006 Julius project team, Nagoya Institute of Technology
 * All rights reserved
 */

#include <sent/stddefs.h>
#include <sent/mfcc.h>


#ifdef MFCC_SINCOS_TABLE
/* sin/cos tables */

/* sin/cos table for hamming window */
static double *costbl_hamming; ///< Cos table for hamming window
static int costbl_hamming_len = 0; ///< Length of above
/* cos/-sin table for FFT */
static double *costbl_fft; ///< Cos table for FFT
static double *sintbl_fft; ///< Sin table for FFT
static int tbllen = 0; ///< Length of above
/* cos table for MakeMFCC */
static double *costbl_makemfcc; ///< Cos table for DCT
static int costbl_makemfcc_len = 0; ///< Length of above
/* sin table for WeightCepstrum */
static double *sintbl_wcep; ///< Sin table for cepstrum weighting
static int sintbl_wcep_len = 0; ///< Length of above

#endif /* MFCC_SINCOS_TABLE */

static float sqrt2var;		///< Work area that holds value of sqrt(2.0) / fbank_num

#ifdef MFCC_SINCOS_TABLE
/* prepare tables */

/** 
 * Generate table for hamming window.
 * 
 * @param framesize [in] window size
 */
void
make_costbl_hamming(int framesize)
{
  int i;
  float a;

  if (costbl_hamming_len  == framesize) return;
  if (costbl_hamming_len != 0) {
    j_printerr("WARNING!! frame size changed!! (%d -> %d)\n",
	       costbl_hamming_len, framesize);
    free(costbl_hamming);
  }
  costbl_hamming = (double *)mymalloc(sizeof(double) * framesize);

  a = 2.0 * PI / (framesize - 1);
  for(i=1;i<=framesize;i++) {
    /*costbl_hamming[i-1] = 0.54 - 0.46 * cos(2 * PI * (i - 1) / (float)(framesize - 1));*/
    costbl_hamming[i-1] = 0.54 - 0.46 * cos(a * (i - 1));
  }
  costbl_hamming_len = framesize;
#ifdef MFCC_TABLE_DEBUG
  j_printerr("generated Hamming cos table (%d bytes)\n",
 	     costbl_hamming_len * sizeof(double));
#endif
}

/** 
 * Build tables for FFT.
 * 
 * @param n [in] 2^n = FFT point
 */
void
make_fft_table(int n)
{
  int m;
  int me, me1;
  
  if (tbllen == n) return;	/* need not update */
  
  if (tbllen != 0) {
    j_printerr("WARNING!! fft size changed!! (%d -> %d)\n", tbllen, n);
    free(costbl_fft);
    free(sintbl_fft);
  }    
  costbl_fft = (double *)mymalloc(sizeof(double) * n);
  sintbl_fft = (double *)mymalloc(sizeof(double) * n);
  for (m = 1; m <= n; m++) {
    me = 1 << m;
    me1 = me / 2;
    costbl_fft[m-1] =  cos(PI / me1);
    sintbl_fft[m-1] = -sin(PI / me1);
  }
  tbllen = n;
#ifdef MFCC_TABLE_DEBUG
  j_printerr("generated FFT sin/cos table (%d bytes)\n", tbllen * sizeof(double));
#endif
}

/** 
 * Generate table for DCT operation to make mfcc from fbank.
 * 
 * @param fbank_num [in] number of filer banks
 * @param mfcc_dim [in] number of dimensions in MFCC
 */
void
make_costbl_makemfcc(int fbank_num, int mfcc_dim)
{
  int size;
  int i, j, k;
  float B, C;

  size = fbank_num * mfcc_dim;
  if (costbl_makemfcc_len  == size) return;
  if (costbl_makemfcc_len != 0) {
    j_printerr("WARNING!! fbank num or mfcc dim changed!! (%d -> %d)\n",
	       costbl_makemfcc_len, size);
    free(costbl_makemfcc);
  }
  costbl_makemfcc = (double *)mymalloc(sizeof(double) * size);

  B = PI / fbank_num;
  k = 0;
  for(i=1;i<=mfcc_dim;i++) {
    C = i * B;
    for(j=1;j<=fbank_num;j++) {
      costbl_makemfcc[k] = cos(C * (j - 0.5));
      k++;
    }
  }
  costbl_makemfcc_len = size;
#ifdef MFCC_TABLE_DEBUG
  j_printerr("generated MakeMFCC cos table (%d bytes)\n",
 	     costbl_makemfcc_len * sizeof(double));
#endif
}

/** 
 * Generate table for weighing cepstrum.
 * 
 * @param lifter [in] cepstral liftering coefficient
 * @param mfcc_dim [in] number of dimensions in MFCC
 */
void
make_sintbl_wcep(int lifter, int mfcc_dim)
{
  int i;
  float a, b;

  if (sintbl_wcep_len  == mfcc_dim) return;
  if (sintbl_wcep_len != 0) {
    j_printerr("WARNING!! mfcc dim changed!! (%d -> %d)\n",
	       sintbl_wcep_len, mfcc_dim);
    free(sintbl_wcep);
  }
  sintbl_wcep = (double *)mymalloc(sizeof(double) * mfcc_dim);
  a = PI / lifter;
  b = lifter / 2.0;
  for(i=0;i<mfcc_dim;i++) {
    sintbl_wcep[i] = 1.0 + b * sin((i+1) * a);
  }
  sintbl_wcep_len = mfcc_dim;
#ifdef MFCC_TABLE_DEBUG
  j_printerr("generated WeightCepstrum sin table (%d bytes)\n",
 	     sintbl_wcep_len * sizeof(double));
#endif
}

#endif /* MFCC_SINCOS_TABLE */

/** 
 * Return mel-frequency.
 * 
 * @param k [in] channel number of filter bank
 * @param fres [in] constant value computed by "1.0E7 / (para.smp_period * fb.fftN * 700.0)"
 * 
 * @return the mel frequency.
 */
float Mel(int k, float fres)
{
  return(1127 * log(1 + (k-1) * fres));
}

/** 
 * Build filterbank information and generate tables for MFCC comptutation.
 * 
 * @param para [in] configuration parameters
 * 
 * @return the generated filterbank information. 
 */
FBankInfo InitFBank(Value para)
{
  FBankInfo fb;
  float mlo, mhi, ms, melk;
  int k, chan, maxChan, nv2;

  /* Calculate FFT size */
  fb.fftN = 2;  fb.n = 1;
  while(para.framesize > fb.fftN){
    fb.fftN *= 2; fb.n++;
  }

  nv2 = fb.fftN / 2;
  fb.fres = 1.0E7 / (para.smp_period * fb.fftN * 700.0);
  maxChan = para.fbank_num + 1;
  fb.klo = 2;   fb.khi = nv2;
  mlo = 0;      mhi = Mel(nv2 + 1, fb.fres);

  /* lo pass filter */
  if (para.lopass >= 0) {
    mlo = 1127*log(1+(float)para.lopass/700.0);
    fb.klo = ((float)para.lopass * para.smp_period * 1.0e-7 * fb.fftN) + 2.5;
    if (fb.klo<2) fb.klo = 2;
  }
  /* hi pass filter */
  if (para.hipass >= 0) {
    mhi = 1127*log(1+(float)para.hipass/700.0);
    fb.khi = ((float)para.hipass * para.smp_period * 1.0e-7 * fb.fftN) + 0.5;
    if (fb.khi>nv2) fb.khi = nv2;
  }

  /* Create vector of fbank centre frequencies */
  if((fb.cf = (float *)mymalloc((maxChan + 1) * sizeof(float))) == NULL){
    j_error("InitFBank: failed to malloc\n");
  }
  ms = mhi - mlo;
  for (chan = 1; chan <= maxChan; chan++) 
    fb.cf[chan] = ((float)chan / maxChan)*ms + mlo;

  /* Create loChan map, loChan[fftindex] -> lower channel index */
  if((fb.loChan = (short *)mymalloc((nv2 + 1) * sizeof(short))) == NULL){
    j_error("InitFBank: failed to malloc\n");
  }
  for(k = 1, chan = 1; k <= nv2; k++){
    if (k < fb.klo || k > fb.khi) fb.loChan[k] = -1;
    else {
      melk = Mel(k, fb.fres);
      while (fb.cf[chan] < melk && chan <= maxChan) ++chan;
      fb.loChan[k] = chan - 1;
    }
  }

  /* Create vector of lower channel weights */   
  if((fb.loWt = (float *)mymalloc((nv2 + 1) * sizeof(float))) == NULL){
    j_error("InitFBank: failed to malloc\n");
  }
  for(k = 1; k <= nv2; k++) {
    chan = fb.loChan[k];
    if (k < fb.klo || k > fb.khi) fb.loWt[k] = 0.0;
    else {
      if (chan > 0) 
	fb.loWt[k] = (fb.cf[chan + 1] - Mel(k, fb.fres)) / (fb.cf[chan + 1] - fb.cf[chan]);
      else
	fb.loWt[k] = (fb.cf[1] - Mel(k, fb.fres)) / (fb.cf[1] - mlo);
    }
  }
  
  /* Create workspace for fft */
  if((fb.Re = (float *)mymalloc((fb.fftN + 1) * sizeof(float))) == NULL){
    j_error("InitFBank: failed to malloc\n");
  }
  if((fb.Im = (float *)mymalloc((fb.fftN + 1) * sizeof(float))) == NULL){
    j_error("InitFBank: failed to malloc\n");
  }

#ifdef MFCC_SINCOS_TABLE
  /* generate tables */
  make_costbl_hamming(para.framesize);
  make_fft_table(fb.n);
  make_costbl_makemfcc(para.fbank_num, para.mfcc_dim);
  make_sintbl_wcep(para.lifter, para.mfcc_dim);
#endif

  sqrt2var = sqrt(2.0 / para.fbank_num);

  return(fb);
}

/** 
 * Free FBankInfo.
 * 
 * @param fb [in] filterbank information
 */
void
FreeFBank(FBankInfo fb)
{
  free(fb.cf);
  free(fb.loChan);
  free(fb.loWt);
  free(fb.Re);
  free(fb.Im);
}

/** 
 * Remove DC offset per frame
 * 
 * @param wave [i/o] waveform data in the current frame
 * @param framesize [in] frame size
 * 
 */
void
ZMeanFrame(float *wave, int framesize)
{		   
  int i;
  float mean;

  mean = 0.0;
  for(i = 1; i <= framesize; i++) mean += wave[i];
  mean /= framesize;
  for(i = 1; i <= framesize; i++) wave[i] -= mean;
}

/** 
 * Calculate Log Raw Energy.
 * 
 * @param wave [in] waveform data in the current frame
 * @param framesize [in] frame size
 * 
 * @return the calculated log raw energy.
 */