c@251: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@251: 
c@251: /*
c@251:     QM DSP Library
c@251: 
c@251:     Centre for Digital Music, Queen Mary, University of London.
c@251:     This file copyright 2005 Nicolas Chetry, copyright 2008 QMUL.
c@251:     All rights reserved.
c@251: */
c@251: 
c@251: #include <cmath>
c@251: #include <cstdlib>
c@251: 
c@251: #include "MFCC.h"
c@251: #include "dsp/transforms/FFT.h"
c@251: #include "base/Window.h"
c@251: 
c@251: MFCC::MFCC(MFCCConfig config)
c@251: {
c@251:     int i,j;
c@251: 
c@251:     /* Calculate at startup */
c@251:     double *freqs, *lower, *center, *upper, *triangleHeight, *fftFreqs;
c@251:   
c@251:     lowestFrequency   = 66.6666666;
c@251:     linearFilters     = 13;
c@251:     linearSpacing     = 66.66666666;
c@251:     logFilters        = 27;
c@251:     logSpacing        = 1.0711703;
c@251:   
c@251:     /* FFT and analysis window sizes */
c@251:     fftSize           = config.fftsize;
c@251: 
c@251:     totalFilters      = linearFilters + logFilters;
c@251:   
c@251:     samplingRate      = config.FS;
c@251:   
c@251:     /* The number of cepstral componenents */
c@251:     nceps             = config.nceps;
c@251: 
c@251:     /* Set if user want C0 */
c@251:     WANT_C0           = (config.want_c0 ? 1 : 0);
c@251:   
c@251:     /* Allocate space for feature vector */
c@251:     if (WANT_C0 == 1) {
c@251:         ceps              = (double*)calloc(nceps+1, sizeof(double));
c@251:     } else {
c@251:         ceps              = (double*)calloc(nceps, sizeof(double));
c@251:     }
c@251:  
c@251:     /* Allocate space for local vectors */
c@251:     mfccDCTMatrix     = (double**)calloc(nceps+1, sizeof(double*));
c@251:     for (i=0;i<nceps+1; i++) {
c@251:         mfccDCTMatrix[i]= (double*)calloc(totalFilters, sizeof(double)); 
c@251:     }
c@251: 
c@251:     mfccFilterWeights = (double**)calloc(totalFilters, sizeof(double*));
c@251:     for (i=0;i<totalFilters; i++) {
c@251:         mfccFilterWeights[i] = (double*)calloc(fftSize, sizeof(double)); 
c@251:     }
c@251:     
c@251:     freqs  = (double*)calloc(totalFilters+2,sizeof(double));
c@251:     
c@251:     lower  = (double*)calloc(totalFilters,sizeof(double));
c@251:     center = (double*)calloc(totalFilters,sizeof(double));
c@251:     upper  = (double*)calloc(totalFilters,sizeof(double));
c@251:     
c@251:     triangleHeight = (double*)calloc(totalFilters,sizeof(double));
c@251:     fftFreqs       = (double*)calloc(fftSize,sizeof(double));
c@251:   
c@251:     for (i=0;i<linearFilters;i++) {
c@251:         freqs[i] = lowestFrequency + ((double)i) * linearSpacing;
c@251:     }
c@251:   
c@251:     for (i=linearFilters; i<totalFilters+2; i++) {
c@251:         freqs[i] = freqs[linearFilters-1] * 
c@251:             pow(logSpacing, (double)(i-linearFilters+1));
c@251:     }
c@251:   
c@251:     /* Define lower, center and upper */
c@251:     memcpy(lower,  freqs,totalFilters*sizeof(double));
c@251:     memcpy(center, &freqs[1],totalFilters*sizeof(double));
c@251:     memcpy(upper,  &freqs[2],totalFilters*sizeof(double));
c@251:     
c@251:     for (i=0;i<totalFilters;i++){
c@251:         triangleHeight[i] = 2./(upper[i]-lower[i]);
c@251:     }
c@251:   
c@251:     for (i=0;i<fftSize;i++){
c@251:         fftFreqs[i] = ((double) i / ((double) fftSize ) * 
c@251:                        (double) samplingRate);
c@251:     }
c@251: 
c@251:     /* Build now the mccFilterWeight matrix */
c@251:     for (i=0;i<totalFilters;i++){
c@251: 
c@251:         for (j=0;j<fftSize;j++) {
c@251:       
c@251:             if ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i])) {
c@251:           
c@251:                 mfccFilterWeights[i][j] = triangleHeight[i] * 
c@251:                     (fftFreqs[j]-lower[i]) / (center[i]-lower[i]); 
c@251:           
c@251:             }
c@251:             else
c@251:             {
c@251:                 mfccFilterWeights[i][j] = 0.0;
c@251:             }
c@251: 
c@251:             if ((fftFreqs[j]>center[i]) && (fftFreqs[j]<upper[i])) {
c@251: 
c@251:                 mfccFilterWeights[i][j] = mfccFilterWeights[i][j] + triangleHeight[i] * (upper[i]-fftFreqs[j]) 
c@251:                     / (upper[i]-center[i]);
c@251:             }
c@251:             else
c@251:             {
c@251:                 mfccFilterWeights[i][j] = mfccFilterWeights[i][j] + 0.0;
c@251:             }
c@251:         }
c@251: 
c@251:     }
c@251: 
c@251:     /*
c@251:      * We calculate now mfccDCT matrix 
c@251:      * NB: +1 because of the DC component
c@251:      */
c@251:   
c@251:     for (i=0; i<nceps+1; i++) {
c@251:         for (j=0; j<totalFilters; j++) {
c@251:             mfccDCTMatrix[i][j] = (1./sqrt((double) totalFilters / 2.))  
c@251:                 * cos((double) i * ((double) j + 0.5) / (double) totalFilters * M_PI);
c@251:         }
c@251:     }
c@251: 
c@251:     for (j=0;j<totalFilters;j++){
c@251:         mfccDCTMatrix[0][j]     =  (sqrt(2.)/2.) * mfccDCTMatrix[0][j];
c@251:     }
c@251:    
c@251:     /* The analysis window */
c@251:     window = new Window<double>(HammingWindow, fftSize);
c@251: 
c@251:     /* Allocate memory for the FFT */
c@251:     imagIn      = (double*)calloc(fftSize,sizeof(double));
c@251:     realOut     = (double*)calloc(fftSize,sizeof(double));
c@251:     imagOut     = (double*)calloc(fftSize,sizeof(double));
c@251:   
c@251:     free(freqs);
c@251:     free(lower);
c@251:     free(center);
c@251:     free(upper);
c@251:     free(triangleHeight);
c@251:     free(fftFreqs);
c@251: }
c@251: 
c@251: MFCC::~MFCC()
c@251: {
c@251:     int i;
c@251:   
c@251:     /* Free the structure */
c@251:     for (i=0;i<nceps+1;i++) {
c@251:         free(mfccDCTMatrix[i]);
c@251:     }
c@251:     free(mfccDCTMatrix);
c@251:     
c@251:     for (i=0;i<totalFilters; i++) {
c@251:         free(mfccFilterWeights[i]);
c@251:     }
c@251:     free(mfccFilterWeights);
c@251:     
c@251:     /* Free the feature vector */
c@251:     free(ceps);
c@251:     
c@251:     /* The analysis window */
c@251:     delete window;
c@251:     
c@251:     /* Free the FFT */
c@251:     free(imagIn);
c@251:     free(realOut);
c@251:     free(imagOut);
c@251: }
c@251: 
c@251: 
c@251: /*
c@251:  * 
c@251:  * Extract the MFCC on the input frame 
c@251:  *
c@251:  * looks like we have to have length = fftSize ??????
c@251:  * 
c@251:  */ 
c@251: int MFCC::process(int length, double *inframe, double *outceps)
c@251: {
c@251:     int i,j;
c@251:   
c@251:     double *fftMag;
c@251:     double *earMag;
c@251: 
c@251:     double *inputData;
c@251:   
c@251:     double tmp;
c@251: 
c@251:     earMag    = (double*)calloc(totalFilters, sizeof(double));
c@251:     inputData = (double*)calloc(fftSize, sizeof(double)); 
c@251:     fftMag    = (double*)calloc(fftSize, sizeof(double)); 
c@251:     
c@251:     /* Zero-pad if needed */
c@251:     memcpy(inputData, inframe, length*sizeof(double));
c@251: 
c@251:     window->cut(inputData);
c@251:   
c@251:     /* Calculate the fft on the input frame */
c@251:     FFT::process(fftSize, 0, inputData, imagIn, realOut, imagOut);
c@251: 
c@251:     for (i = 0; i < fftSize; ++i) {
c@251:         fftMag[i] = sqrt(realOut[i] * realOut[i] +
c@251:                          imagOut[i] * imagOut[i]);
c@251:     }
c@251: 
c@251:     /* Multiply by mfccFilterWeights */
c@251:     for (i=0;i<totalFilters;i++) {
c@251:         tmp = 0.;
c@251:         for(j=0;j<fftSize/2; j++) {
c@251:             tmp = tmp + (mfccFilterWeights[i][j]*fftMag[j]);
c@251:         }
c@251:         if (tmp>0)
c@251:             earMag[i] = log10(tmp);
c@251: 	else
c@251:             earMag[i] = 0.0;
c@251:     }
c@251: 
c@251:     /*
c@251:      * 
c@251:      * Calculate now the cepstral coefficients 
c@251:      * with or without the DC component
c@251:      *
c@251:      */
c@251:   
c@251:     if (WANT_C0==1) {
c@251:      
c@251:         for (i=0;i<nceps+1;i++) {
c@251:             tmp = 0.;
c@251:             for (j=0;j<totalFilters;j++){
c@251:                 tmp = tmp + mfccDCTMatrix[i][j]*earMag[j];
c@251:             }
c@251:             outceps[i] = tmp;
c@251:         }
c@251:     }
c@251:     else 
c@251:     {  
c@251:         for (i=1;i<nceps+1;i++) {
c@251:             tmp = 0.;
c@251:             for (j=0;j<totalFilters;j++){
c@251:                 tmp = tmp + mfccDCTMatrix[i][j]*earMag[j];
c@251:             }
c@251:             outceps[i-1] = tmp;
c@251:         }
c@251:     }
c@251:     
c@251:     free(fftMag);
c@251:     free(earMag);
c@251:     free(inputData);
c@251: 
c@251:     return nceps;
c@251: }
c@251: