view Source/MFCC.cpp @ 0:25bf17994ef1

First commit. VS2013, Codeblocks and Mac OSX configuration

/*
  ==============================================================================

    MFCC.cpp
    Created: 2 Sep 2014 3:30:50pm
    Author:  david.ronan

  ==============================================================================
*/

#include "MFCC.h"
#include <math.h>
#include <windows.h>

//-----------------------------------------------------------------------------  ComputeMFCC
void MFCC::ComputeMFCC(float* magnitude, std::vector<float> &mfccs)
{
  //Apply the Mel Filters to the spectrum magnitude:
  for(int i=0; i<m_iTotalMFCCFilters; i++)
  {  
	//Multiply spectrum with spectral mask
	vDSP_vmul(magnitude, 1,	m_ppMFCCFilters[i],	1, m_pfTempMelFilterResult, 1, m_pMFCCFilterLength[i]);

    //Sum the values of each bins
    m_fMelFilteredFFT[i] = vDSP_sve(m_pfTempMelFilterResult, 1, m_pMFCCFilterLength[i]);

    //Log compression
    float filterOut = log10(m_fMelFilteredFFT[i]*10.f+1.f);
    m_fMelFilteredFFT[i]=filterOut;
  }  

  for(int j = 0; j < m_iNumMFCCCoefs; j++)
  {
	  //Cosine Transform to reduce dimensionality:
	  vDSP_mmul(m_ppMFFC_DCT, 1, m_fMelFilteredFFT, 1, m_pMFCC, 1, m_iNumMFCCCoefs, 1, m_iTotalMFCCFilters);  
	  mfccs.push_back(m_pMFCC[j]);
  }

}


//-----------------------------------------------------------------------------  initMFFCvariables
void MFCC::initMFFCvariables(int NCoeffs, int Nfft, float fSampleRate)
{
  m_iNumMFCCCoefs = NCoeffs;

  m_pMFCC = new float[m_iNumMFCCCoefs];

  //Add a value to take into account the 0 coefficient
  int iStartMfccCoeff = 1;   

  //Mel FilterBank parameters
  float fLowestFreq       = 133.3333f;
  int   iNumLinearFilter  = 13;
  float fLinearSpacing    = 66.66666666f;
  int   iNumLogFilters    = 27;
  float fLogSpacing       = 1.0711703f;
  float fFreqPerBin       = fSampleRate / (float)(Nfft);
  m_iTotalMFCCFilters	  = iNumLinearFilter + iNumLogFilters;

  m_fMelFilteredFFT   = new float[m_iTotalMFCCFilters];
  m_pMFCCFilterStart  = new int[m_iTotalMFCCFilters];
  m_pMFCCFilterLength = new int[m_iTotalMFCCFilters];
  m_ppMFCCFilters     = new float*[m_iTotalMFCCFilters];

  float FilterLimits[3];
  int iFilterWidthMax = 0;
  for(int i = 0; i < m_iTotalMFCCFilters; i++)
  {
    for(int j = 0; j < 3; j++)
    {
      if(i + j < iNumLinearFilter)
      {
        FilterLimits[j] = (i + j) * fLinearSpacing + fLowestFreq;
      }
      else
      {
        float fLowestLogFreq = (iNumLinearFilter-1) * fLinearSpacing + fLowestFreq;
        FilterLimits[j] = fLowestLogFreq * powf(fLogSpacing, (float)((i + j) - iNumLinearFilter + 1));
      }
    }

    float fTriangleHeight = 2.f / (FilterLimits[2] - FilterLimits[0]);

    int iStartBin = (int)(ceil(FilterLimits[0] / fFreqPerBin));
    int iStopBin  = (int)(floor(FilterLimits[2] / fFreqPerBin));

    int iFilterWidth       = iStopBin-iStartBin+1;
    m_pMFCCFilterStart[i]  = iStartBin;
    m_pMFCCFilterLength[i] = iFilterWidth;
    m_ppMFCCFilters[i]     = new float[iFilterWidth];

    if(iFilterWidth > iFilterWidthMax)
	{
		iFilterWidthMax = iFilterWidth;
	}

    for(int n=iStartBin; n<iStopBin+1; n++)
    {
      float fFreq = n * fFreqPerBin;
      if(fFreq <= FilterLimits[1])
      {
        float factor =(fFreq-FilterLimits[0]) / (FilterLimits[1] - FilterLimits[0]);
        float filterVal = fTriangleHeight * factor;
        m_ppMFCCFilters[i][n - iStartBin] = filterVal;
      }
      else
      {
        float factor =(FilterLimits[2] - fFreq) / (FilterLimits[2] - FilterLimits[1]);
        float filterVal = fTriangleHeight * factor;
        m_ppMFCCFilters[i][n - iStartBin] = filterVal ;
      }
    }
  }

  m_pfTempMelFilterResult = new float[iFilterWidthMax];

  //Compute the DCT reduction matrix
  m_ppMFFC_DCT = new float[NCoeffs * m_iTotalMFCCFilters];

  float fScaleDCTMatrix = 1.f / sqrtf(m_iTotalMFCCFilters / 2.f);

  for(int n = iStartMfccCoeff; n < NCoeffs+iStartMfccCoeff; n++)
  {
    for(int m = 0; m < m_iTotalMFCCFilters; m++)
    {
      m_ppMFFC_DCT[(n - iStartMfccCoeff) * m_iTotalMFCCFilters + m] = fScaleDCTMatrix * cosf(n * (m + 0.5f) * 3.141592653589793f / m_iTotalMFCCFilters);
    }
  }
}

//-----------------------------------------------------------------------------  freeMFCCmemory
void MFCC::freeMFCCmemory()
{
	if(m_pMFCC != NULL)
	{				        
		delete(m_pMFCC);					        
		m_pMFCC=nullptr;
	}
  
	if(m_ppMFFC_DCT != NULL)
	{						
		delete(m_ppMFFC_DCT);				      
		m_ppMFFC_DCT=nullptr;
	}
  
	if(m_fMelFilteredFFT != NULL)
	{				
		delete(m_fMelFilteredFFT);		    
		m_fMelFilteredFFT=nullptr;
	}
  
	if(m_pMFCCFilterStart != NULL)
	{				
		delete(m_pMFCCFilterStart);		    
		m_pMFCCFilterStart=nullptr;
	}
  
	if(m_pMFCCFilterLength != NULL)
	{				
		delete(m_pMFCCFilterLength);		  
		m_pMFCCFilterLength=nullptr;
	}
  
	if(m_pfTempMelFilterResult != NULL)
	{
		delete(m_pfTempMelFilterResult);  
		m_pfTempMelFilterResult=nullptr;
	}

	if(m_ppMFCCFilters != NULL)
	{
		for(int i=0; i<m_iTotalMFCCFilters; i++)
		{
			if(m_ppMFCCFilters[i] != NULL)
			{ 
				delete(m_ppMFCCFilters[i]);
				m_ppMFCCFilters[i]=nullptr;
			}
		}

		delete (m_ppMFCCFilters); 
		m_ppMFCCFilters=nullptr;
	}
}

//-----------------------------------------------------------------------------  vDSP_vmul
inline void MFCC::vDSP_vmul(float* v1,int stride1, float* v2, int stride2, float* out, int outstride, int iNumSamples)
{
  for(int i=0; i<iNumSamples; i++)
  {
    out[i*outstride] = v1[i*stride1] * v2[i*stride2];
  }
}

//-----------------------------------------------------------------------------  vDSP_sve
inline float MFCC::vDSP_sve(float* v, int , int iNumSamples)
{
  float fSum = 0;

  for(int i=0; i < iNumSamples; i++)
  {
    fSum += v[i];
  }

  return fSum;
}

//-----------------------------------------------------------------------------  vDSP_mmul
inline void MFCC::vDSP_mmul(float* v1, int , float* v2, int , float* &vout, int , int M, int , int P)
{
  for(int m = 0; m < M; m++)
  {
    float fProdSum = 0.f;

    for(int p = 0; p < P; p++)
    {
      fProdSum+=v1[m * P + p] * v2[p];
    }
    
	vout[m]=fProdSum;
  }
}