view Source/MFCC.cpp @ 1:e86e9c111b29

Updates stuff that potentially fixes the memory leak and also makes it work on Windows and Linux (Need to test). Still have to fix fftw include for linux in Jucer.

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

    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 != nullptr)
	{
		delete(m_pMFCC);
		m_pMFCC=nullptr;
	}

	if(m_ppMFFC_DCT != nullptr)
	{
		delete(m_ppMFFC_DCT);
		m_ppMFFC_DCT=nullptr;
	}

	if(m_fMelFilteredFFT != nullptr)
	{
		delete(m_fMelFilteredFFT);
		m_fMelFilteredFFT=nullptr;
	}

	if(m_pMFCCFilterStart != nullptr)
	{
		delete(m_pMFCCFilterStart);
		m_pMFCCFilterStart=nullptr;
	}

	if(m_pMFCCFilterLength != nullptr)
	{
		delete(m_pMFCCFilterLength);
		m_pMFCCFilterLength=nullptr;
	}

	if(m_pfTempMelFilterResult != nullptr)
	{
		delete(m_pfTempMelFilterResult);
		m_pfTempMelFilterResult=nullptr;
	}

	if(m_ppMFCCFilters != nullptr)
	{
		for(int i=0; i<m_iTotalMFCCFilters; i++)
		{
			if(m_ppMFCCFilters[i] != nullptr)
			{
				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;
  }
}