function [cps,dct]=mfcc(s,CPS_COEFFS,icsi_yes);
% MFCC, [cps,dct]=mfcc(s,CPS_COEFFS,icsi_yes)
% everything hard-wired inside. Reads mel fb from ~/matlab/fbanks/
% can eventually output also the dct
% CPS_COEFFS is the number of cps. coeffs including the 0-th one. 
% if icsi_yes is !=0, will do repetition of 1st and last channels and 
% icsi-style dct 
%
% error affecting the scaling corrected Mon May 21 17:35:06 PDT 2001
% should be 1/N to give same results as the ifft. 

load /u/honza/matlab/fbanks/mfb
%load /home/cernocky/OGI/PLP/mfb
N=size(mfb,2);
CPS_ORDER=15;

S = abs(specgram(s,256,8000,hamming(200),120)).^2;
Nfr = size(S,2);

fbout = mfb' * S;       

% log
fboutlog = log(fbout);

if (icsi_yes == 0)   % NORMAL CASE
  % dct basis in cols - first define basis:
  dct=[]; k=1:23;
  for n=0:(CPS_COEFFS-1);
    bas= 1/N * cos(pi * n / N * (k-0.5)); %plot (bas); pause; 
    dct = [dct bas'];
  end
  
  % multiply by the dct
  R = dct' * fboutlog;
  
else %%%%%%%%%%%%% ICSI case 
  % define basis
  dct=[]; k=0:24; N=24;
  for n=0:(CPS_COEFFS-1);
    bas= 1 / N*cos(pi * n / N * k);  % was 2 / N !!!!!!!!!!!!!!!!!!!!! 
    % but 1st and last should be there just once.
    bas(1) = bas(1) / 2; bas(N+1) = bas(N+1) /2;
    %plot (bas); pause;
    dct = [dct bas'];
  end
  % repeat channel, multiply by the dct
  R = dct' * [fboutlog(1,:); fboutlog; fboutlog(23,:)];
end
  
% nothing more, R is the result ! 
cps = R;