function [a,g,e] = lpcfit(x,p,h,w,ov)% [a,g,e] = lpcfit(x,p,h,w,ov)  Fit LPC to short-time segments%    x is a stretch of signal.  Using w point (2*h) windows every %    h points (128), fit order p LPC models.  Return the successive %    all-pole coefficients as rows of a, the per-frame gains in g %    and the residual excitation in e.%    ov nonzero selects overlap-add of window-length%    residuals, otherwise successive hop-sized residuals are concatenated%    for independent near-perfect reconstruction with lpcsynth.%    (default is 1)if nargin < 2  p = 12;endif nargin < 3  h = 128;endif nargin < 4  w = 2*h;endif nargin < 5  ov = 1;endif (size(x,2) == 1)  x = x';  % Convert X from column to rowendnpts = length(x);nhops = floor(npts/h);% Pad x with zeros so that we can extract complete w-length windows% from itx = [zeros(1,(w-h)/2),x,zeros(1,(w-h/2))];a = zeros(nhops, p+1);g = zeros(nhops, 1);if ov == 0  e = zeros(1, npts);else  e = zeros(1, (nhops-1)*h+w);end% Pre-emphasispre = [1 -0.9];x = filter(pre,1,x);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);  % extract just the points we need (middle p+1 points)  rxx = rxx(w+[0:p]);  % Setup the normal equations  R = toeplitz(rxx(1:p));  % Solve for a (horribly inefficient to use full inv())  an = inv(R)*rxx(2:(p+1))';  % Calculate residual by filtering windowed xx  aa = [1 -an'];  if ov == 0    rs = filter(aa, 1, xx((w-h)/2 + [1:h]));  else    rs = filter(aa,1,wxx);  end  G = sqrt(mean(rs.^2));  % Save filter, gain and residual  a(hop,:) = aa;  g(hop) = G;  if ov == 0    e((hop - 1)*h + [1:h]) = rs'/G;  else    e((hop - 1)*h + [1:w]) =  e((hop - 1)*h + [1:w]) + rs/G;  endend% Throw away first (win-hop)/2 pts if in overlap mode% for proper synchronization of resynthif ov ~= 0  e = e((1+((w-h)/2)):end);end