function [a,g]=lpc(h,N);
%LPC  Linear Predictive Coefficients.
%   [A,G] = LPC(X,N) finds the coefficients and gain of an Nth order 
%   auto-regressive process that models the time series X as follows:
%      X(n) = G*( -A(2)*X(n-1) - A(3)*X(n-2) - ... - A(N+1)*X(n-N) )
%   Here, X is the real input time series (a vector), and N is the order of 
%   denominator polynomial A(z), i.e. A = [ 1 A(2) ... A(N+1) ].
%
%   If you leave N unspecified, LPC uses a default N = LENGTH(X)-1.
%
%   LPC uses the autocorrelation method, also known as the maximum
%   entropy method (MEM) of spectral estimation. It computes a biased
%   estimate of the autocorrelation function of X, and then calls the
%   LEVINSON function.
%
%   If X is a matrix containing a separate signal in each column, LPC 
%   returns a model estimate for each column in the rows of A, and a row 
%   vector of gains G.
%
%   See also PRONY, STMCB, LEVINSON.

%   R is the auto correlation vector, R(1) = E(h(t)h*(t)), 
%       R(2) = E(h(t+1)h*(t)),...
%   This method is sometimes called the MEM method of spectral estimation.

%   Author(s): T. Krauss, 9-21-93
%   Copyright (c) 1988-97 by The MathWorks, Inc.
%   $Revision: 1.11 $  $Date: 1997/04/14 20:01:27 $

    error(nargchk(1,2,nargin))

    [r,c] = size(h);
    if (c>1)&(r==1)
        h = h(:);
    end
    numsigs = size(h,2);
    if nargin<2, N = size(h,1)-1; end
    if (N>size(h,1)-1),
        % disp('Warning: zero-padding short input sequence')
        h(N+1,:)=zeros(1,numsigs);
    end

    R = flipud(fftfilt(h,conj(flipud(h))))/length(h);
    % R is the autocorrelation vector.  
    % Equivalent code (for the single signal case):
    %   R = xcorr(h,'biased');
    %   M = length(h);
    %   R(1:M-1) = [];

    a = levinson(R,N);
    g = real(sqrt( sum((a.').*R(1:N+1,:)) ));

