function [k,var] = aplatest(x,P)
% Estimation of All-Pole lattice parameters using LS Forward/Backward
%  error minimization using Burg approach
% [k,var] = aplatest(x,P)
%  Inputs:
%   x : signal data values 
%   P : model-order
% Outputs:
%   k : lattice (Burg) parameters for m = 1,2,...,P
% var : variance of input white noise
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------

% Function Initialization
N = length(x); ef = zeros(P,N); eb = ef;
Ef = zeros(P); Eb = Ef; Bfb = Ef;
kFP = zeros(P); kBP = kFP; kB = kFP;

% Input data                            % Step-1
x = reshape(x,N,1);

% Algorithm Initialization              % Step-2
ef0 = x; eb0 = x; Ni = 0+1; Nf = N-1+1; % Step-2a
Ef0 = real(ef0(Ni:Nf)'*ef0(Ni:Nf));     % Step-2b: forward error
Eb0 = real(eb0(Ni:Nf)'*eb0(Ni:Nf));     % Step-2b: backward error
Bfb0 = ef0(Ni+1:Nf)'*eb0(Ni:Nf-1);      % Step-2b: f/b error
kFP0 = -Bfb0/Eb0; kBP0 = -Bfb0/Ef0;     % Step-2c
kB0 = 2*kFP0*kBP0/(kFP0+kBP0);          % Step-2d
ef(1,Ni) = ef0(Ni);                     % Step-2e: Begin error computation
eb(1,Ni) = kB0*ef0(Ni);
for n = Ni+1:Nf
   ef(1,n) = ef0(n) + conj(kB0)*eb0(n-1);
   eb(1,n) = eb0(n-1) + kB0*ef0(n-1);
end                                     % Step-2e: End error computation

% Algorithm Recursion                   % Step-3
for m = 2:P
   Ef(m-1)  = real(ef(m-1,m:end)*ef(m-1,m:end)');
   Eb(m-1)  = real(eb(m-1,m:end)*eb(m-1,m:end)');
   Bfb(m-1) = ef(m-1,m:end)*eb(m-1,m-1:end-1)';
   kFP(m-1) = -Bfb(m-1)/Eb(m-1); 
   kBP(m-1) = -Bfb(m-1)/Ef(m-1);
   kB(m-1)  = 2*kFP(m-1)*kBP(m-1)/(kFP(m-1)+kBP(m-1));
   for n = m:N
      ef(m,n) = ef(m-1,n) + conj(kB(m-1))*eb(m-1,n-1);
      eb(m,n) = eb(m-1,n-1) + kB(m-1)*ef(m-1,n-1);
   end
end

% Output                                % Step-4
k = [kB0,kB(1:P-1)];
var = x'*x/N;
for m = 1:P; var = var*(1-k(m)*conj(k(m))); end;