# HG changeset patch # User cannam # Date 1199978168 0 # Node ID 54a9627272716ca9c0b0f01e9c56ef40033b3650 # Parent 2b74bd60c61f2ebd1de4f3d1d91581460434cb14 ... diff -r 2b74bd60c61f -r 54a962727271 dsp/mfcc/MFCC.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/mfcc/MFCC.c Thu Jan 10 15:16:08 2008 +0000 @@ -0,0 +1,292 @@ +#include +#include +#include +#include + +#include "mfcc.h" +#include "SBFFT.h" +#include "windows.h" + +/* + * + * Initialise the MFCC structure and return a pointer to a + * feature vector + * + */ + +extern mfcc_t* init_mfcc(int fftSize, int nceps , int samplingRate, int WANT_C0){ + + int i,j; + /* Calculate at startup */ + double *freqs, *lower, *center, *upper, *triangleHeight, *fftFreqs; + + /* Allocate space for the structure */ + mfcc_t* mfcc_p = (mfcc_t*)malloc(sizeof(mfcc_t)); + + mfcc_p->lowestFrequency = 66.6666666; + mfcc_p->linearFilters = 13; + mfcc_p->linearSpacing = 66.66666666; + mfcc_p->logFilters = 27; + mfcc_p->logSpacing = 1.0711703; + + /* FFT and analysis window sizes */ + mfcc_p->fftSize = fftSize; + + mfcc_p->totalFilters = mfcc_p->linearFilters + mfcc_p->logFilters; + + mfcc_p->samplingRate = samplingRate; + + /* The number of cepstral componenents */ + mfcc_p->nceps = nceps; + + /* Set if user want C0 */ + mfcc_p->WANT_C0 = WANT_C0; + + /* Allocate space for feature vector */ + if (mfcc_p->WANT_C0==1) { + mfcc_p->ceps = (double*)calloc(nceps+1,sizeof(double)); + } else { + mfcc_p->ceps = (double*)calloc(nceps,sizeof(double)); + } + + /* Allocate space for local vectors */ + mfcc_p->mfccDCTMatrix = (double**)calloc(mfcc_p->nceps+1, sizeof(double*)); + for (i=0;inceps+1; i++) { + mfcc_p->mfccDCTMatrix[i]= (double*)calloc(mfcc_p->totalFilters, sizeof(double)); + } + + mfcc_p->mfccFilterWeights = (double**)calloc(mfcc_p->totalFilters, sizeof(double*)); + for (i=0;itotalFilters; i++) { + mfcc_p->mfccFilterWeights[i] = (double*)calloc(mfcc_p->fftSize, sizeof(double)); + } + + freqs = (double*)calloc(mfcc_p->totalFilters+2,sizeof(double)); + + lower = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); + center = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); + upper = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); + + triangleHeight = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); + fftFreqs = (double*)calloc(mfcc_p->fftSize,sizeof(double)); + + for (i=0;ilinearFilters;i++) { + freqs[i] = mfcc_p->lowestFrequency + ((double)i) * mfcc_p->linearSpacing; + } + + for (i=mfcc_p->linearFilters; itotalFilters+2; i++) { + freqs[i] = freqs[mfcc_p->linearFilters-1] * + pow(mfcc_p->logSpacing, (double)(i-mfcc_p->linearFilters+1)); + } + + /* Define lower, center and upper */ + memcpy(lower, freqs,mfcc_p->totalFilters*sizeof(double)); + memcpy(center, &freqs[1],mfcc_p->totalFilters*sizeof(double)); + memcpy(upper, &freqs[2],mfcc_p->totalFilters*sizeof(double)); + + for (i=0;itotalFilters;i++){ + triangleHeight[i] = 2./(upper[i]-lower[i]); + } + + for (i=0;ifftSize;i++){ + fftFreqs[i] = ((double) i / ((double) mfcc_p->fftSize ) * + (double) mfcc_p->samplingRate); + } + + /* Build now the mccFilterWeight matrix */ + for (i=0;itotalFilters;i++){ + + for (j=0;jfftSize;j++) { + + if ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i])) { + + mfcc_p->mfccFilterWeights[i][j] = triangleHeight[i] * + (fftFreqs[j]-lower[i]) / (center[i]-lower[i]); + + } + else + { + + mfcc_p->mfccFilterWeights[i][j] = 0.0; + + } + + if ((fftFreqs[j]>center[i]) && (fftFreqs[j]mfccFilterWeights[i][j] = mfcc_p->mfccFilterWeights[i][j] + triangleHeight[i] * (upper[i]-fftFreqs[j]) + / (upper[i]-center[i]); + } + else + { + + mfcc_p->mfccFilterWeights[i][j] = mfcc_p->mfccFilterWeights[i][j] + 0.0; + + } + } + + } + +#ifndef PI +#define PI 3.14159265358979323846264338327950288 +#endif + + /* + * + * We calculate now mfccDCT matrix + * NB: +1 because of the DC component + * + */ + + for (i=0; itotalFilters; j++) { + mfcc_p->mfccDCTMatrix[i][j] = (1./sqrt((double) mfcc_p->totalFilters / 2.)) + * cos((double) i * ((double) j + 0.5) / (double) mfcc_p->totalFilters * PI); + } + } + + for (j=0;jtotalFilters;j++){ + mfcc_p->mfccDCTMatrix[0][j] = (sqrt(2.)/2.) * mfcc_p->mfccDCTMatrix[0][j]; + } + + /* The analysis window */ + mfcc_p->window = hamming(mfcc_p->fftSize); + + /* Allocate memory for the FFT */ + mfcc_p->imagIn = (double*)calloc(mfcc_p->fftSize,sizeof(double)); + mfcc_p->realOut = (double*)calloc(mfcc_p->fftSize,sizeof(double)); + mfcc_p->imagOut = (double*)calloc(mfcc_p->fftSize,sizeof(double)); + + free(freqs); + free(lower); + free(center); + free(upper); + free(triangleHeight); + free(fftFreqs); + + return mfcc_p; + +} + +/* + * + * Free the memory that has been allocated + * + */ + +extern void close_mfcc(mfcc_t* mfcc_p) { + + int i; + + /* Free the structure */ + for (i=0;inceps+1;i++) { + free(mfcc_p->mfccDCTMatrix[i]); + } + free(mfcc_p->mfccDCTMatrix); + + for (i=0;itotalFilters; i++) { + free(mfcc_p->mfccFilterWeights[i]); + } + free(mfcc_p->mfccFilterWeights); + + /* Free the feature vector */ + free(mfcc_p->ceps); + + /* The analysis window */ + free(mfcc_p->window); + + /* Free the FFT */ + free(mfcc_p->imagIn); + free(mfcc_p->realOut); + free(mfcc_p->imagOut); + + /* Free the structure itself */ + free(mfcc_p); + mfcc_p = NULL; + +} + +/* + * + * Extract the MFCC on the input frame + * + */ + + +// looks like we have to have length = mfcc_p->fftSize ?????? + +extern int do_mfcc(mfcc_t* mfcc_p, double* frame, int length){ + + int i,j; + + double *fftMag; + double *earMag; + + double *inputData; + + double tmp; + + earMag = (double*)calloc(mfcc_p->totalFilters, sizeof(double)); + inputData = (double*)calloc(mfcc_p->fftSize, sizeof(double)); + + /* Zero-pad if needed */ + memcpy(inputData, frame, length*sizeof(double)); + + /* Calculate the fft on the input frame */ + fft_process(mfcc_p->fftSize, 0, inputData, mfcc_p->imagIn, mfcc_p->realOut, mfcc_p->imagOut); + + /* Get the magnitude */ + fftMag = abs_fft(mfcc_p->realOut, mfcc_p->imagOut, mfcc_p->fftSize); + + /* Multiply by mfccFilterWeights */ + for (i=0;itotalFilters;i++) { + tmp = 0.; + for(j=0;jfftSize/2; j++) { + tmp = tmp + (mfcc_p->mfccFilterWeights[i][j]*fftMag[j]); + } + if (tmp>0) + earMag[i] = log10(tmp); + else + earMag[i] = 0.0; + } + + /* + * + * Calculate now the ceptral coefficients + * with or without the DC component + * + */ + + if (mfcc_p->WANT_C0==1) { + + for (i=0;inceps+1;i++) { + tmp = 0.; + for (j=0;jtotalFilters;j++){ + tmp = tmp + mfcc_p->mfccDCTMatrix[i][j]*earMag[j]; + } + /* Send to workspace */ + mfcc_p->ceps[i] = tmp; + } + + } + else + { + for (i=1;inceps+1;i++) { + tmp = 0.; + for (j=0;jtotalFilters;j++){ + tmp = tmp + mfcc_p->mfccDCTMatrix[i][j]*earMag[j]; + } + /* Send to workspace */ + mfcc_p->ceps[i-1] = tmp; + } + } + + free(fftMag); + free(earMag); + free(inputData); + + return mfcc_p->nceps; + +} + + + + diff -r 2b74bd60c61f -r 54a962727271 dsp/mfcc/MFCC.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/mfcc/MFCC.h Thu Jan 10 15:16:08 2008 +0000 @@ -0,0 +1,51 @@ +#ifndef _LIB_MFCC_H +#define _LIB_MFCC_H + +#define MFCC 6 + +typedef struct mfcc_t { + + /* Filter bank parameters */ + double lowestFrequency; + int linearFilters; + double linearSpacing; + int logFilters; + double logSpacing; + + /* FFT length */ + int fftSize; + + /* Analysis window length*/ + int windowSize; + + int totalFilters; + + /* Misc. */ + int samplingRate; + int nceps; + + /* MFCC vector */ + double *ceps; + + double **mfccDCTMatrix; + double **mfccFilterWeights; + + /* The analysis window */ + double *window; + + /* For the FFT */ + double* imagIn; // always zero + double* realOut; + double* imagOut; + + /* Set if user want C0 */ + int WANT_C0; + +} mfcc_t; + +extern mfcc_t* init_mfcc(int fftSize, int nceps , int samplingRate, int WANT_C0); +extern int do_mfcc(mfcc_t* mfcc_p, double* frame, int length); +extern void close_mfcc(mfcc_t* mfcc_p); + +#endif + diff -r 2b74bd60c61f -r 54a962727271 dsp/mfcc/mfcc.c --- a/dsp/mfcc/mfcc.c Thu Jan 10 15:14:53 2008 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,292 +0,0 @@ -#include -#include -#include -#include - -#include "mfcc.h" -#include "SBFFT.h" -#include "windows.h" - -/* - * - * Initialise the MFCC structure and return a pointer to a - * feature vector - * - */ - -extern mfcc_t* init_mfcc(int fftSize, int nceps , int samplingRate, int WANT_C0){ - - int i,j; - /* Calculate at startup */ - double *freqs, *lower, *center, *upper, *triangleHeight, *fftFreqs; - - /* Allocate space for the structure */ - mfcc_t* mfcc_p = (mfcc_t*)malloc(sizeof(mfcc_t)); - - mfcc_p->lowestFrequency = 66.6666666; - mfcc_p->linearFilters = 13; - mfcc_p->linearSpacing = 66.66666666; - mfcc_p->logFilters = 27; - mfcc_p->logSpacing = 1.0711703; - - /* FFT and analysis window sizes */ - mfcc_p->fftSize = fftSize; - - mfcc_p->totalFilters = mfcc_p->linearFilters + mfcc_p->logFilters; - - mfcc_p->samplingRate = samplingRate; - - /* The number of cepstral componenents */ - mfcc_p->nceps = nceps; - - /* Set if user want C0 */ - mfcc_p->WANT_C0 = WANT_C0; - - /* Allocate space for feature vector */ - if (mfcc_p->WANT_C0==1) { - mfcc_p->ceps = (double*)calloc(nceps+1,sizeof(double)); - } else { - mfcc_p->ceps = (double*)calloc(nceps,sizeof(double)); - } - - /* Allocate space for local vectors */ - mfcc_p->mfccDCTMatrix = (double**)calloc(mfcc_p->nceps+1, sizeof(double*)); - for (i=0;inceps+1; i++) { - mfcc_p->mfccDCTMatrix[i]= (double*)calloc(mfcc_p->totalFilters, sizeof(double)); - } - - mfcc_p->mfccFilterWeights = (double**)calloc(mfcc_p->totalFilters, sizeof(double*)); - for (i=0;itotalFilters; i++) { - mfcc_p->mfccFilterWeights[i] = (double*)calloc(mfcc_p->fftSize, sizeof(double)); - } - - freqs = (double*)calloc(mfcc_p->totalFilters+2,sizeof(double)); - - lower = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); - center = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); - upper = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); - - triangleHeight = (double*)calloc(mfcc_p->totalFilters,sizeof(double)); - fftFreqs = (double*)calloc(mfcc_p->fftSize,sizeof(double)); - - for (i=0;ilinearFilters;i++) { - freqs[i] = mfcc_p->lowestFrequency + ((double)i) * mfcc_p->linearSpacing; - } - - for (i=mfcc_p->linearFilters; itotalFilters+2; i++) { - freqs[i] = freqs[mfcc_p->linearFilters-1] * - pow(mfcc_p->logSpacing, (double)(i-mfcc_p->linearFilters+1)); - } - - /* Define lower, center and upper */ - memcpy(lower, freqs,mfcc_p->totalFilters*sizeof(double)); - memcpy(center, &freqs[1],mfcc_p->totalFilters*sizeof(double)); - memcpy(upper, &freqs[2],mfcc_p->totalFilters*sizeof(double)); - - for (i=0;itotalFilters;i++){ - triangleHeight[i] = 2./(upper[i]-lower[i]); - } - - for (i=0;ifftSize;i++){ - fftFreqs[i] = ((double) i / ((double) mfcc_p->fftSize ) * - (double) mfcc_p->samplingRate); - } - - /* Build now the mccFilterWeight matrix */ - for (i=0;itotalFilters;i++){ - - for (j=0;jfftSize;j++) { - - if ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i])) { - - mfcc_p->mfccFilterWeights[i][j] = triangleHeight[i] * - (fftFreqs[j]-lower[i]) / (center[i]-lower[i]); - - } - else - { - - mfcc_p->mfccFilterWeights[i][j] = 0.0; - - } - - if ((fftFreqs[j]>center[i]) && (fftFreqs[j]mfccFilterWeights[i][j] = mfcc_p->mfccFilterWeights[i][j] + triangleHeight[i] * (upper[i]-fftFreqs[j]) - / (upper[i]-center[i]); - } - else - { - - mfcc_p->mfccFilterWeights[i][j] = mfcc_p->mfccFilterWeights[i][j] + 0.0; - - } - } - - } - -#ifndef PI -#define PI 3.14159265358979323846264338327950288 -#endif - - /* - * - * We calculate now mfccDCT matrix - * NB: +1 because of the DC component - * - */ - - for (i=0; itotalFilters; j++) { - mfcc_p->mfccDCTMatrix[i][j] = (1./sqrt((double) mfcc_p->totalFilters / 2.)) - * cos((double) i * ((double) j + 0.5) / (double) mfcc_p->totalFilters * PI); - } - } - - for (j=0;jtotalFilters;j++){ - mfcc_p->mfccDCTMatrix[0][j] = (sqrt(2.)/2.) * mfcc_p->mfccDCTMatrix[0][j]; - } - - /* The analysis window */ - mfcc_p->window = hamming(mfcc_p->fftSize); - - /* Allocate memory for the FFT */ - mfcc_p->imagIn = (double*)calloc(mfcc_p->fftSize,sizeof(double)); - mfcc_p->realOut = (double*)calloc(mfcc_p->fftSize,sizeof(double)); - mfcc_p->imagOut = (double*)calloc(mfcc_p->fftSize,sizeof(double)); - - free(freqs); - free(lower); - free(center); - free(upper); - free(triangleHeight); - free(fftFreqs); - - return mfcc_p; - -} - -/* - * - * Free the memory that has been allocated - * - */ - -extern void close_mfcc(mfcc_t* mfcc_p) { - - int i; - - /* Free the structure */ - for (i=0;inceps+1;i++) { - free(mfcc_p->mfccDCTMatrix[i]); - } - free(mfcc_p->mfccDCTMatrix); - - for (i=0;itotalFilters; i++) { - free(mfcc_p->mfccFilterWeights[i]); - } - free(mfcc_p->mfccFilterWeights); - - /* Free the feature vector */ - free(mfcc_p->ceps); - - /* The analysis window */ - free(mfcc_p->window); - - /* Free the FFT */ - free(mfcc_p->imagIn); - free(mfcc_p->realOut); - free(mfcc_p->imagOut); - - /* Free the structure itself */ - free(mfcc_p); - mfcc_p = NULL; - -} - -/* - * - * Extract the MFCC on the input frame - * - */ - - -// looks like we have to have length = mfcc_p->fftSize ?????? - -extern int do_mfcc(mfcc_t* mfcc_p, double* frame, int length){ - - int i,j; - - double *fftMag; - double *earMag; - - double *inputData; - - double tmp; - - earMag = (double*)calloc(mfcc_p->totalFilters, sizeof(double)); - inputData = (double*)calloc(mfcc_p->fftSize, sizeof(double)); - - /* Zero-pad if needed */ - memcpy(inputData, frame, length*sizeof(double)); - - /* Calculate the fft on the input frame */ - fft_process(mfcc_p->fftSize, 0, inputData, mfcc_p->imagIn, mfcc_p->realOut, mfcc_p->imagOut); - - /* Get the magnitude */ - fftMag = abs_fft(mfcc_p->realOut, mfcc_p->imagOut, mfcc_p->fftSize); - - /* Multiply by mfccFilterWeights */ - for (i=0;itotalFilters;i++) { - tmp = 0.; - for(j=0;jfftSize/2; j++) { - tmp = tmp + (mfcc_p->mfccFilterWeights[i][j]*fftMag[j]); - } - if (tmp>0) - earMag[i] = log10(tmp); - else - earMag[i] = 0.0; - } - - /* - * - * Calculate now the ceptral coefficients - * with or without the DC component - * - */ - - if (mfcc_p->WANT_C0==1) { - - for (i=0;inceps+1;i++) { - tmp = 0.; - for (j=0;jtotalFilters;j++){ - tmp = tmp + mfcc_p->mfccDCTMatrix[i][j]*earMag[j]; - } - /* Send to workspace */ - mfcc_p->ceps[i] = tmp; - } - - } - else - { - for (i=1;inceps+1;i++) { - tmp = 0.; - for (j=0;jtotalFilters;j++){ - tmp = tmp + mfcc_p->mfccDCTMatrix[i][j]*earMag[j]; - } - /* Send to workspace */ - mfcc_p->ceps[i-1] = tmp; - } - } - - free(fftMag); - free(earMag); - free(inputData); - - return mfcc_p->nceps; - -} - - - - diff -r 2b74bd60c61f -r 54a962727271 dsp/mfcc/mfcc.h --- a/dsp/mfcc/mfcc.h Thu Jan 10 15:14:53 2008 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,51 +0,0 @@ -#ifndef _LIB_MFCC_H -#define _LIB_MFCC_H - -#define MFCC 6 - -typedef struct mfcc_t { - - /* Filter bank parameters */ - double lowestFrequency; - int linearFilters; - double linearSpacing; - int logFilters; - double logSpacing; - - /* FFT length */ - int fftSize; - - /* Analysis window length*/ - int windowSize; - - int totalFilters; - - /* Misc. */ - int samplingRate; - int nceps; - - /* MFCC vector */ - double *ceps; - - double **mfccDCTMatrix; - double **mfccFilterWeights; - - /* The analysis window */ - double *window; - - /* For the FFT */ - double* imagIn; // always zero - double* realOut; - double* imagOut; - - /* Set if user want C0 */ - int WANT_C0; - -} mfcc_t; - -extern mfcc_t* init_mfcc(int fftSize, int nceps , int samplingRate, int WANT_C0); -extern int do_mfcc(mfcc_t* mfcc_p, double* frame, int length); -extern void close_mfcc(mfcc_t* mfcc_p); - -#endif -