function mf = mfcc(sig, MFCC_or, Fs);
%
% mf = mfcc(sig, MFCC_or, Fs);
% stylianou
% 

if(nargin<3)
     Fs = 16000;
end

% definitions
% Definition of filter bank for mfcc:
LIM_BF = [0 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 ...
          3150 3700 4400 5300 6400 7700 9500 12000 15500];

% plot the filterbank
if 0
	plot(LIM_BF(1:3),[0 1 0],'k'); hold on;
	for m=1:21
		plot(LIM_BF([0:2]+m),[0 1 0],'k');
	end
	hold off;
	return
end

bf_max = 25;
CM = cos((((1:bf_max-1)-1/2)'*(1:MFCC_or))*pi/(bf_max-1));

sig = sig(:);
L = length(sig);
Win = hanning(L);

sig = sig.*Win;
a = sqrt(1/sum(sig.^2)); % normalization
sig = a*sig;

nfft = 1024;

S = fft(sig,nfft);
S = abs(S.^2/(sum(Win)/2)^2);
ndsp = zeros(1,(bf_max-1));
f =  zeros(1,(bf_max-1));
for Bf_k = 1:(bf_max-1)
    ndsp(Bf_k-1+1) = sum( S(1+round(nfft*LIM_BF(Bf_k)/Fs) : 1+round(nfft*LIM_BF(Bf_k+1)/Fs)) );
    % normalize by the number of points accumulated for the computation of this band
    ndsp(Bf_k-1+1) = ndsp(Bf_k-1+1)/ ( (LIM_BF(Bf_k+1)-LIM_BF(Bf_k)) / Fs*nfft );
    f(Bf_k-1+1) = ((LIM_BF(Bf_k+1) + LIM_BF(Bf_k))/2)/Fs;
end

% mel-scaled frequencies cepstrum coefficients
mf = log10(ndsp)*CM;