function [h,k,kh,Ph]=levinson(r,d,Py,M)
%
% Programmed by: Dimitris Manolakis, 1998
%
%-----------------------------------------------------------
% 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.
%-----------------------------------------------------------



a=zeros(M-1,1); b=zeros(M-1,1);
k=zeros(M-1,1); kh=zeros(M,1);
h=zeros(M,1);   Ph=zeros(M,1);

% Initialization

P=r(1); beta=r(2);
k(1)=-beta/P;
a(1)=k(1);  b(1)=a(1);
P=P+beta*K(1);
betah=-d(1);
c(1)=d(1)/r(1);
kh(1)=h(1);
Ph(1)=Py+betah*kh(1);

% Recursion

for m=2:M
  betah=-b(1:m-1)'*d(1:m-1)-d(m);
  kh(m)=-betah/P;
  h(1:m)=[(h(1:m-1))' 0]'+[b(1:m-1)' 1]'*kh(m);
  Ph(m)=Ph(m-1)+betah*kh(m);
  if m<M
  beta=(r(2:m))'*flipud(a(1:m-1))+r(m+1);
  k(m)=-beta/P;
  a(1:m)=[(a(1:m-1))' 0]'+[b(1:m-1)' 1]'*k(m);
  b(1:m)=flipud(a(1:m));
  P=P+beta*k(m);
  else
  end
end