function [cps,dct,R]=all_pole(s,LPC_ORDER,CPS_ORDER,icsi_yes);
% ALL_POLE modeling, [cps,dct]=all_pole(s,LPC_ORDER,icsi_yes)
% everything hard-wired inside. Reads mel fb from ~/matlab/fbanks/
% can eventually output also the dct
% 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
N=size(mfb,2);

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

% tests with preemphasis on the spec:
%qq = abs(freqz([1 -0.97],1,129)).^2;  
%S = S .* (qq * ones(1,Nfr)); 

fbout = mfb' * S;       

% cube root
fbout3 = fbout .^ (1/3);
%fbout3 = log(fbout); 

if (icsi_yes == 0)   % NORMAL CASE
  % dct basis in cols - first define basis:
  dct=[]; k=1:23;
  for n=0:22; % was 14 
    bas= 1/N * cos(pi * n / N * (k-0.5)); %plot (bas); pause; 
    dct = [dct bas'];
  end
  
  % multiply by the dct
  R = dct' * fbout3;
  
else %%%%%%%%%%%%% ICSI case 
  % define basis
  dct=[]; k=0:24; N=24;
  for n=0:22; % was 14 
    bas= 1/N *cos(pi * n / N * k); 
    % 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' * [fbout3(1,:); fbout3; fbout3(23,:)];
end

% now all-pole, do also gain acc to definition 
cps = zeros(CPS_ORDER,Nfr); %aux=[];
for ii=1:Nfr;
  r = R(1:LPC_ORDER+1,ii); %figure(1); plot (r); grid
  a = levinson(r); 
%  aux = [aux a(2:LPC_ORDER+1)'];
  c = a_to_cepst (a',CPS_ORDER-1);
  gain2 = r'*a';
  c(1) = log(gain2);
%  c(1) = log(r(1));
  cps(:,ii) = c;
%  pause
end