wav2mfcc-buffer.c

Go to the documentation of this file.

/************************************************************************/
/*    wav2mfcc.c   Convert Speech file to MFCC_E_D_(Z) file             */
/*----------------------------------------------------------------------*/
/*    Author    : Yuichiro Nakano                                       */
/*                                                                      */
/*    Copyright(C) Yuichiro Nakano 1996-1998                            */
/*----------------------------------------------------------------------*/
/************************************************************************/


#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; 
static int costbl_hamming_len = 0; 
/* cos/-sin table for FFT */
static double *costbl_fft; 
static double *sintbl_fft; 
static int tbllen = 0; 
/* cos table for MakeMFCC */
static double *costbl_makemfcc; 
static int costbl_makemfcc_len = 0; 
/* sin table for WeightCepstrum */
static double *sintbl_wcep; 
static int sintbl_wcep_len = 0; 

#endif /* MFCC_SINCOS_TABLE */

static float sqrt2var;          

#ifdef MFCC_SINCOS_TABLE
/* prepare tables */

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
}

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
}

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
}

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 */



int Wav2MFCC(SP16 *wave, float **mfcc, Value para, int nSamples, float *ssbuf, int ssbuflen)
{
  float *bf;                        /* Work space for FFT */
  double *fbank;                    /* Filterbank */
  int i, k, t;
  int end = 0, start = 1;
  int frame_num;                    /* Number of samples in output file */
  FBankInfo fb;
  float energy = 0.0;
  float c0 = 0.0;
  int p;

  /* Get filterbank information */
  fb = InitFBank(para);
  
  if((fbank = (double *)mymalloc((para.fbank_num+1)*sizeof(double))) == NULL){
    j_error("Error: Wav2MFCC: failed to malloc\n");
  }
  if((bf = (float *)mymalloc(fb.fftN * sizeof(float))) == NULL){
    j_error("Error: Wav2MFCC: failed to malloc\n");
  }

  if (ssbuf != NULL) {
    /* check ssbuf length */
    if (fb.fftN != ssbuflen) {
      j_error("Error: Wav2MFCC: noise spectrum length not match\n");
    }
  }

  frame_num = (int)((nSamples - para.framesize) / para.frameshift) + 1;

  for(t = 0; t < frame_num; t++){
    if(end != 0) start = end - (para.framesize - para.frameshift) - 1;

    k = 1;
    for(i = start; i <= start + para.framesize; i++){
      bf[k] = (float)wave[i - 1];  k++;
    }
    end = i;
    
    if (para.zmeanframe) {
      ZMeanFrame(bf, para.framesize);
    }

    /* Calculate base MFCC coefficients */

    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[t], para);
    /* weight cepstrum */
    WeightCepstrum(mfcc[t], para);
    /* set energy to mfcc */
    p = para.mfcc_dim;
    if (para.c0) mfcc[t][p++] = c0;
    if (para.energy) mfcc[t][p++] = energy;
  }
  
  /* Normalise Log Energy */
  if (para.energy && para.enormal) NormaliseLogE(mfcc, frame_num, para);
  
  /* Delta (consider energy suppress) */
  if (para.delta) Delta(mfcc, frame_num, para);

  /* Acceleration */
  if (para.acc) Accel(mfcc, frame_num, para);

  /* Cepstrum Mean Normalization */
  if(para.cmn) CMN(mfcc, frame_num, para.mfcc_dim);

  FreeFBank(fb);
  free(fbank);
  free(bf);
  return(frame_num);
}


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

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

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
}

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

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
}

float Mel(int k, float fres)
{
  return(1127 * log(1 + (k-1) * fres));
}

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

void
FreeFBank(FBankInfo fb)
{
  free(fb.cf);
  free(fb.loChan);
  free(fb.loWt);
  free(fb.Re);
  free(fb.Im);
}

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


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
}


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


void NormaliseLogE(float **mfcc, int frame_num, Value para)
{  
  float max, min, f;
  int t;
  int l;

  l = para.mfcc_dim;
  if (para.c0) l++;

  /* find max log energy */
  max = mfcc[0][l];
  for(t = 0; t < frame_num; t++)
    if(mfcc[t][l] > max) max = mfcc[t][l];

  /* set the silence floor */
  min = max - (para.silFloor * log(10.0)) / 10.0;  

  /* normalise */
  for(t = 0; t < frame_num; t++){
    f = mfcc[t][l];
    if (f < min) f = min;
    mfcc[t][l] = 1.0 - (max - f) * para.escale;
  }
}

void Delta(float **c, int frame, Value para)
{
  int theta, t, n, B = 0;
  float A1, A2, sum;
  int offset;
  float *ed;

  for(theta = 1; theta <= para.delWin; theta++)
    B += theta * theta;

  if (para.absesup) ed = (float *)mymalloc(sizeof(float) * frame);

  for(t = 0; t < frame; t++){
    for(n = 0; n < para.baselen; n++){
      sum = 0;
      for(theta = 1; theta <= para.delWin; theta++){
        /* Replicate the first or last vector */
        /* at the beginning and end of speech */
        if (t - theta < 0) A1 = c[0][n];
        else A1 = c[t - theta][n];
        if (t + theta >= frame) A2 = c[frame - 1][n];
        else A2 = c[t + theta][n];
        sum += theta * (A2 - A1);
      }
      if (para.absesup && n == para.baselen-1) {
        ed[t] = sum / (2 * B);
      } else {
        c[t][para.baselen + n] = sum / (2 * B);
      }
    }
  }
  
  if (para.absesup) {
    for (t=0;t<frame;t++) {
      memmove(&(c[t][para.baselen-1]), &(c[t][para.baselen]), sizeof(float) * para.baselen - 1);
      c[t][para.baselen * 2 - 2] = ed[t];
    }
    free(ed);
  }
}


void Accel(float **c, int frame, Value para)
{
  int theta, t, n, B = 0;
  int src, dst;
  float A1, A2, sum;

  for(theta = 1; theta <= para.accWin; theta++)
    B += theta * theta;

  for(t = 0; t < frame; t++){
    src = para.baselen * 2 - 1;
    if (para.absesup) src--;
    dst = src + para.baselen;
    for(n = 0; n < para.baselen; n++){
      sum = 0;
      for(theta = 1; theta <= para.accWin; theta++){
        /* Replicate the first or last vector */
        /* at the beginning and end of speech */
        if (t - theta < 0) A1 = c[0][src];
        else A1 = c[t - theta][src];
        if (t + theta >= frame) A2 = c[frame - 1][src];
        else A2 = c[t + theta][src];
        sum += theta * (A2 - A1);
      }
      c[t][dst] = sum / (2 * B);
      src--;
      dst--;
    }
  }
}

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


void CMN(float **mfcc, int frame_num, int dim)
{
  int i, t;
  float *mfcc_ave, *sum;

  mfcc_ave = (float *)mycalloc(dim, sizeof(float));
  sum = (float *)mycalloc(dim, sizeof(float));

  for(i = 0; i < dim; i++){
    sum[i] = 0.0;
    for(t = 0; t < frame_num; t++)
      sum[i] += mfcc[t][i];
    mfcc_ave[i] = sum[i] / frame_num;
  }
  for(t = 0; t < frame_num; t++){
    for(i = 0; i < dim; i++)
      mfcc[t][i] = mfcc[t][i] - mfcc_ave[i];
  }
  free(sum);
  free(mfcc_ave);
}

---