d@0: /*
d@0:   ==============================================================================
d@0: 
d@0:     MFCC.cpp
d@0:     Created: 2 Sep 2014 3:30:50pm
d@0:     Author:  david.ronan
d@0: 
d@0:   ==============================================================================
d@0: */
d@0: 
d@0: #include "MFCC.h"
d@0: #include <math.h>
d@0: #include <windows.h>
d@0: 
d@0: //-----------------------------------------------------------------------------  ComputeMFCC
d@0: void MFCC::ComputeMFCC(float* magnitude, std::vector<float> &mfccs)
d@0: {
d@0:   //Apply the Mel Filters to the spectrum magnitude:
d@0:   for(int i=0; i<m_iTotalMFCCFilters; i++)
d@0:   {  
d@0: 	//Multiply spectrum with spectral mask
d@0: 	vDSP_vmul(magnitude, 1,	m_ppMFCCFilters[i],	1, m_pfTempMelFilterResult, 1, m_pMFCCFilterLength[i]);
d@0: 
d@0:     //Sum the values of each bins
d@0:     m_fMelFilteredFFT[i] = vDSP_sve(m_pfTempMelFilterResult, 1, m_pMFCCFilterLength[i]);
d@0: 
d@0:     //Log compression
d@0:     float filterOut = log10(m_fMelFilteredFFT[i]*10.f+1.f);
d@0:     m_fMelFilteredFFT[i]=filterOut;
d@0:   }  
d@0: 
d@0:   for(int j = 0; j < m_iNumMFCCCoefs; j++)
d@0:   {
d@0: 	  //Cosine Transform to reduce dimensionality:
d@0: 	  vDSP_mmul(m_ppMFFC_DCT, 1, m_fMelFilteredFFT, 1, m_pMFCC, 1, m_iNumMFCCCoefs, 1, m_iTotalMFCCFilters);  
d@0: 	  mfccs.push_back(m_pMFCC[j]);
d@0:   }
d@0: 
d@0: }
d@0: 
d@0: 
d@0: //-----------------------------------------------------------------------------  initMFFCvariables
d@0: void MFCC::initMFFCvariables(int NCoeffs, int Nfft, float fSampleRate)
d@0: {
d@0:   m_iNumMFCCCoefs = NCoeffs;
d@0: 
d@0:   m_pMFCC = new float[m_iNumMFCCCoefs];
d@0: 
d@0:   //Add a value to take into account the 0 coefficient
d@0:   int iStartMfccCoeff = 1;   
d@0: 
d@0:   //Mel FilterBank parameters
d@0:   float fLowestFreq       = 133.3333f;
d@0:   int   iNumLinearFilter  = 13;
d@0:   float fLinearSpacing    = 66.66666666f;
d@0:   int   iNumLogFilters    = 27;
d@0:   float fLogSpacing       = 1.0711703f;
d@0:   float fFreqPerBin       = fSampleRate / (float)(Nfft);
d@0:   m_iTotalMFCCFilters	  = iNumLinearFilter + iNumLogFilters;
d@0: 
d@0:   m_fMelFilteredFFT   = new float[m_iTotalMFCCFilters];
d@0:   m_pMFCCFilterStart  = new int[m_iTotalMFCCFilters];
d@0:   m_pMFCCFilterLength = new int[m_iTotalMFCCFilters];
d@0:   m_ppMFCCFilters     = new float*[m_iTotalMFCCFilters];
d@0: 
d@0:   float FilterLimits[3];
d@0:   int iFilterWidthMax = 0;
d@0:   for(int i = 0; i < m_iTotalMFCCFilters; i++)
d@0:   {
d@0:     for(int j = 0; j < 3; j++)
d@0:     {
d@0:       if(i + j < iNumLinearFilter)
d@0:       {
d@0:         FilterLimits[j] = (i + j) * fLinearSpacing + fLowestFreq;
d@0:       }
d@0:       else
d@0:       {
d@0:         float fLowestLogFreq = (iNumLinearFilter-1) * fLinearSpacing + fLowestFreq;
d@0:         FilterLimits[j] = fLowestLogFreq * powf(fLogSpacing, (float)((i + j) - iNumLinearFilter + 1));
d@0:       }
d@0:     }
d@0: 
d@0:     float fTriangleHeight = 2.f / (FilterLimits[2] - FilterLimits[0]);
d@0: 
d@0:     int iStartBin = (int)(ceil(FilterLimits[0] / fFreqPerBin));
d@0:     int iStopBin  = (int)(floor(FilterLimits[2] / fFreqPerBin));
d@0: 
d@0:     int iFilterWidth       = iStopBin-iStartBin+1;
d@0:     m_pMFCCFilterStart[i]  = iStartBin;
d@0:     m_pMFCCFilterLength[i] = iFilterWidth;
d@0:     m_ppMFCCFilters[i]     = new float[iFilterWidth];
d@0: 
d@0:     if(iFilterWidth > iFilterWidthMax)
d@0: 	{
d@0: 		iFilterWidthMax = iFilterWidth;
d@0: 	}
d@0: 
d@0:     for(int n=iStartBin; n<iStopBin+1; n++)
d@0:     {
d@0:       float fFreq = n * fFreqPerBin;
d@0:       if(fFreq <= FilterLimits[1])
d@0:       {
d@0:         float factor =(fFreq-FilterLimits[0]) / (FilterLimits[1] - FilterLimits[0]);
d@0:         float filterVal = fTriangleHeight * factor;
d@0:         m_ppMFCCFilters[i][n - iStartBin] = filterVal;
d@0:       }
d@0:       else
d@0:       {
d@0:         float factor =(FilterLimits[2] - fFreq) / (FilterLimits[2] - FilterLimits[1]);
d@0:         float filterVal = fTriangleHeight * factor;
d@0:         m_ppMFCCFilters[i][n - iStartBin] = filterVal ;
d@0:       }
d@0:     }
d@0:   }
d@0: 
d@0:   m_pfTempMelFilterResult = new float[iFilterWidthMax];
d@0: 
d@0:   //Compute the DCT reduction matrix
d@0:   m_ppMFFC_DCT = new float[NCoeffs * m_iTotalMFCCFilters];
d@0: 
d@0:   float fScaleDCTMatrix = 1.f / sqrtf(m_iTotalMFCCFilters / 2.f);
d@0: 
d@0:   for(int n = iStartMfccCoeff; n < NCoeffs+iStartMfccCoeff; n++)
d@0:   {
d@0:     for(int m = 0; m < m_iTotalMFCCFilters; m++)
d@0:     {
d@0:       m_ppMFFC_DCT[(n - iStartMfccCoeff) * m_iTotalMFCCFilters + m] = fScaleDCTMatrix * cosf(n * (m + 0.5f) * 3.141592653589793f / m_iTotalMFCCFilters);
d@0:     }
d@0:   }
d@0: }
d@0: 
d@0: //-----------------------------------------------------------------------------  freeMFCCmemory
d@0: void MFCC::freeMFCCmemory()
d@0: {
d@0: 	if(m_pMFCC != NULL)
d@0: 	{				        
d@0: 		delete(m_pMFCC);					        
d@0: 		m_pMFCC=nullptr;
d@0: 	}
d@0:   
d@0: 	if(m_ppMFFC_DCT != NULL)
d@0: 	{						
d@0: 		delete(m_ppMFFC_DCT);				      
d@0: 		m_ppMFFC_DCT=nullptr;
d@0: 	}
d@0:   
d@0: 	if(m_fMelFilteredFFT != NULL)
d@0: 	{				
d@0: 		delete(m_fMelFilteredFFT);		    
d@0: 		m_fMelFilteredFFT=nullptr;
d@0: 	}
d@0:   
d@0: 	if(m_pMFCCFilterStart != NULL)
d@0: 	{				
d@0: 		delete(m_pMFCCFilterStart);		    
d@0: 		m_pMFCCFilterStart=nullptr;
d@0: 	}
d@0:   
d@0: 	if(m_pMFCCFilterLength != NULL)
d@0: 	{				
d@0: 		delete(m_pMFCCFilterLength);		  
d@0: 		m_pMFCCFilterLength=nullptr;
d@0: 	}
d@0:   
d@0: 	if(m_pfTempMelFilterResult != NULL)
d@0: 	{
d@0: 		delete(m_pfTempMelFilterResult);  
d@0: 		m_pfTempMelFilterResult=nullptr;
d@0: 	}
d@0: 
d@0: 	if(m_ppMFCCFilters != NULL)
d@0: 	{
d@0: 		for(int i=0; i<m_iTotalMFCCFilters; i++)
d@0: 		{
d@0: 			if(m_ppMFCCFilters[i] != NULL)
d@0: 			{ 
d@0: 				delete(m_ppMFCCFilters[i]);
d@0: 				m_ppMFCCFilters[i]=nullptr;
d@0: 			}
d@0: 		}
d@0: 
d@0: 		delete (m_ppMFCCFilters); 
d@0: 		m_ppMFCCFilters=nullptr;
d@0: 	}
d@0: }
d@0: 
d@0: //-----------------------------------------------------------------------------  vDSP_vmul
d@0: inline void MFCC::vDSP_vmul(float* v1,int stride1, float* v2, int stride2, float* out, int outstride, int iNumSamples)
d@0: {
d@0:   for(int i=0; i<iNumSamples; i++)
d@0:   {
d@0:     out[i*outstride] = v1[i*stride1] * v2[i*stride2];
d@0:   }
d@0: }
d@0: 
d@0: //-----------------------------------------------------------------------------  vDSP_sve
d@0: inline float MFCC::vDSP_sve(float* v, int , int iNumSamples)
d@0: {
d@0:   float fSum = 0;
d@0: 
d@0:   for(int i=0; i < iNumSamples; i++)
d@0:   {
d@0:     fSum += v[i];
d@0:   }
d@0: 
d@0:   return fSum;
d@0: }
d@0: 
d@0: //-----------------------------------------------------------------------------  vDSP_mmul
d@0: inline void MFCC::vDSP_mmul(float* v1, int , float* v2, int , float* &vout, int , int M, int , int P)
d@0: {
d@0:   for(int m = 0; m < M; m++)
d@0:   {
d@0:     float fProdSum = 0.f;
d@0: 
d@0:     for(int p = 0; p < P; p++)
d@0:     {
d@0:       fProdSum+=v1[m * P + p] * v2[p];
d@0:     }
d@0:     
d@0: 	vout[m]=fProdSum;
d@0:   }
d@0: }