cannam@25: #include <stdio.h>
cannam@25: #include <stdlib.h>
cannam@25: #include <string.h>
cannam@25: #include <math.h>
cannam@25: 
cannam@25: #include "mfcc.h"
cannam@25: #include "SBFFT.h"
cannam@25: #include "windows.h"
cannam@25: 
cannam@25: /*
cannam@25:  * 
cannam@25:  *  Initialise the MFCC structure and return a pointer to a 
cannam@25:  *  feature vector 
cannam@25:  * 
cannam@25:  */
cannam@25: 
cannam@25: extern mfcc_t* init_mfcc(int fftSize, int nceps , int samplingRate, int WANT_C0){
cannam@25: 
cannam@25:   int i,j;
cannam@25:   /* Calculate at startup */
cannam@25:   double *freqs, *lower, *center, *upper, *triangleHeight, *fftFreqs;
cannam@25:   
cannam@25:   /* Allocate space for the structure */
cannam@25:   mfcc_t* mfcc_p = (mfcc_t*)malloc(sizeof(mfcc_t));
cannam@25:   
cannam@25:   mfcc_p->lowestFrequency   = 66.6666666;
cannam@25:   mfcc_p->linearFilters     = 13;
cannam@25:   mfcc_p->linearSpacing     = 66.66666666;
cannam@25:   mfcc_p->logFilters        = 27;
cannam@25:   mfcc_p->logSpacing        = 1.0711703;
cannam@25:   
cannam@25:   /* FFT and analysis window sizes */
cannam@25:   mfcc_p->fftSize           = fftSize;
cannam@25: 
cannam@25:   mfcc_p->totalFilters      = mfcc_p->linearFilters + mfcc_p->logFilters;
cannam@25:   
cannam@25:   mfcc_p->samplingRate      = samplingRate;
cannam@25:   
cannam@25:   /* The number of cepstral componenents */
cannam@25:   mfcc_p->nceps             = nceps;
cannam@25: 
cannam@25:   /* Set if user want C0 */
cannam@25:   mfcc_p->WANT_C0           = WANT_C0;
cannam@25:   
cannam@25:   /* Allocate space for feature vector */
cannam@25:   if (mfcc_p->WANT_C0==1) {
cannam@25:     mfcc_p->ceps              = (double*)calloc(nceps+1,sizeof(double));
cannam@25:   } else {
cannam@25:     mfcc_p->ceps              = (double*)calloc(nceps,sizeof(double));
cannam@25:   }
cannam@25:  
cannam@25:   /* Allocate space for local vectors */
cannam@25:   mfcc_p->mfccDCTMatrix     = (double**)calloc(mfcc_p->nceps+1, sizeof(double*));
cannam@25:   for (i=0;i<mfcc_p->nceps+1; i++) {
cannam@25:     mfcc_p->mfccDCTMatrix[i]= (double*)calloc(mfcc_p->totalFilters, sizeof(double)); 
cannam@25:   }
cannam@25: 
cannam@25:   mfcc_p->mfccFilterWeights = (double**)calloc(mfcc_p->totalFilters, sizeof(double*));
cannam@25:   for (i=0;i<mfcc_p->totalFilters; i++) {
cannam@25:     mfcc_p->mfccFilterWeights[i] = (double*)calloc(mfcc_p->fftSize, sizeof(double)); 
cannam@25:   }
cannam@25: 
cannam@25:   freqs  = (double*)calloc(mfcc_p->totalFilters+2,sizeof(double));
cannam@25:   
cannam@25:   lower  = (double*)calloc(mfcc_p->totalFilters,sizeof(double));
cannam@25:   center = (double*)calloc(mfcc_p->totalFilters,sizeof(double));
cannam@25:   upper  = (double*)calloc(mfcc_p->totalFilters,sizeof(double));
cannam@25: 
cannam@25:   triangleHeight = (double*)calloc(mfcc_p->totalFilters,sizeof(double));
cannam@25:   fftFreqs       = (double*)calloc(mfcc_p->fftSize,sizeof(double));
cannam@25:   
cannam@25:   for (i=0;i<mfcc_p->linearFilters;i++) {
cannam@25:     freqs[i] = mfcc_p->lowestFrequency + ((double)i) * mfcc_p->linearSpacing;
cannam@25:   }
cannam@25:   
cannam@25:   for (i=mfcc_p->linearFilters; i<mfcc_p->totalFilters+2; i++) {
cannam@25:     freqs[i] = freqs[mfcc_p->linearFilters-1] * 
cannam@25:                pow(mfcc_p->logSpacing, (double)(i-mfcc_p->linearFilters+1));
cannam@25:   }
cannam@25:   
cannam@25:   /* Define lower, center and upper */
cannam@25:   memcpy(lower,  freqs,mfcc_p->totalFilters*sizeof(double));
cannam@25:   memcpy(center, &freqs[1],mfcc_p->totalFilters*sizeof(double));
cannam@25:   memcpy(upper,  &freqs[2],mfcc_p->totalFilters*sizeof(double));
cannam@25:     
cannam@25:   for (i=0;i<mfcc_p->totalFilters;i++){
cannam@25:       triangleHeight[i] = 2./(upper[i]-lower[i]);
cannam@25:   }
cannam@25:   
cannam@25:   for (i=0;i<mfcc_p->fftSize;i++){
cannam@25:       fftFreqs[i] = ((double) i / ((double) mfcc_p->fftSize ) * 
cannam@25:                      (double) mfcc_p->samplingRate);
cannam@25:   }
cannam@25: 
cannam@25:   /* Build now the mccFilterWeight matrix */
cannam@25:   for (i=0;i<mfcc_p->totalFilters;i++){
cannam@25: 
cannam@25:     for (j=0;j<mfcc_p->fftSize;j++) {
cannam@25:       
cannam@25:         if ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i])) {
cannam@25:           
cannam@25:           mfcc_p->mfccFilterWeights[i][j] = triangleHeight[i] * 
cannam@25:                                             (fftFreqs[j]-lower[i]) / (center[i]-lower[i]); 
cannam@25:           
cannam@25:         }
cannam@25:         else
cannam@25:         {
cannam@25:           
cannam@25:           mfcc_p->mfccFilterWeights[i][j] = 0.0;
cannam@25:         
cannam@25:         }
cannam@25: 
cannam@25:         if ((fftFreqs[j]>center[i]) && (fftFreqs[j]<upper[i])) {
cannam@25: 
cannam@25:           mfcc_p->mfccFilterWeights[i][j] = mfcc_p->mfccFilterWeights[i][j] + triangleHeight[i] * (upper[i]-fftFreqs[j]) 
cannam@25:                                             / (upper[i]-center[i]);
cannam@25:         }
cannam@25:         else
cannam@25:         {
cannam@25: 
cannam@25:           mfcc_p->mfccFilterWeights[i][j] = mfcc_p->mfccFilterWeights[i][j] + 0.0;
cannam@25:             
cannam@25:         }
cannam@25:     }
cannam@25: 
cannam@25:   }
cannam@25: 
cannam@25: #ifndef PI    
cannam@25: #define PI 3.14159265358979323846264338327950288
cannam@25: #endif 
cannam@25: 
cannam@25:   /*
cannam@25:    * 
cannam@25:    * We calculate now mfccDCT matrix 
cannam@25:    * NB: +1 because of the DC component
cannam@25:    * 
cannam@25:    */
cannam@25:   
cannam@25:   for (i=0; i<nceps+1; i++) {
cannam@25:     for (j=0; j<mfcc_p->totalFilters; j++) {
cannam@25:         mfcc_p->mfccDCTMatrix[i][j] = (1./sqrt((double) mfcc_p->totalFilters / 2.))  
cannam@25:                                       * cos((double) i * ((double) j + 0.5) / (double) mfcc_p->totalFilters * PI);
cannam@25:     }
cannam@25:   }
cannam@25: 
cannam@25:   for (j=0;j<mfcc_p->totalFilters;j++){
cannam@25:     mfcc_p->mfccDCTMatrix[0][j]     =  (sqrt(2.)/2.) * mfcc_p->mfccDCTMatrix[0][j];
cannam@25:   }
cannam@25:    
cannam@25:   /* The analysis window */
cannam@25:   mfcc_p->window = hamming(mfcc_p->fftSize);
cannam@25: 
cannam@25:   /* Allocate memory for the FFT */
cannam@25:   mfcc_p->imagIn      = (double*)calloc(mfcc_p->fftSize,sizeof(double));
cannam@25:   mfcc_p->realOut     = (double*)calloc(mfcc_p->fftSize,sizeof(double));
cannam@25:   mfcc_p->imagOut     = (double*)calloc(mfcc_p->fftSize,sizeof(double));
cannam@25:   
cannam@25:   free(freqs);
cannam@25:   free(lower);
cannam@25:   free(center);
cannam@25:   free(upper);
cannam@25:   free(triangleHeight);
cannam@25:   free(fftFreqs);
cannam@25: 
cannam@25:   return mfcc_p;
cannam@25: 
cannam@25: }
cannam@25: 
cannam@25: /*
cannam@25:  *
cannam@25:  *  Free the memory that has been allocated 
cannam@25:  *
cannam@25:  */ 
cannam@25: 
cannam@25: extern void close_mfcc(mfcc_t* mfcc_p) {
cannam@25: 
cannam@25:   int i;
cannam@25:   
cannam@25:   /* Free the structure */
cannam@25:   for (i=0;i<mfcc_p->nceps+1;i++) {
cannam@25:     free(mfcc_p->mfccDCTMatrix[i]);
cannam@25:   }
cannam@25:   free(mfcc_p->mfccDCTMatrix);
cannam@25: 
cannam@25:   for (i=0;i<mfcc_p->totalFilters; i++) {
cannam@25:     free(mfcc_p->mfccFilterWeights[i]);
cannam@25:   }
cannam@25:   free(mfcc_p->mfccFilterWeights);
cannam@25: 
cannam@25:   /* Free the feature vector */
cannam@25:   free(mfcc_p->ceps);
cannam@25: 
cannam@25:   /* The analysis window */
cannam@25:   free(mfcc_p->window);
cannam@25: 
cannam@25:   /* Free the FFT */
cannam@25:   free(mfcc_p->imagIn);
cannam@25:   free(mfcc_p->realOut);
cannam@25:   free(mfcc_p->imagOut);
cannam@25: 
cannam@25:   /* Free the structure itself */
cannam@25:   free(mfcc_p);
cannam@25:   mfcc_p = NULL;
cannam@25: 
cannam@25: }
cannam@25: 
cannam@25: /*
cannam@25:  * 
cannam@25:  * Extract the MFCC on the input frame 
cannam@25:  * 
cannam@25:  */ 
cannam@25: 
cannam@25: 
cannam@25: // looks like we have to have length = mfcc_p->fftSize ??????
cannam@25: 
cannam@25: extern int do_mfcc(mfcc_t* mfcc_p, double* frame, int length){
cannam@25: 
cannam@25:   int i,j;
cannam@25:   
cannam@25:   double *fftMag;
cannam@25:   double *earMag;
cannam@25: 
cannam@25:   double *inputData;
cannam@25:   
cannam@25:   double tmp;
cannam@25: 
cannam@25:   earMag    = (double*)calloc(mfcc_p->totalFilters, sizeof(double));
cannam@25:   inputData = (double*)calloc(mfcc_p->fftSize, sizeof(double)); 
cannam@25:   
cannam@25:   /* Zero-pad if needed */
cannam@25:   memcpy(inputData, frame, length*sizeof(double));
cannam@25:   
cannam@25:   /* Calculate the fft on the input frame */
cannam@25:   fft_process(mfcc_p->fftSize, 0, inputData, mfcc_p->imagIn, mfcc_p->realOut, mfcc_p->imagOut);
cannam@25: 
cannam@25:   /* Get the magnitude */
cannam@25:   fftMag = abs_fft(mfcc_p->realOut, mfcc_p->imagOut, mfcc_p->fftSize);
cannam@25: 
cannam@25:   /* Multiply by mfccFilterWeights */
cannam@25:   for (i=0;i<mfcc_p->totalFilters;i++) {
cannam@25:     tmp = 0.;
cannam@25:     for(j=0;j<mfcc_p->fftSize/2; j++) {
cannam@25:       tmp = tmp + (mfcc_p->mfccFilterWeights[i][j]*fftMag[j]);
cannam@25:     }
cannam@25:     if (tmp>0)
cannam@25: 		earMag[i] = log10(tmp);
cannam@25: 	else
cannam@25: 		earMag[i] = 0.0;
cannam@25:   }
cannam@25: 
cannam@25:   /*
cannam@25:    * 
cannam@25:    * Calculate now the ceptral coefficients 
cannam@25:    * with or without the DC component
cannam@25:    *
cannam@25:    */
cannam@25:   
cannam@25:    if (mfcc_p->WANT_C0==1) {
cannam@25:      
cannam@25:     for (i=0;i<mfcc_p->nceps+1;i++) {
cannam@25:       tmp = 0.;
cannam@25:       for (j=0;j<mfcc_p->totalFilters;j++){
cannam@25:         tmp = tmp + mfcc_p->mfccDCTMatrix[i][j]*earMag[j];
cannam@25:       }
cannam@25:       /* Send to workspace */
cannam@25:       mfcc_p->ceps[i] = tmp;
cannam@25:     }
cannam@25: 
cannam@25:    }
cannam@25:     else 
cannam@25:    {  
cannam@25:       for (i=1;i<mfcc_p->nceps+1;i++) {
cannam@25:         tmp = 0.;
cannam@25:         for (j=0;j<mfcc_p->totalFilters;j++){
cannam@25:           tmp = tmp + mfcc_p->mfccDCTMatrix[i][j]*earMag[j];
cannam@25:         }
cannam@25:         /* Send to workspace */
cannam@25:         mfcc_p->ceps[i-1] = tmp;
cannam@25:       }
cannam@25:   }
cannam@25:     
cannam@25:   free(fftMag);
cannam@25:   free(earMag);
cannam@25:   free(inputData);
cannam@25: 
cannam@25:   return mfcc_p->nceps;
cannam@25: 
cannam@25: }
cannam@25: 
cannam@25: 
cannam@25: 
cannam@25: