d@0
|
1 /*
|
d@0
|
2 ==============================================================================
|
d@0
|
3
|
d@0
|
4 MFCC.cpp
|
d@0
|
5 Created: 2 Sep 2014 3:30:50pm
|
d@0
|
6 Author: david.ronan
|
d@0
|
7
|
d@0
|
8 ==============================================================================
|
d@0
|
9 */
|
d@0
|
10
|
d@0
|
11 #include "MFCC.h"
|
d@0
|
12 #include <math.h>
|
d@0
|
13 #include <windows.h>
|
d@0
|
14
|
d@0
|
15 //----------------------------------------------------------------------------- ComputeMFCC
|
d@0
|
16 void MFCC::ComputeMFCC(float* magnitude, std::vector<float> &mfccs)
|
d@0
|
17 {
|
d@0
|
18 //Apply the Mel Filters to the spectrum magnitude:
|
d@0
|
19 for(int i=0; i<m_iTotalMFCCFilters; i++)
|
d@0
|
20 {
|
d@0
|
21 //Multiply spectrum with spectral mask
|
d@0
|
22 vDSP_vmul(magnitude, 1, m_ppMFCCFilters[i], 1, m_pfTempMelFilterResult, 1, m_pMFCCFilterLength[i]);
|
d@0
|
23
|
d@0
|
24 //Sum the values of each bins
|
d@0
|
25 m_fMelFilteredFFT[i] = vDSP_sve(m_pfTempMelFilterResult, 1, m_pMFCCFilterLength[i]);
|
d@0
|
26
|
d@0
|
27 //Log compression
|
d@0
|
28 float filterOut = log10(m_fMelFilteredFFT[i]*10.f+1.f);
|
d@0
|
29 m_fMelFilteredFFT[i]=filterOut;
|
d@0
|
30 }
|
d@0
|
31
|
d@0
|
32 for(int j = 0; j < m_iNumMFCCCoefs; j++)
|
d@0
|
33 {
|
d@0
|
34 //Cosine Transform to reduce dimensionality:
|
d@0
|
35 vDSP_mmul(m_ppMFFC_DCT, 1, m_fMelFilteredFFT, 1, m_pMFCC, 1, m_iNumMFCCCoefs, 1, m_iTotalMFCCFilters);
|
d@0
|
36 mfccs.push_back(m_pMFCC[j]);
|
d@0
|
37 }
|
d@0
|
38
|
d@0
|
39 }
|
d@0
|
40
|
d@0
|
41
|
d@0
|
42 //----------------------------------------------------------------------------- initMFFCvariables
|
d@0
|
43 void MFCC::initMFFCvariables(int NCoeffs, int Nfft, float fSampleRate)
|
d@0
|
44 {
|
d@0
|
45 m_iNumMFCCCoefs = NCoeffs;
|
d@0
|
46
|
d@0
|
47 m_pMFCC = new float[m_iNumMFCCCoefs];
|
d@0
|
48
|
d@0
|
49 //Add a value to take into account the 0 coefficient
|
d@0
|
50 int iStartMfccCoeff = 1;
|
d@0
|
51
|
d@0
|
52 //Mel FilterBank parameters
|
d@0
|
53 float fLowestFreq = 133.3333f;
|
d@0
|
54 int iNumLinearFilter = 13;
|
d@0
|
55 float fLinearSpacing = 66.66666666f;
|
d@0
|
56 int iNumLogFilters = 27;
|
d@0
|
57 float fLogSpacing = 1.0711703f;
|
d@0
|
58 float fFreqPerBin = fSampleRate / (float)(Nfft);
|
d@0
|
59 m_iTotalMFCCFilters = iNumLinearFilter + iNumLogFilters;
|
d@0
|
60
|
d@0
|
61 m_fMelFilteredFFT = new float[m_iTotalMFCCFilters];
|
d@0
|
62 m_pMFCCFilterStart = new int[m_iTotalMFCCFilters];
|
d@0
|
63 m_pMFCCFilterLength = new int[m_iTotalMFCCFilters];
|
d@0
|
64 m_ppMFCCFilters = new float*[m_iTotalMFCCFilters];
|
d@0
|
65
|
d@0
|
66 float FilterLimits[3];
|
d@0
|
67 int iFilterWidthMax = 0;
|
d@0
|
68 for(int i = 0; i < m_iTotalMFCCFilters; i++)
|
d@0
|
69 {
|
d@0
|
70 for(int j = 0; j < 3; j++)
|
d@0
|
71 {
|
d@0
|
72 if(i + j < iNumLinearFilter)
|
d@0
|
73 {
|
d@0
|
74 FilterLimits[j] = (i + j) * fLinearSpacing + fLowestFreq;
|
d@0
|
75 }
|
d@0
|
76 else
|
d@0
|
77 {
|
d@0
|
78 float fLowestLogFreq = (iNumLinearFilter-1) * fLinearSpacing + fLowestFreq;
|
d@0
|
79 FilterLimits[j] = fLowestLogFreq * powf(fLogSpacing, (float)((i + j) - iNumLinearFilter + 1));
|
d@0
|
80 }
|
d@0
|
81 }
|
d@0
|
82
|
d@0
|
83 float fTriangleHeight = 2.f / (FilterLimits[2] - FilterLimits[0]);
|
d@0
|
84
|
d@0
|
85 int iStartBin = (int)(ceil(FilterLimits[0] / fFreqPerBin));
|
d@0
|
86 int iStopBin = (int)(floor(FilterLimits[2] / fFreqPerBin));
|
d@0
|
87
|
d@0
|
88 int iFilterWidth = iStopBin-iStartBin+1;
|
d@0
|
89 m_pMFCCFilterStart[i] = iStartBin;
|
d@0
|
90 m_pMFCCFilterLength[i] = iFilterWidth;
|
d@0
|
91 m_ppMFCCFilters[i] = new float[iFilterWidth];
|
d@0
|
92
|
d@0
|
93 if(iFilterWidth > iFilterWidthMax)
|
d@0
|
94 {
|
d@0
|
95 iFilterWidthMax = iFilterWidth;
|
d@0
|
96 }
|
d@0
|
97
|
d@0
|
98 for(int n=iStartBin; n<iStopBin+1; n++)
|
d@0
|
99 {
|
d@0
|
100 float fFreq = n * fFreqPerBin;
|
d@0
|
101 if(fFreq <= FilterLimits[1])
|
d@0
|
102 {
|
d@0
|
103 float factor =(fFreq-FilterLimits[0]) / (FilterLimits[1] - FilterLimits[0]);
|
d@0
|
104 float filterVal = fTriangleHeight * factor;
|
d@0
|
105 m_ppMFCCFilters[i][n - iStartBin] = filterVal;
|
d@0
|
106 }
|
d@0
|
107 else
|
d@0
|
108 {
|
d@0
|
109 float factor =(FilterLimits[2] - fFreq) / (FilterLimits[2] - FilterLimits[1]);
|
d@0
|
110 float filterVal = fTriangleHeight * factor;
|
d@0
|
111 m_ppMFCCFilters[i][n - iStartBin] = filterVal ;
|
d@0
|
112 }
|
d@0
|
113 }
|
d@0
|
114 }
|
d@0
|
115
|
d@0
|
116 m_pfTempMelFilterResult = new float[iFilterWidthMax];
|
d@0
|
117
|
d@0
|
118 //Compute the DCT reduction matrix
|
d@0
|
119 m_ppMFFC_DCT = new float[NCoeffs * m_iTotalMFCCFilters];
|
d@0
|
120
|
d@0
|
121 float fScaleDCTMatrix = 1.f / sqrtf(m_iTotalMFCCFilters / 2.f);
|
d@0
|
122
|
d@0
|
123 for(int n = iStartMfccCoeff; n < NCoeffs+iStartMfccCoeff; n++)
|
d@0
|
124 {
|
d@0
|
125 for(int m = 0; m < m_iTotalMFCCFilters; m++)
|
d@0
|
126 {
|
d@0
|
127 m_ppMFFC_DCT[(n - iStartMfccCoeff) * m_iTotalMFCCFilters + m] = fScaleDCTMatrix * cosf(n * (m + 0.5f) * 3.141592653589793f / m_iTotalMFCCFilters);
|
d@0
|
128 }
|
d@0
|
129 }
|
d@0
|
130 }
|
d@0
|
131
|
d@0
|
132 //----------------------------------------------------------------------------- freeMFCCmemory
|
d@0
|
133 void MFCC::freeMFCCmemory()
|
d@0
|
134 {
|
d@0
|
135 if(m_pMFCC != NULL)
|
d@0
|
136 {
|
d@0
|
137 delete(m_pMFCC);
|
d@0
|
138 m_pMFCC=nullptr;
|
d@0
|
139 }
|
d@0
|
140
|
d@0
|
141 if(m_ppMFFC_DCT != NULL)
|
d@0
|
142 {
|
d@0
|
143 delete(m_ppMFFC_DCT);
|
d@0
|
144 m_ppMFFC_DCT=nullptr;
|
d@0
|
145 }
|
d@0
|
146
|
d@0
|
147 if(m_fMelFilteredFFT != NULL)
|
d@0
|
148 {
|
d@0
|
149 delete(m_fMelFilteredFFT);
|
d@0
|
150 m_fMelFilteredFFT=nullptr;
|
d@0
|
151 }
|
d@0
|
152
|
d@0
|
153 if(m_pMFCCFilterStart != NULL)
|
d@0
|
154 {
|
d@0
|
155 delete(m_pMFCCFilterStart);
|
d@0
|
156 m_pMFCCFilterStart=nullptr;
|
d@0
|
157 }
|
d@0
|
158
|
d@0
|
159 if(m_pMFCCFilterLength != NULL)
|
d@0
|
160 {
|
d@0
|
161 delete(m_pMFCCFilterLength);
|
d@0
|
162 m_pMFCCFilterLength=nullptr;
|
d@0
|
163 }
|
d@0
|
164
|
d@0
|
165 if(m_pfTempMelFilterResult != NULL)
|
d@0
|
166 {
|
d@0
|
167 delete(m_pfTempMelFilterResult);
|
d@0
|
168 m_pfTempMelFilterResult=nullptr;
|
d@0
|
169 }
|
d@0
|
170
|
d@0
|
171 if(m_ppMFCCFilters != NULL)
|
d@0
|
172 {
|
d@0
|
173 for(int i=0; i<m_iTotalMFCCFilters; i++)
|
d@0
|
174 {
|
d@0
|
175 if(m_ppMFCCFilters[i] != NULL)
|
d@0
|
176 {
|
d@0
|
177 delete(m_ppMFCCFilters[i]);
|
d@0
|
178 m_ppMFCCFilters[i]=nullptr;
|
d@0
|
179 }
|
d@0
|
180 }
|
d@0
|
181
|
d@0
|
182 delete (m_ppMFCCFilters);
|
d@0
|
183 m_ppMFCCFilters=nullptr;
|
d@0
|
184 }
|
d@0
|
185 }
|
d@0
|
186
|
d@0
|
187 //----------------------------------------------------------------------------- vDSP_vmul
|
d@0
|
188 inline void MFCC::vDSP_vmul(float* v1,int stride1, float* v2, int stride2, float* out, int outstride, int iNumSamples)
|
d@0
|
189 {
|
d@0
|
190 for(int i=0; i<iNumSamples; i++)
|
d@0
|
191 {
|
d@0
|
192 out[i*outstride] = v1[i*stride1] * v2[i*stride2];
|
d@0
|
193 }
|
d@0
|
194 }
|
d@0
|
195
|
d@0
|
196 //----------------------------------------------------------------------------- vDSP_sve
|
d@0
|
197 inline float MFCC::vDSP_sve(float* v, int , int iNumSamples)
|
d@0
|
198 {
|
d@0
|
199 float fSum = 0;
|
d@0
|
200
|
d@0
|
201 for(int i=0; i < iNumSamples; i++)
|
d@0
|
202 {
|
d@0
|
203 fSum += v[i];
|
d@0
|
204 }
|
d@0
|
205
|
d@0
|
206 return fSum;
|
d@0
|
207 }
|
d@0
|
208
|
d@0
|
209 //----------------------------------------------------------------------------- vDSP_mmul
|
d@0
|
210 inline void MFCC::vDSP_mmul(float* v1, int , float* v2, int , float* &vout, int , int M, int , int P)
|
d@0
|
211 {
|
d@0
|
212 for(int m = 0; m < M; m++)
|
d@0
|
213 {
|
d@0
|
214 float fProdSum = 0.f;
|
d@0
|
215
|
d@0
|
216 for(int p = 0; p < P; p++)
|
d@0
|
217 {
|
d@0
|
218 fProdSum+=v1[m * P + p] * v2[p];
|
d@0
|
219 }
|
d@0
|
220
|
d@0
|
221 vout[m]=fProdSum;
|
d@0
|
222 }
|
d@0
|
223 } |