function [a, g, e, k] = lpcfitdurb(x,p,h,w)
% [a,g,e,k] = lpcfitdurb(x,p,h,w) Fit short-time LPC model using Durbin recursion
%     x is a signal.  Fit pth-order LPC models to every frame of w 
%     points (default 2*h), hopping h points each time (default 128).
%     Each row of a is returned as the LPC polynomial [1 a1 a2 .. ap]
%     and gain for each frame in g.
%     e is returned as the residual excitation.  k returns reflect coefs.
% 2001-02-27 dpwe@ee.columbia.edu

if nargin < 2
  p = 12;
end
if nargin < 3
  h = 128;
end
if nargin < 4
  w = 2*h;
end

if (size(x,2) == 1)
  x = x';  % Convert X from column to row
end

npts = length(x);

nhops = 1+ floor((npts-w)/h);

a = zeros(nhops, p+1);
g = zeros(nhops, 1);
e = zeros(1, npts);
k = zeros(1,p);

% Pre-emphasis
pre = [1 -0.9];
x = filter(pre,1,x);

an = zeros(1,p);

for hop = 1:nhops

  % Extract segment of signal
  xx = x((hop - 1)*h + [1:w]);
  % Apply hanning window
  wxx = xx .* hanning(w)';
  % Form autocorrelation (calculates *way* too many points)
  rxx = xcorr(wxx,wxx);
  % extract just the points we need (middle p+1 points)
  r = rxx(w+[0:p])';
  
  % Do the Levinson-Durbin recursion
  E = r(1);
  ki = r(2);
  % Unwrap first iteration
  ki = ki/E;
  an(1) = ki;
  E = (1-ki.^2)*E;
  k(hop,1) = ki;
  % Rest of loop
  for i = 2:p
    ki = (r(1+i) - an(1:(i-1)) * r(1+((i-1):-1:1)))/E;
    an(i) = ki;
    an(1:(i-1)) = an(1:(i-1)) - ki*an((i-1):-1:1);
    E = (1-ki.^2)*E;
    k(hop,i) = ki;
  end
  
  % Calculate residual by filtering xx (after windowing)
%  rs = filter([1 -an],1,wxx);
%  G = sqrt(mean(rs.^2));
%  G = sqrt(E);
  % Save them
%  a(hop,:) = [1 -an];
%  g(hop) = G;
%  e((hop - 1)*h + [1:w]) =  e((hop - 1)*h + [1:w]) + rs/G;

  aa = [1 -an];
  rs = filter(aa, 1, xx((w-h)/2 + [1:h]));
  G = sqrt(mean(rs.^2));
  % Save filter, gain and residual
  a(hop,:) = aa;
  g(hop) = G;
  %e((hop - 1)*h + [1:w]) =  e((hop - 1)*h + [1:w]) + rs'/G;
  e((hop - 1)*h + [1:h]) = rs'/G;
  
end