|
wolffd@0
|
1 % mfcc - Mel frequency cepstrum coefficient analysis.
|
|
wolffd@0
|
2 % [ceps,freqresp,fb,fbrecon,freqrecon] = ...
|
|
wolffd@0
|
3 % mfcc(input, samplingRate, [frameRate])
|
|
wolffd@0
|
4 % Find the cepstral coefficients (ceps) corresponding to the
|
|
wolffd@0
|
5 % input. Four other quantities are optionally returned that
|
|
wolffd@0
|
6 % represent:
|
|
wolffd@0
|
7 % the detailed fft magnitude (freqresp) used in MFCC calculation,
|
|
wolffd@0
|
8 % the mel-scale filter bank output (fb)
|
|
wolffd@0
|
9 % the filter bank output by inverting the cepstrals with a cosine
|
|
wolffd@0
|
10 % transform (fbrecon),
|
|
wolffd@0
|
11 % the smooth frequency response by interpolating the fb reconstruction
|
|
wolffd@0
|
12 % (freqrecon)
|
|
wolffd@0
|
13 % -- Malcolm Slaney, August 1993
|
|
wolffd@0
|
14 % Modified a bit to make testing an algorithm easier... 4/15/94
|
|
wolffd@0
|
15 % Fixed Cosine Transform (indices of cos() were swapped) - 5/26/95
|
|
wolffd@0
|
16 % Added optional frameRate argument - 6/8/95
|
|
wolffd@0
|
17 % Added proper filterbank reconstruction using inverse DCT - 10/27/95
|
|
wolffd@0
|
18 % Added filterbank inversion to reconstruct spectrum - 11/1/95
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 % (c) 1998 Interval Research Corporation
|
|
wolffd@0
|
21
|
|
wolffd@0
|
22 function [ceps,freqresp,fb,fbrecon,freqrecon] = ...
|
|
wolffd@0
|
23 mfcc(input, samplingRate, frameRate)
|
|
wolffd@0
|
24 global mfccDCTMatrix mfccFilterWeights
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 [r c] = size(input);
|
|
wolffd@0
|
27 if (r > c)
|
|
wolffd@0
|
28 input=input';
|
|
wolffd@0
|
29 end
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 % Filter bank parameters
|
|
wolffd@0
|
32 lowestFrequency = 133.3333;
|
|
wolffd@0
|
33 linearFilters = 13;
|
|
wolffd@0
|
34 linearSpacing = 66.66666666;
|
|
wolffd@0
|
35 logFilters = 27;
|
|
wolffd@0
|
36 logSpacing = 1.0711703;
|
|
wolffd@0
|
37 fftSize = 512;
|
|
wolffd@0
|
38 cepstralCoefficients = 13;
|
|
wolffd@0
|
39 windowSize = 400;
|
|
wolffd@0
|
40 windowSize = 256; % Standard says 400, but 256 makes more sense
|
|
wolffd@0
|
41 % Really should be a function of the sample
|
|
wolffd@0
|
42 % rate (and the lowestFrequency) and the
|
|
wolffd@0
|
43 % frame rate.
|
|
wolffd@0
|
44 if (nargin < 2) samplingRate = 16000; end;
|
|
wolffd@0
|
45 if (nargin < 3) frameRate = 100; end;
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 % Keep this around for later....
|
|
wolffd@0
|
48 totalFilters = linearFilters + logFilters;
|
|
wolffd@0
|
49
|
|
wolffd@0
|
50 % Now figure the band edges. Interesting frequencies are spaced
|
|
wolffd@0
|
51 % by linearSpacing for a while, then go logarithmic. First figure
|
|
wolffd@0
|
52 % all the interesting frequencies. Lower, center, and upper band
|
|
wolffd@0
|
53 % edges are all consequtive interesting frequencies.
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 freqs = lowestFrequency + (0:linearFilters-1)*linearSpacing;
|
|
wolffd@0
|
56 freqs(linearFilters+1:totalFilters+2) = ...
|
|
wolffd@0
|
57 freqs(linearFilters) * logSpacing.^(1:logFilters+2);
|
|
wolffd@0
|
58
|
|
wolffd@0
|
59 lower = freqs(1:totalFilters);
|
|
wolffd@0
|
60 center = freqs(2:totalFilters+1);
|
|
wolffd@0
|
61 upper = freqs(3:totalFilters+2);
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 % We now want to combine FFT bins so that each filter has unit
|
|
wolffd@0
|
64 % weight, assuming a triangular weighting function. First figure
|
|
wolffd@0
|
65 % out the height of the triangle, then we can figure out each
|
|
wolffd@0
|
66 % frequencies contribution
|
|
wolffd@0
|
67 mfccFilterWeights = zeros(totalFilters,fftSize);
|
|
wolffd@0
|
68 triangleHeight = 2./(upper-lower);
|
|
wolffd@0
|
69 fftFreqs = (0:fftSize-1)/fftSize*samplingRate;
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71 for chan=1:totalFilters
|
|
wolffd@0
|
72 mfccFilterWeights(chan,:) = ...
|
|
wolffd@0
|
73 (fftFreqs > lower(chan) & fftFreqs <= center(chan)).* ...
|
|
wolffd@0
|
74 triangleHeight(chan).*(fftFreqs-lower(chan))/(center(chan)-lower(chan)) + ...
|
|
wolffd@0
|
75 (fftFreqs > center(chan) & fftFreqs < upper(chan)).* ...
|
|
wolffd@0
|
76 triangleHeight(chan).*(upper(chan)-fftFreqs)/(upper(chan)-center(chan));
|
|
wolffd@0
|
77 end
|
|
wolffd@0
|
78 %semilogx(fftFreqs,mfccFilterWeights')
|
|
wolffd@0
|
79 %axis([lower(1) upper(totalFilters) 0 max(max(mfccFilterWeights))])
|
|
wolffd@0
|
80
|
|
wolffd@0
|
81 hamWindow = 0.54 - 0.46*cos(2*pi*(0:windowSize-1)/windowSize);
|
|
wolffd@0
|
82
|
|
wolffd@0
|
83 if 0 % Window it like ComplexSpectrum
|
|
wolffd@0
|
84 windowStep = samplingRate/frameRate;
|
|
wolffd@0
|
85 a = .54;
|
|
wolffd@0
|
86 b = -.46;
|
|
wolffd@0
|
87 wr = sqrt(windowStep/windowSize);
|
|
wolffd@0
|
88 phi = pi/windowSize;
|
|
wolffd@0
|
89 hamWindow = 2*wr/sqrt(4*a*a+2*b*b)* ...
|
|
wolffd@0
|
90 (a + b*cos(2*pi*(0:windowSize-1)/windowSize + phi));
|
|
wolffd@0
|
91 end
|
|
wolffd@0
|
92
|
|
wolffd@0
|
93 % Figure out Discrete Cosine Transform. We want a matrix
|
|
wolffd@0
|
94 % dct(i,j) which is totalFilters x cepstralCoefficients in size.
|
|
wolffd@0
|
95 % The i,j component is given by
|
|
wolffd@0
|
96 % cos( i * (j+0.5)/totalFilters pi )
|
|
wolffd@0
|
97 % where we have assumed that i and j start at 0.
|
|
wolffd@0
|
98 mfccDCTMatrix = 1/sqrt(totalFilters/2)*cos((0:(cepstralCoefficients-1))' * ...
|
|
wolffd@0
|
99 (2*(0:(totalFilters-1))+1) * pi/2/totalFilters);
|
|
wolffd@0
|
100 mfccDCTMatrix(1,:) = mfccDCTMatrix(1,:) * sqrt(2)/2;
|
|
wolffd@0
|
101
|
|
wolffd@0
|
102 %imagesc(mfccDCTMatrix);
|
|
wolffd@0
|
103
|
|
wolffd@0
|
104 % Filter the input with the preemphasis filter. Also figure how
|
|
wolffd@0
|
105 % many columns of data we will end up with.
|
|
wolffd@0
|
106 if 1
|
|
wolffd@0
|
107 preEmphasized = filter([1 -.97], 1, input);
|
|
wolffd@0
|
108 else
|
|
wolffd@0
|
109 preEmphasized = input;
|
|
wolffd@0
|
110 end
|
|
wolffd@0
|
111 windowStep = samplingRate/frameRate;
|
|
wolffd@0
|
112 cols = fix((length(input)-windowSize)/windowStep);
|
|
wolffd@0
|
113
|
|
wolffd@0
|
114 % Allocate all the space we need for the output arrays.
|
|
wolffd@0
|
115 ceps = zeros(cepstralCoefficients, cols);
|
|
wolffd@0
|
116 if (nargout > 1) freqresp = zeros(fftSize/2, cols); end;
|
|
wolffd@0
|
117 if (nargout > 2) fb = zeros(totalFilters, cols); end;
|
|
wolffd@0
|
118
|
|
wolffd@0
|
119 % Invert the filter bank center frequencies. For each FFT bin
|
|
wolffd@0
|
120 % we want to know the exact position in the filter bank to find
|
|
wolffd@0
|
121 % the original frequency response. The next block of code finds the
|
|
wolffd@0
|
122 % integer and fractional sampling positions.
|
|
wolffd@0
|
123 if (nargout > 4)
|
|
wolffd@0
|
124 fr = (0:(fftSize/2-1))'/(fftSize/2)*samplingRate/2;
|
|
wolffd@0
|
125 j = 1;
|
|
wolffd@0
|
126 for i=1:(fftSize/2)
|
|
wolffd@0
|
127 if fr(i) > center(j+1)
|
|
wolffd@0
|
128 j = j + 1;
|
|
wolffd@0
|
129 end
|
|
wolffd@0
|
130 if j > totalFilters-1
|
|
wolffd@0
|
131 j = totalFilters-1;
|
|
wolffd@0
|
132 end
|
|
wolffd@0
|
133 fr(i) = min(totalFilters-.0001, ...
|
|
wolffd@0
|
134 max(1,j + (fr(i)-center(j))/(center(j+1)-center(j))));
|
|
wolffd@0
|
135 end
|
|
wolffd@0
|
136 fri = fix(fr);
|
|
wolffd@0
|
137 frac = fr - fri;
|
|
wolffd@0
|
138
|
|
wolffd@0
|
139 freqrecon = zeros(fftSize/2, cols);
|
|
wolffd@0
|
140 end
|
|
wolffd@0
|
141
|
|
wolffd@0
|
142 % Ok, now let's do the processing. For each chunk of data:
|
|
wolffd@0
|
143 % * Window the data with a hamming window,
|
|
wolffd@0
|
144 % * Shift it into FFT order,
|
|
wolffd@0
|
145 % * Find the magnitude of the fft,
|
|
wolffd@0
|
146 % * Convert the fft data into filter bank outputs,
|
|
wolffd@0
|
147 % * Find the log base 10,
|
|
wolffd@0
|
148 % * Find the cosine transform to reduce dimensionality.
|
|
wolffd@0
|
149 for start=0:cols-1
|
|
wolffd@0
|
150 first = start*windowStep + 1;
|
|
wolffd@0
|
151 last = first + windowSize-1;
|
|
wolffd@0
|
152 fftData = zeros(1,fftSize);
|
|
wolffd@0
|
153 fftData(1:windowSize) = preEmphasized(first:last).*hamWindow;
|
|
wolffd@0
|
154 fftMag = abs(fft(fftData));
|
|
wolffd@0
|
155 earMag = log10(mfccFilterWeights * fftMag');
|
|
wolffd@0
|
156
|
|
wolffd@0
|
157 ceps(:,start+1) = mfccDCTMatrix * earMag;
|
|
wolffd@0
|
158 if (nargout > 1) freqresp(:,start+1) = fftMag(1:fftSize/2)'; end;
|
|
wolffd@0
|
159 if (nargout > 2) fb(:,start+1) = earMag; end
|
|
wolffd@0
|
160 if (nargout > 3)
|
|
wolffd@0
|
161 fbrecon(:,start+1) = ...
|
|
wolffd@0
|
162 mfccDCTMatrix(1:cepstralCoefficients,:)' * ...
|
|
wolffd@0
|
163 ceps(:,start+1);
|
|
wolffd@0
|
164 end
|
|
wolffd@0
|
165 if (nargout > 4)
|
|
wolffd@0
|
166 f10 = 10.^fbrecon(:,start+1);
|
|
wolffd@0
|
167 freqrecon(:,start+1) = samplingRate/fftSize * ...
|
|
wolffd@0
|
168 (f10(fri).*(1-frac) + f10(fri+1).*frac);
|
|
wolffd@0
|
169 end
|
|
wolffd@0
|
170 end
|
|
wolffd@0
|
171
|
|
wolffd@0
|
172 % OK, just to check things, let's also reconstruct the original FB
|
|
wolffd@0
|
173 % output. We do this by multiplying the cepstral data by the transpose
|
|
wolffd@0
|
174 % of the original DCT matrix. This all works because we were careful to
|
|
wolffd@0
|
175 % scale the DCT matrix so it was orthonormal.
|
|
wolffd@0
|
176 if 1 & (nargout > 3)
|
|
wolffd@0
|
177 fbrecon = mfccDCTMatrix(1:cepstralCoefficients,:)' * ceps;
|
|
wolffd@0
|
178 % imagesc(mt(:,1:cepstralCoefficients)*mfccDCTMatrix);
|
|
wolffd@0
|
179 end;
|
|
wolffd@0
|
180
|
|
wolffd@0
|
181
|
