%LRLS1 Problem 1.1.1.1.2.4
%
%   'ifile.mat' - input file containing:
%      I - members of ensemble
%      K - iterations
%      sigmax - standard deviation of input
%      Wo - coefficient vector of plant
%      sigman - standard deviation of measurement noise
%      epsilon - small auxiliary constant
%      lambda - forgetting factor
%      b - bits in decimal part
% 
%   'ofile.mat' - output file containing:
%      ind - sample indexes 
%      MSE - mean-square error
 
clear all	% clear memory
load ifile;	% read input variables
L=length(Wo);		% plant and filter length
N=L-1;			% plant and filter order
MSE=zeros(K,1);		% prepare to accumulate MSE*I

for i=1:I,		% ensemble
   X=zeros(L,1);    	% initial memory
   x=randn(K,1)*sigmax;		% input 
   n=randn(K,1)*sigman;		% measurement noise 
   V=zeros(L,1);	% inicial output coefficients
   %
   % initial parameters
   delta=zeros(1,L);
   deltaD=zeros(1,L);
   oldgamma=ones(1,L+1);	
   oldeb=zeros(1,L);
   oldxibmin=epsilon*ones(1,L);
   oldxifmin=epsilon*ones(1,L);
   for k=1:K,		% iterations
      X=[x(k)		
         X(1:N)];	% new input vector
      d=Wo'*X;		% desired signal sample
      %
      % initialization
      gamma(1)=1;
      eb(1)=x(k);
      ebp(1)=eb(1);
      ef=x(k);
      xibmin(1)=x(k)^2+lambda*oldxifmin(1);
      xibmin(1)=qround(xibmin(1),b);
      xifmin(1)=xibmin(1);
      e=d+n(k);
      ep=e; 
      for j=1:L,	% orders
	 %
         % auxiliary parameters
         delta(j)=lambda*delta(j)+oldeb(j)*ef/oldgamma(j);
         delta(j)=qround(delta(j),b);
         gamma(j+1)=gamma(j)-(eb(j))^2/xibmin(j);
         gamma(j+1)=qround(gamma(j+1),b);
         %
         % reflection coefficients
	 kappab(j)=delta(j)/xifmin(j);
	 kappaf(j)=delta(j)/oldxibmin(j);
         %
         % prediction errors
	 eb(j+1)=oldeb(j)-kappab(j)*ef;
         eb(j+1)=qround(eb(j+1),b);
         ebp(j+1)=eb(j+1)/gamma(j+1);
	 ef=ef-kappaf(j)*oldeb(j);
         ef=qround(ef,b);
         %
         % minimum least-squares solutions
	 aux=(delta(j))^2;
	 xibmin(j+1)=oldxibmin(j)-aux/xifmin(j);
         xibmin(j+1)=qround(xibmin(j+1),b);
	 xifmin(j+1)=xifmin(j)-aux/oldxibmin(j);
         xifmin(j+1)=qround(xifmin(j+1),b);
	 %
         % feedforward filtering
         ep=ep-V(j)*ebp(j);		
         ep=qround(ep,b);
         deltaD(j)=lambda*deltaD(j)+e*eb(j)/gamma(j);
	 deltaD(j)=qround(deltaD(j),b);
         V(j)=deltaD(j)/xibmin(j);
         e=e-V(j)*eb(j);		
         e=qround(e,b);
      end
      MSE(k)=MSE(k)+ep^2;	% accumulate MSE*I
      %
      % updated parameters
      oldgamma=gamma;
      oldeb=eb;
      oldxibmin=xibmin;
      oldxifmin=xifmin;
   end
end

ind=0:(K-1);		% sample indexes
MSE=MSE/I;		% calculate misadjustment
save ofile ind MSE;	% write output variables