wolffd@0: %  mfcc - Mel frequency cepstrum coefficient analysis.
wolffd@0: %   [ceps,freqresp,fb,fbrecon,freqrecon] = ...
wolffd@0: %			mfcc(input, samplingRate, [frameRate])
wolffd@0: % Find the cepstral coefficients (ceps) corresponding to the
wolffd@0: % input.  Four other quantities are optionally returned that
wolffd@0: % represent:
wolffd@0: %	the detailed fft magnitude (freqresp) used in MFCC calculation, 
wolffd@0: %	the mel-scale filter bank output (fb)
wolffd@0: %	the filter bank output by inverting the cepstrals with a cosine 
wolffd@0: %		transform (fbrecon),
wolffd@0: %	the smooth frequency response by interpolating the fb reconstruction 
wolffd@0: %		(freqrecon)
wolffd@0: %  -- Malcolm Slaney, August 1993
wolffd@0: % Modified a bit to make testing an algorithm easier... 4/15/94
wolffd@0: % Fixed Cosine Transform (indices of cos() were swapped) - 5/26/95
wolffd@0: % Added optional frameRate argument - 6/8/95
wolffd@0: % Added proper filterbank reconstruction using inverse DCT - 10/27/95
wolffd@0: % Added filterbank inversion to reconstruct spectrum - 11/1/95
wolffd@0: 
wolffd@0: % (c) 1998 Interval Research Corporation  
wolffd@0: 
wolffd@0: function [ceps,freqresp,fb,fbrecon,freqrecon] = ...
wolffd@0: 		mfcc(input, samplingRate, frameRate)
wolffd@0: global mfccDCTMatrix mfccFilterWeights
wolffd@0: 
wolffd@0: [r c] = size(input);
wolffd@0: if (r > c) 
wolffd@0: 	input=input';
wolffd@0: end
wolffd@0: 
wolffd@0: %	Filter bank parameters
wolffd@0: lowestFrequency = 133.3333;
wolffd@0: linearFilters = 13;
wolffd@0: linearSpacing = 66.66666666;
wolffd@0: logFilters = 27;
wolffd@0: logSpacing = 1.0711703;
wolffd@0: fftSize = 512;
wolffd@0: cepstralCoefficients = 13;
wolffd@0: windowSize = 400;
wolffd@0: windowSize = 256;		% Standard says 400, but 256 makes more sense
wolffd@0: 				% Really should be a function of the sample
wolffd@0: 				% rate (and the lowestFrequency) and the
wolffd@0: 				% frame rate.
wolffd@0: if (nargin < 2) samplingRate = 16000; end;
wolffd@0: if (nargin < 3) frameRate = 100; end;
wolffd@0: 
wolffd@0: % Keep this around for later....
wolffd@0: totalFilters = linearFilters + logFilters;
wolffd@0: 
wolffd@0: % Now figure the band edges.  Interesting frequencies are spaced
wolffd@0: % by linearSpacing for a while, then go logarithmic.  First figure
wolffd@0: % all the interesting frequencies.  Lower, center, and upper band
wolffd@0: % edges are all consequtive interesting frequencies. 
wolffd@0: 
wolffd@0: freqs = lowestFrequency + (0:linearFilters-1)*linearSpacing;
wolffd@0: freqs(linearFilters+1:totalFilters+2) = ...
wolffd@0: 		      freqs(linearFilters) * logSpacing.^(1:logFilters+2);
wolffd@0: 
wolffd@0: lower = freqs(1:totalFilters);
wolffd@0: center = freqs(2:totalFilters+1);
wolffd@0: upper = freqs(3:totalFilters+2);
wolffd@0: 
wolffd@0: % We now want to combine FFT bins so that each filter has unit
wolffd@0: % weight, assuming a triangular weighting function.  First figure
wolffd@0: % out the height of the triangle, then we can figure out each 
wolffd@0: % frequencies contribution
wolffd@0: mfccFilterWeights = zeros(totalFilters,fftSize);
wolffd@0: triangleHeight = 2./(upper-lower);
wolffd@0: fftFreqs = (0:fftSize-1)/fftSize*samplingRate;
wolffd@0: 
wolffd@0: for chan=1:totalFilters
wolffd@0: 	mfccFilterWeights(chan,:) = ...
wolffd@0:   (fftFreqs > lower(chan) & fftFreqs <= center(chan)).* ...
wolffd@0:    triangleHeight(chan).*(fftFreqs-lower(chan))/(center(chan)-lower(chan)) + ...
wolffd@0:   (fftFreqs > center(chan) & fftFreqs < upper(chan)).* ...
wolffd@0:    triangleHeight(chan).*(upper(chan)-fftFreqs)/(upper(chan)-center(chan));
wolffd@0: end
wolffd@0: %semilogx(fftFreqs,mfccFilterWeights')
wolffd@0: %axis([lower(1) upper(totalFilters) 0 max(max(mfccFilterWeights))])
wolffd@0: 
wolffd@0: hamWindow = 0.54 - 0.46*cos(2*pi*(0:windowSize-1)/windowSize);
wolffd@0: 
wolffd@0: if 0					% Window it like ComplexSpectrum
wolffd@0: 	windowStep = samplingRate/frameRate;
wolffd@0: 	a = .54;
wolffd@0: 	b = -.46;
wolffd@0: 	wr = sqrt(windowStep/windowSize);
wolffd@0: 	phi = pi/windowSize;
wolffd@0: 	hamWindow = 2*wr/sqrt(4*a*a+2*b*b)* ...
wolffd@0: 		(a + b*cos(2*pi*(0:windowSize-1)/windowSize + phi));
wolffd@0: end
wolffd@0: 
wolffd@0: % Figure out Discrete Cosine Transform.  We want a matrix
wolffd@0: % dct(i,j) which is totalFilters x cepstralCoefficients in size.
wolffd@0: % The i,j component is given by 
wolffd@0: %                cos( i * (j+0.5)/totalFilters pi )
wolffd@0: % where we have assumed that i and j start at 0.
wolffd@0: mfccDCTMatrix = 1/sqrt(totalFilters/2)*cos((0:(cepstralCoefficients-1))' * ...
wolffd@0: 				(2*(0:(totalFilters-1))+1) * pi/2/totalFilters);
wolffd@0: mfccDCTMatrix(1,:) = mfccDCTMatrix(1,:) * sqrt(2)/2;
wolffd@0: 
wolffd@0: %imagesc(mfccDCTMatrix);
wolffd@0: 
wolffd@0: % Filter the input with the preemphasis filter.  Also figure how
wolffd@0: % many columns of data we will end up with.
wolffd@0: if 1
wolffd@0: 	preEmphasized = filter([1 -.97], 1, input);
wolffd@0: else
wolffd@0: 	preEmphasized = input;
wolffd@0: end
wolffd@0: windowStep = samplingRate/frameRate;
wolffd@0: cols = fix((length(input)-windowSize)/windowStep);
wolffd@0: 
wolffd@0: % Allocate all the space we need for the output arrays.
wolffd@0: ceps = zeros(cepstralCoefficients, cols);
wolffd@0: if (nargout > 1) freqresp = zeros(fftSize/2, cols); end;
wolffd@0: if (nargout > 2) fb = zeros(totalFilters, cols); end;
wolffd@0: 
wolffd@0: % Invert the filter bank center frequencies.  For each FFT bin
wolffd@0: % we want to know the exact position in the filter bank to find
wolffd@0: % the original frequency response.  The next block of code finds the
wolffd@0: % integer and fractional sampling positions.
wolffd@0: if (nargout > 4)
wolffd@0: 	fr = (0:(fftSize/2-1))'/(fftSize/2)*samplingRate/2;
wolffd@0: 	j = 1;
wolffd@0: 	for i=1:(fftSize/2)
wolffd@0: 		if fr(i) > center(j+1)
wolffd@0: 			j = j + 1;
wolffd@0: 		end
wolffd@0: 		if j > totalFilters-1
wolffd@0: 			j = totalFilters-1;
wolffd@0: 		end
wolffd@0: 		fr(i) = min(totalFilters-.0001, ...
wolffd@0: 		    max(1,j + (fr(i)-center(j))/(center(j+1)-center(j))));
wolffd@0: 	end
wolffd@0: 	fri = fix(fr);
wolffd@0: 	frac = fr - fri;
wolffd@0: 
wolffd@0: 	freqrecon = zeros(fftSize/2, cols);
wolffd@0: end
wolffd@0: 
wolffd@0: % Ok, now let's do the processing.  For each chunk of data:
wolffd@0: %    * Window the data with a hamming window,
wolffd@0: %    * Shift it into FFT order,
wolffd@0: %    * Find the magnitude of the fft,
wolffd@0: %    * Convert the fft data into filter bank outputs,
wolffd@0: %    * Find the log base 10,
wolffd@0: %    * Find the cosine transform to reduce dimensionality.
wolffd@0: for start=0:cols-1
wolffd@0:     first = start*windowStep + 1;
wolffd@0:     last = first + windowSize-1;
wolffd@0:     fftData = zeros(1,fftSize);
wolffd@0:     fftData(1:windowSize) = preEmphasized(first:last).*hamWindow;
wolffd@0:     fftMag = abs(fft(fftData));
wolffd@0:     earMag = log10(mfccFilterWeights * fftMag');
wolffd@0: 
wolffd@0:     ceps(:,start+1) = mfccDCTMatrix * earMag;
wolffd@0:     if (nargout > 1) freqresp(:,start+1) = fftMag(1:fftSize/2)'; end;
wolffd@0:     if (nargout > 2) fb(:,start+1) = earMag; end
wolffd@0: 	if (nargout > 3) 
wolffd@0: 		fbrecon(:,start+1) = ...
wolffd@0: 			mfccDCTMatrix(1:cepstralCoefficients,:)' * ...
wolffd@0: 			ceps(:,start+1);
wolffd@0: 	end
wolffd@0: 	if (nargout > 4)
wolffd@0: 		f10 = 10.^fbrecon(:,start+1);
wolffd@0: 		freqrecon(:,start+1) = samplingRate/fftSize * ...
wolffd@0: 			(f10(fri).*(1-frac) + f10(fri+1).*frac);
wolffd@0: 	end
wolffd@0: end
wolffd@0: 
wolffd@0: % OK, just to check things, let's also reconstruct the original FB
wolffd@0: % output.  We do this by multiplying the cepstral data by the transpose
wolffd@0: % of the original DCT matrix.  This all works because we were careful to
wolffd@0: % scale the DCT matrix so it was orthonormal.
wolffd@0: if 1 & (nargout > 3) 
wolffd@0: 	fbrecon = mfccDCTMatrix(1:cepstralCoefficients,:)' * ceps;
wolffd@0: %	imagesc(mt(:,1:cepstralCoefficients)*mfccDCTMatrix);
wolffd@0: end;
wolffd@0: 
wolffd@0: