cannam@26: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@25: 
cannam@25: /*
cannam@26:     QM DSP Library
cannam@25: 
cannam@26:     Centre for Digital Music, Queen Mary, University of London.
cannam@26:     This file copyright 2005 Nicolas Chetry, copyright 2008 QMUL.
Chris@84: 
Chris@84:     This program is free software; you can redistribute it and/or
Chris@84:     modify it under the terms of the GNU General Public License as
Chris@84:     published by the Free Software Foundation; either version 2 of the
Chris@84:     License, or (at your option) any later version.  See the file
Chris@84:     COPYING included with this distribution for more information.
cannam@26: */
cannam@25: 
cannam@26: #include <cmath>
cannam@26: #include <cstdlib>
cannam@47: #include <cstring>
cannam@26: 
cannam@26: #include "MFCC.h"
cannam@26: #include "dsp/transforms/FFT.h"
cannam@26: #include "base/Window.h"
cannam@26: 
cannam@26: MFCC::MFCC(MFCCConfig config)
cannam@26: {
cannam@26:     int i,j;
cannam@26: 
cannam@26:     /* Calculate at startup */
cannam@26:     double *freqs, *lower, *center, *upper, *triangleHeight, *fftFreqs;
cannam@25:   
cannam@26:     lowestFrequency   = 66.6666666;
cannam@26:     linearFilters     = 13;
cannam@26:     linearSpacing     = 66.66666666;
cannam@26:     logFilters        = 27;
cannam@26:     logSpacing        = 1.0711703;
cannam@25:   
cannam@26:     /* FFT and analysis window sizes */
cannam@26:     fftSize           = config.fftsize;
cannam@64:     fft               = new FFTReal(fftSize);
cannam@26: 
cannam@26:     totalFilters      = linearFilters + logFilters;
cannam@30:     logPower          = config.logpower;
cannam@25:   
cannam@26:     samplingRate      = config.FS;
cannam@26:   
cannam@26:     /* The number of cepstral componenents */
cannam@26:     nceps             = config.nceps;
cannam@25: 
cannam@26:     /* Set if user want C0 */
cannam@26:     WANT_C0           = (config.want_c0 ? 1 : 0);
cannam@25:   
cannam@26:     /* Allocate space for feature vector */
cannam@26:     if (WANT_C0 == 1) {
cannam@26:         ceps              = (double*)calloc(nceps+1, sizeof(double));
cannam@26:     } else {
cannam@26:         ceps              = (double*)calloc(nceps, sizeof(double));
cannam@26:     }
cannam@26:  
cannam@26:     /* Allocate space for local vectors */
cannam@26:     mfccDCTMatrix     = (double**)calloc(nceps+1, sizeof(double*));
cannam@30:     for (i = 0; i < nceps+1; i++) {
cannam@26:         mfccDCTMatrix[i]= (double*)calloc(totalFilters, sizeof(double)); 
cannam@26:     }
cannam@26: 
cannam@26:     mfccFilterWeights = (double**)calloc(totalFilters, sizeof(double*));
cannam@30:     for (i = 0; i < totalFilters; i++) {
cannam@26:         mfccFilterWeights[i] = (double*)calloc(fftSize, sizeof(double)); 
cannam@26:     }
cannam@26:     
cannam@26:     freqs  = (double*)calloc(totalFilters+2,sizeof(double));
cannam@26:     
cannam@26:     lower  = (double*)calloc(totalFilters,sizeof(double));
cannam@26:     center = (double*)calloc(totalFilters,sizeof(double));
cannam@26:     upper  = (double*)calloc(totalFilters,sizeof(double));
cannam@26:     
cannam@26:     triangleHeight = (double*)calloc(totalFilters,sizeof(double));
cannam@26:     fftFreqs       = (double*)calloc(fftSize,sizeof(double));
cannam@25:   
cannam@30:     for (i = 0; i < linearFilters; i++) {
cannam@26:         freqs[i] = lowestFrequency + ((double)i) * linearSpacing;
cannam@26:     }
cannam@26:   
cannam@30:     for (i = linearFilters; i < totalFilters+2; i++) {
cannam@26:         freqs[i] = freqs[linearFilters-1] * 
cannam@26:             pow(logSpacing, (double)(i-linearFilters+1));
cannam@26:     }
cannam@26:   
cannam@26:     /* Define lower, center and upper */
cannam@26:     memcpy(lower,  freqs,totalFilters*sizeof(double));
cannam@26:     memcpy(center, &freqs[1],totalFilters*sizeof(double));
cannam@26:     memcpy(upper,  &freqs[2],totalFilters*sizeof(double));
cannam@26:     
cannam@26:     for (i=0;i<totalFilters;i++){
cannam@26:         triangleHeight[i] = 2./(upper[i]-lower[i]);
cannam@26:     }
cannam@26:   
cannam@26:     for (i=0;i<fftSize;i++){
cannam@26:         fftFreqs[i] = ((double) i / ((double) fftSize ) * 
cannam@26:                        (double) samplingRate);
cannam@26:     }
cannam@25: 
cannam@26:     /* Build now the mccFilterWeight matrix */
cannam@26:     for (i=0;i<totalFilters;i++){
cannam@25: 
cannam@26:         for (j=0;j<fftSize;j++) {
cannam@26:       
cannam@26:             if ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i])) {
cannam@26:           
cannam@26:                 mfccFilterWeights[i][j] = triangleHeight[i] * 
cannam@26:                     (fftFreqs[j]-lower[i]) / (center[i]-lower[i]); 
cannam@26:           
cannam@26:             }
cannam@26:             else
cannam@26:             {
cannam@26:                 mfccFilterWeights[i][j] = 0.0;
cannam@26:             }
cannam@25: 
cannam@26:             if ((fftFreqs[j]>center[i]) && (fftFreqs[j]<upper[i])) {
cannam@25: 
cannam@30:                 mfccFilterWeights[i][j] = mfccFilterWeights[i][j]
cannam@30:                     + triangleHeight[i] * (upper[i]-fftFreqs[j]) 
cannam@26:                     / (upper[i]-center[i]);
cannam@26:             }
cannam@26:             else
cannam@26:             {
cannam@26:                 mfccFilterWeights[i][j] = mfccFilterWeights[i][j] + 0.0;
cannam@26:             }
cannam@25:         }
cannam@25: 
cannam@26:     }
cannam@25: 
cannam@26:     /*
cannam@26:      * We calculate now mfccDCT matrix 
cannam@26:      * NB: +1 because of the DC component
cannam@26:      */
cannam@29: 
cannam@29:     const double pi = 3.14159265358979323846264338327950288;
cannam@26:   
cannam@30:     for (i = 0; i < nceps+1; i++) {
cannam@30:         for (j = 0; j < totalFilters; j++) {
cannam@26:             mfccDCTMatrix[i][j] = (1./sqrt((double) totalFilters / 2.))  
cannam@29:                 * cos((double) i * ((double) j + 0.5) / (double) totalFilters * pi);
cannam@25:         }
cannam@25:     }
cannam@25: 
cannam@30:     for (j = 0; j < totalFilters; j++){
cannam@30:         mfccDCTMatrix[0][j] = (sqrt(2.)/2.) * mfccDCTMatrix[0][j];
cannam@26:     }
cannam@26:    
cannam@26:     /* The analysis window */
cannam@32:     window      = new Window<double>(config.window, fftSize);
cannam@25: 
cannam@26:     /* Allocate memory for the FFT */
cannam@30:     realOut     = (double*)calloc(fftSize, sizeof(double));
cannam@30:     imagOut     = (double*)calloc(fftSize, sizeof(double));
cannam@30: 
cannam@30:     earMag      = (double*)calloc(totalFilters, sizeof(double));
cannam@30:     fftMag      = (double*)calloc(fftSize/2, sizeof(double));
cannam@25:   
cannam@26:     free(freqs);
cannam@26:     free(lower);
cannam@26:     free(center);
cannam@26:     free(upper);
cannam@26:     free(triangleHeight);
cannam@26:     free(fftFreqs);
cannam@25: }
cannam@25: 
cannam@26: MFCC::~MFCC()
cannam@26: {
cannam@26:     int i;
cannam@26:   
cannam@26:     /* Free the structure */
cannam@30:     for (i = 0; i < nceps+1; i++) {
cannam@26:         free(mfccDCTMatrix[i]);
cannam@26:     }
cannam@26:     free(mfccDCTMatrix);
cannam@26:     
cannam@30:     for (i = 0; i < totalFilters; i++) {
cannam@26:         free(mfccFilterWeights[i]);
cannam@26:     }
cannam@26:     free(mfccFilterWeights);
cannam@26:     
cannam@26:     /* Free the feature vector */
cannam@26:     free(ceps);
cannam@26:     
cannam@26:     /* The analysis window */
cannam@26:     delete window;
cannam@30: 
cannam@30:     free(earMag);
cannam@30:     free(fftMag);
cannam@26:     
cannam@26:     /* Free the FFT */
cannam@26:     free(realOut);
cannam@26:     free(imagOut);
cannam@64: 
cannam@64:     delete fft;
cannam@26: }
cannam@25: 
cannam@25: 
cannam@25: /*
cannam@25:  * 
cannam@25:  * Extract the MFCC on the input frame 
cannam@25:  * 
cannam@25:  */ 
cannam@30: int MFCC::process(const double *inframe, double *outceps)
cannam@26: {
cannam@30:     double *inputData = (double *)malloc(fftSize * sizeof(double));
cannam@30:     for (int i = 0; i < fftSize; ++i) inputData[i] = inframe[i];
cannam@25: 
cannam@30:     window->cut(inputData);
cannam@26:   
cannam@26:     /* Calculate the fft on the input frame */
cannam@64:     fft->process(0, inputData, realOut, imagOut);
cannam@25: 
cannam@30:     free(inputData);
cannam@30: 
cannam@30:     return process(realOut, imagOut, outceps);
cannam@30: }
cannam@30: 
cannam@30: int MFCC::process(const double *real, const double *imag, double *outceps)
cannam@30: {
cannam@30:     int i, j;
cannam@30: 
cannam@30:     for (i = 0; i < fftSize/2; ++i) {
cannam@30:         fftMag[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);
cannam@30:     }
cannam@30: 
cannam@30:     for (i = 0; i < totalFilters; ++i) {
cannam@30:         earMag[i] = 0.0;
cannam@25:     }
cannam@25: 
cannam@26:     /* Multiply by mfccFilterWeights */
cannam@30:     for (i = 0; i < totalFilters; i++) {
cannam@30:         double tmp = 0.0;
cannam@30:         for (j = 0; j < fftSize/2; j++) {
cannam@30:             tmp = tmp + (mfccFilterWeights[i][j] * fftMag[j]);
cannam@26:         }
cannam@30:         if (tmp > 0) earMag[i] = log10(tmp);
cannam@30: 	else earMag[i] = 0.0;
cannam@30: 
cannam@30:         if (logPower != 1.0) {
cannam@30:             earMag[i] = pow(earMag[i], logPower);
cannam@30:         }
cannam@26:     }
cannam@26: 
cannam@26:     /*
cannam@26:      * 
cannam@26:      * Calculate now the cepstral coefficients 
cannam@26:      * with or without the DC component
cannam@26:      *
cannam@26:      */
cannam@26:   
cannam@30:     if (WANT_C0 == 1) {
cannam@26:      
cannam@30:         for (i = 0; i < nceps+1; i++) {
cannam@30:             double tmp = 0.;
cannam@30:             for (j = 0; j < totalFilters; j++){
cannam@30:                 tmp = tmp + mfccDCTMatrix[i][j] * earMag[j];
cannam@26:             }
cannam@26:             outceps[i] = tmp;
cannam@26:         }
cannam@26:     }
cannam@25:     else 
cannam@26:     {  
cannam@30:         for (i = 1; i < nceps+1; i++) {
cannam@30:             double tmp = 0.;
cannam@30:             for (j = 0; j < totalFilters; j++){
cannam@30:                 tmp = tmp + mfccDCTMatrix[i][j] * earMag[j];
cannam@26:             }
cannam@26:             outceps[i-1] = tmp;
cannam@25:         }
cannam@26:     }
cannam@25:     
cannam@26:     return nceps;
cannam@25: }
cannam@25: