Templated MFCCs

4 files changed, 210 insertions, 177 deletions

diff --git a/ofxMaxim/ofxMaxim/libs/maxiFFT.h b/ofxMaxim/ofxMaxim/libs/maxiFFT.h
index 8a34c0d..d099e84 100644
--- a/ofxMaxim/ofxMaxim/libs/maxiFFT.h
+++ b/ofxMaxim/ofxMaxim/libs/maxiFFT.h
@@ -37,6 +37,7 @@
 
 
 #include "fft.h"
+#include "stddef.h"
 
 class maxiFFT {
 	
diff --git a/ofxMaxim/ofxMaxim/libs/maxiMFCC.cpp b/ofxMaxim/ofxMaxim/libs/maxiMFCC.cpp
index 6317ebb..e72fd75 100644
--- a/ofxMaxim/ofxMaxim/libs/maxiMFCC.cpp
+++ b/ofxMaxim/ofxMaxim/libs/maxiMFCC.cpp
@@ -8,51 +8,53 @@
  */
 
 #include "maxiMFCC.h"
-#include <math.h>
-#include <iostream>
-#ifdef __APPLE_CC__
-#include <Accelerate/Accelerate.h>
-#endif
 
-using namespace std;
 
-void maxiMFCC::setup(unsigned int numBins, unsigned int numFilters, unsigned int numCoeffs, double minFreq, double maxFreq, unsigned int sampleRate) 
-{
-	this->numFilters = numFilters;
-	this->numCoeffs = numCoeffs;
-	this->minFreq = minFreq;
-	this->maxFreq = maxFreq;
-	this->sampleRate = sampleRate;
-	this->numBins = numBins;
-	melFilters = NULL;
-	melBands = (double*) malloc(sizeof(double) * numFilters);
-#ifdef __APPLE_CC__
-    doubleSpec = (double*)malloc(sizeof(double) * numBins);
-#endif
-    //create new matrix
-    dctMatrix = (double*)malloc(sizeof(double) * numCoeffs * numFilters);
-	calcMelFilterBank(sampleRate, numBins);
-	createDCTCoeffs();
-}
 
-void maxiMFCC::mfcc(float* powerSpectrum, double *mfccs) {
-	melFilterAndLogSquare(powerSpectrum);
-	dct(mfccs);
+template <>
+void maxiMFCCAnalyser<double>::dct(double *mfccs) {
+	vDSP_mmulD(melBands, 1, dctMatrix, 1, mfccs, 1, 1, numCoeffs, numFilters);
+	double n = (double) numCoeffs;
+	vDSP_vsdivD(mfccs, 1, &n, mfccs, 1, numCoeffs);
 }
 
+template <>
+void maxiMFCCAnalyser<float>::dct(float *mfccs) {
+	vDSP_mmul(melBands, 1, dctMatrix, 1, mfccs, 1, 1, numCoeffs, numFilters);
+	float n = (float) numCoeffs;
+	vDSP_vsdiv(mfccs, 1, &n, mfccs, 1, numCoeffs);
+}
 
-void maxiMFCC::melFilterAndLogSquare(float* powerSpectrum){
+template <>
+void maxiMFCCAnalyser<double>::melFilterAndLogSquare(float* powerSpectrum) {
 #ifdef __APPLE_CC__
 	//conv to double
 	vDSP_vspdp(powerSpectrum, 1, doubleSpec, 1, numBins);
-	vDSP_mmulD(doubleSpec, 1, melFilters, 1, melBands, 1, 1, 42, 512);
+	//	vDSP_mmulD(doubleSpec, 1, melFilters, 1, melBands, 1, 1, 42, 512);
+	vDSP_mmulD(doubleSpec, 1, melFilters, 1, melBands, 1, 1, numFilters, numBins);
+#endif
+	melFilterAndLogSq_Part2();
+}
+
+
+template <>
+void maxiMFCCAnalyser<float>::melFilterAndLogSquare(float* powerSpectrum) {
+#ifdef __APPLE_CC__
+	vDSP_mmul(powerSpectrum, 1, melFilters, 1, melBands, 1, 1, numFilters, numBins);
+#endif
+	melFilterAndLogSq_Part2();
+}
+
+template <class T>
+void maxiMFCCAnalyser<T>::melFilterAndLogSq_Part2() {
+#ifdef __APPLE_CC__
 #else
 	for (int filter = 0;filter < numFilters;filter++) {
 		melBands[filter] = 0.0;
 		for (int bin=0;bin<numBins;bin++) {
-//			int idx = (numBins * filter) + bin;
+			//			int idx = (numBins * filter) + bin;
 			int idx = filter + (bin * numFilters);
-//			melBands[filter] += (melFilters[filter][bin] * powerSpectrum[bin]);
+			//			melBands[filter] += (melFilters[filter][bin] * powerSpectrum[bin]);
 			melBands[filter] += (melFilters[idx] * powerSpectrum[bin]);
 		}
 	}
@@ -63,144 +65,14 @@ void maxiMFCC::melFilterAndLogSquare(float* powerSpectrum){
 	}
 }
 
-void maxiMFCC::dct(double *mfccs){
-#ifdef __APPLE_CC__
-	vDSP_mmulD(melBands, 1, dctMatrix, 1, mfccs, 1, 1, numCoeffs, numFilters);
-	double n = (double) numCoeffs;
-	vDSP_vsdivD(mfccs, 1, &n, mfccs, 1, numCoeffs);
-#else
-	for(int i=0; i < numCoeffs; i++) {
-		mfccs[i] = 0.0;
-	}
-	for(int i=0; i < numCoeffs; i++ ) {
-		for(int j=0; j < numFilters; j++) {
-			int idx = i + (j * numCoeffs);
-//			mfccs[i] += (dctMatrix[i][j] * melBands[j]);
-			mfccs[i] += (dctMatrix[idx] * melBands[j]);
-		}
-	}
-	for(int i=0; i < numCoeffs; i++) {
-		mfccs[i] /= numCoeffs;
-	}
-#endif
-}
-
-void maxiMFCC::createDCTCoeffs() {
-    double k = 3.14159265358979323846/numFilters;
-    double w1 = 1.0/(sqrt(numFilters));
-    double w2 = sqrt(2.0/numFilters);
-	
-	
-    //generate dct matrix
-    for(int i = 0; i < numCoeffs; i++)
-    {
-		for(int j = 0; j < numFilters; j++)
-		{
-			int idx = i + (j * numCoeffs);
-			//			if(i == 0)
-			//				dctMatrix[i][j]= w1 * cos(k * (i+1) * (j + 0.5));
-			//			else
-			//				dctMatrix[i][j] = w2 * cos(k * (i+1) * (j + 0.5));
-			if(i == 0)
-				dctMatrix[idx]= w1 * cos(k * (i+1) * (j + 0.5));
-			else
-				dctMatrix[idx] = w2 * cos(k * (i+1) * (j + 0.5));
-		}
-    }
-	
-	
-}
 
 
 
-// implements this formula:
-// mel = 2595 log10(Hz/700 + 1)
-inline double hzToMel(double hz){
-	return 2595.0 * (log10(hz/700.0 + 1.0));
-}
 
-// implements this formula
-// Hz = 700 (10^(mel/2595) - 1)
-inline double melToHz(double mel){
-	return 700.0 * (pow(10, mel/2595.0) - 1.0);
-}
 
 
-void maxiMFCC::calcMelFilterBank(double sampleRate, int numBins){
-    double mel, dMel, maxMel, minMel, nyquist, binFreq, start, end, thisF, nextF, prevF;
-    int numValidBins;
-    
-    // ignore bins over nyquist
-    numValidBins = numBins;
-	
-    nyquist = sampleRate/2;
-    if (maxFreq > nyquist) {
-		maxFreq = nyquist;
-    }
-	
-    maxMel = hzToMel(maxFreq);
-    minMel = hzToMel(minFreq);
-	
-    dMel = (maxMel - minMel) / (numFilters + 2 - 1);
-    
-    double *filtPos = (double*) malloc(sizeof(double) * (numFilters + 2));
-	
-    // first generate an array of start and end freqs for each triangle
-    mel = minMel;
-    for (int i=0;i<numFilters + 2;i++) {
-		// start of the triangle
-		filtPos[i] = melToHz(mel);
-//		std::cout << "[" << i << "] MFCC: centre is at " <<filtPos[i]<<"hz "<<mel<<" mels" << endl;
-		mel += dMel;
-    }
-    // now generate the coefficients for the mag spectrum
-	melFilters = (double*) malloc(sizeof(double) * numFilters * numValidBins);
-	
-    for (int filter = 1; filter < numFilters; filter++) {
-		for (int bin=0;bin<numValidBins;bin++) {
-			// frequency this bin represents
-			binFreq = (double) sampleRate / (double) numValidBins * (double) bin;
-			thisF = filtPos[filter];
-			nextF = filtPos[filter+1];
-			prevF = filtPos[filter-1];
-			int idx = filter + (bin * numFilters);
-			if (binFreq > nextF || binFreq < prevF) {
-				// outside this filter
-				melFilters[idx] = 0;
-				//cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[filter][bin] << endl;
-			}
-			else {
-				double height = 2.0 / (nextF - prevF);
 
-				if (binFreq < thisF) {
-					// up
-					start = prevF;
-					end = thisF;
-					melFilters[idx] = (binFreq - start) * (height / (thisF - start));	  
-				}
-				else {
-					// down
-					start = thisF;
-					end = nextF;
-					melFilters[idx] = height + ((binFreq - thisF) * (-height /(nextF - thisF)));	  
-				}
-//				cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[filter][bin] << endl;
-				//cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[idx] << endl;
-			}
-		}
-    }
-}
 
 
-maxiMFCC::~maxiMFCC() {
-	if (melFilters) {
-		delete[] melFilters;
-		delete[] melBands;
-		delete[] dctMatrix;
-#ifdef __APPLE_CC__
-		delete doubleSpec;
-#endif
-	}
-}
 
 
diff --git a/ofxMaxim/ofxMaxim/libs/maxiMFCC.h b/ofxMaxim/ofxMaxim/libs/maxiMFCC.h
index 0f3d292..60bfe04 100644
--- a/ofxMaxim/ofxMaxim/libs/maxiMFCC.h
+++ b/ofxMaxim/ofxMaxim/libs/maxiMFCC.h
@@ -13,33 +13,192 @@
 #pragma pack(16)
 
 #include "maxiFFT.h"
+#include <math.h>
+#include <iostream>
+#ifdef __APPLE_CC__
+#include <Accelerate/Accelerate.h>
+#endif
 
-class maxiMFCC {
-public:
-	maxiMFCC():melFilters(NULL),dctMatrix(NULL), melBands(NULL){};
-	~maxiMFCC();
-	void setup(unsigned int numBins, unsigned int numFilters, unsigned int numCoeffs, double minFreq, double maxFreq, unsigned int sampleRate);
-	void mfcc(float* powerSpectrum, double *mfccs);
-	double *melBands;
+using namespace std;
+
+
+// implements this formula:
+// mel = 2595 log10(Hz/700 + 1)
+inline double hzToMel(double hz){
+	return 2595.0 * (log10(hz/700.0 + 1.0));
+}
+
+// implements this formula
+// Hz = 700 (10^(mel/2595) - 1)
+inline double melToHz(double mel){
+	return 700.0 * (pow(10, mel/2595.0) - 1.0);
+}
 
+template <class T>
+class maxiMFCCAnalyser {
+public:
+	T *melBands;
+	maxiMFCCAnalyser():melFilters(NULL),dctMatrix(NULL), melBands(NULL){};
+	~maxiMFCCAnalyser() {
+		if (melFilters) {
+			delete[] melFilters;
+			delete[] melBands;
+			delete[] dctMatrix;
+#ifdef __APPLE_CC__
+			delete doubleSpec;
+#endif
+		}
+	}
+	
+	void setup(unsigned int numBins, unsigned int numFilters, unsigned int numCoeffs, double minFreq, double maxFreq, unsigned int sampleRate) 
+	{
+		this->numFilters = numFilters;
+		this->numCoeffs = numCoeffs;
+		this->minFreq = minFreq;
+		this->maxFreq = maxFreq;
+		this->sampleRate = sampleRate;
+		this->numBins = numBins;
+		melFilters = NULL;
+		melBands = (T*) malloc(sizeof(T) * numFilters);
+#ifdef __APPLE_CC__
+		doubleSpec = (T*)malloc(sizeof(T) * numBins);
+#endif
+		//create new matrix
+		dctMatrix = (T*)malloc(sizeof(T) * numCoeffs * numFilters);
+		calcMelFilterBank(sampleRate, numBins);
+		createDCTCoeffs();
+	}
+	void mfcc(float* powerSpectrum, T *mfccs) {
+		melFilterAndLogSquare(powerSpectrum);
+		dct(mfccs);
+	}
+	
 private:
 	unsigned int numFilters, numCoeffs;
 	double minFreq, maxFreq;
 	unsigned int sampleRate;
-//	double **melFilters;
-	double *melFilters;
+	T *melFilters;
 	unsigned int numBins;
-//	double **dctMatrix;
-	double *dctMatrix;
+	T *dctMatrix;
+#ifdef __APPLE_CC__
+	T *doubleSpec;
+#endif
+	
 #ifdef __APPLE_CC__
-	double *doubleSpec;
+	void dct(T *mfccs); //define later
+#else
+	void dct(T *mfccs) {
+		for(int i=0; i < numCoeffs; i++) {
+			mfccs[i] = 0.0;
+		}
+		for(int i=0; i < numCoeffs; i++ ) {
+			for(int j=0; j < numFilters; j++) {
+				int idx = i + (j * numCoeffs);
+				mfccs[i] += (dctMatrix[idx] * melBands[j]);
+			}
+		}
+		for(int i=0; i < numCoeffs; i++) {
+			mfccs[i] /= numCoeffs;
+		}
+	}
 #endif
 	
 	void melFilterAndLogSquare(float* powerSpectrum);
-	void dct(double *mfccs);
-	void calcMelFilterBank(double sampleRate, int numBins);
-	void createDCTCoeffs();
+	void melFilterAndLogSq_Part2();
 
 	
+	void calcMelFilterBank(double sampleRate, int numBins) {
+
+		double mel, dMel, maxMel, minMel, nyquist, binFreq, start, end, thisF, nextF, prevF;
+		int numValidBins;
+		
+		// ignore bins over nyquist
+		numValidBins = numBins;
+		
+		nyquist = sampleRate/2;
+		if (maxFreq > nyquist) {
+			maxFreq = nyquist;
+		}
+		
+		maxMel = hzToMel(maxFreq);
+		minMel = hzToMel(minFreq);
+		
+		dMel = (maxMel - minMel) / (numFilters + 2 - 1);
+		
+		T *filtPos = (T*) malloc(sizeof(double) * (numFilters + 2));
+		
+		// first generate an array of start and end freqs for each triangle
+		mel = minMel;
+		for (int i=0;i<numFilters + 2;i++) {
+			// start of the triangle
+			filtPos[i] = melToHz(mel);
+			//		std::cout << "[" << i << "] MFCC: centre is at " <<filtPos[i]<<"hz "<<mel<<" mels" << endl;
+			mel += dMel;
+		}
+		// now generate the coefficients for the mag spectrum
+		melFilters = (T*) malloc(sizeof(T) * numFilters * numValidBins);
+		
+		for (int filter = 1; filter < numFilters; filter++) {
+			for (int bin=0;bin<numValidBins;bin++) {
+				// frequency this bin represents
+				binFreq = (T) sampleRate / (T) numValidBins * (T) bin;
+				thisF = filtPos[filter];
+				nextF = filtPos[filter+1];
+				prevF = filtPos[filter-1];
+				int idx = filter + (bin * numFilters);
+				if (binFreq > nextF || binFreq < prevF) {
+					// outside this filter
+					melFilters[idx] = 0;
+					//cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[filter][bin] << endl;
+				}
+				else {
+					T height = 2.0 / (nextF - prevF);
+					
+					if (binFreq < thisF) {
+						// up
+						start = prevF;
+						end = thisF;
+						melFilters[idx] = (binFreq - start) * (height / (thisF - start));	  
+					}
+					else {
+						// down
+						start = thisF;
+						end = nextF;
+						melFilters[idx] = height + ((binFreq - thisF) * (-height /(nextF - thisF)));	  
+					}
+					//				cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[filter][bin] << endl;
+					//cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[idx] << endl;
+				}
+			}
+		}
+	}
+	void createDCTCoeffs() {
+		T k = 3.14159265358979323846/numFilters;
+		T w1 = 1.0/(sqrt(numFilters));
+		T w2 = sqrt(2.0/numFilters);
+		
+		
+		//generate dct matrix
+		for(int i = 0; i < numCoeffs; i++)
+		{
+			for(int j = 0; j < numFilters; j++)
+			{
+				int idx = i + (j * numCoeffs);
+				if(i == 0)
+					dctMatrix[idx]= w1 * cos(k * (i+1) * (j + 0.5));
+				else
+					dctMatrix[idx] = w2 * cos(k * (i+1) * (j + 0.5));
+			}
+		}
+		
+		
+	}
+		
+	
 	
-};
\ No newline at end of file
+};
+
+
+
+typedef maxiMFCCAnalyser<double> maxiMFCC;
+//typedef maxiMFCCAnalyser<float> maxiFloatMFCC;
diff --git a/ofxMaxim/ofxMaxim/libs/maximilian.h b/ofxMaxim/ofxMaxim/libs/maximilian.h
index 6018992..2031940 100755
--- a/ofxMaxim/ofxMaxim/libs/maximilian.h
+++ b/ofxMaxim/ofxMaxim/libs/maximilian.h
@@ -134,6 +134,7 @@ class maxiFilter {
 	double z;//pole
 	double c;//filter coefficient
 public:
+	maxiFilter():x(0.0), y(0.0), z(0.0), c(0.0){};
 	double cutoff;
 	double resonance;
 	double lores(double input,double cutoff1, double resonance);