MFCC

function [cepstra, aspectrum, pspectrum] = melfcc(samples, sr, wintime, steptime, nfilts, numcep, preemph)
% [cepstra, aspectrum, pspctrum] = melfcc(samples, sr, wintime, steptime, numcep, preemph)
%   - take power spectra of the STFT
%   - warp to a mel frequency scale
%   - take the DCT of the log-Mel-spectrum
%   - return the first <numcep> components
%   samples:                        vector of signal
%   sr:                                  sample rate
%   wintime (0.025):             window length in second
%   steptime (0.010):           step between successive windows in second
%   numcep (13):                 number of cepstra to return
%   nfilts (40):                       number of triangle filter to use
%   preemph (0.97):             pre-emphasis filter coefficient

if nargin < 2;      sr = 16000;               end
if nargin < 3;      wintime = 0.025;       end
if nargin < 4;      steptime = 0.010;     end
if nargin < 5;      numcep = 13;           end
if nargin < 6;      nfilts = 40;                end
if nargin < 7;      preemph = 0.97;       end

winpts = round(wintime*sr);
steppts = round(steptime*sr);
NFFT = 2^(ceil(log2(winpts)));

samples = filter([1, -preemph], 1, samples);

%   compute FFT power spectrum
%pspectrum = powspec(samples, sr, wintime, steptime, NFFT);
%pspectrum = abs(spectrogram(samples*32768, winpts, winpts - round(steptime*sr), NFFT)).^2;

[frame, tc] = enframe(samples, hamming(winpts), steppts);
pspectrum = abs(fft(frame', NFFT)).^2;

%   obtain mel filter bank
filterBank = melFilterBank(NFFT, sr, nfilts);

%   auditory spectrum
aspectrum = filterBank * pspectrum;

%   apply DCT to convert aspectrum to mel cepstrum
logAspec = log(aspectrum);
cepstra = dct(logAspec);
cepstra = cepstra(1:numcep,:);

if nargout < 1
    [nf, nc] = size(cepstra);
    imh = imagesc(tc/sr, 1:nf, cepstra);
    axis('xy');
    xlabel('Time (s)');
    ylabel('Mel-cepstrum coefficient');
    map = (0:63)' / 63;
    colormap([map, map, map]);
    colorbar;
end
end