%RLS2 Problem 1.1.1.1.2.3
%
%   '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
%      lambda - forgetting factor
%      b - bits in decimal part
% 
%   'ofile.mat' - output file containing:
%      ind - sample indexes 
%      MSE - mean-square error
%      MSNDW - mean-square norm of coefficient-error vector

clear all	% clear memory
load ifile;	% read input variables
L=length(Wo);		% plant and filter length
N=L-1;			% plant and filter order
ilambda=1/lambda;	% auxiliary constant
delta=1/sigmax^2;	% initialization factor for SD
MSE=zeros(K,1);		% prepare to accumulate MSE*I
MSNDW=zeros(K,1);	% prepare to accumulate MNSDW*I

for i=1:I,		% ensemble
   X=zeros(L,1);    	% initial memory
   SD=delta*eye(L);	
	% initial inverse of deterministic correlation matrix
   W=zeros(L,1);	% initial coefficient vector
   x=randn(K,1)*sigmax;		% input 
   n=randn(K,1)*sigman;		% measurement noise 
   for k=1:K,		% iterations
      X=[x(k)		
         X(1:N)];	% new input vector
      d=Wo'*X;		% desired signal sample
      ep=d+n(k)-W'*X;	
      ep=qround(ep,b);	% a priori error sample
      psi=SD*X;
      psi=qround(psi,b);
      SD=ilambda*qround((SD-psi*psi'/qround((lambda+psi'*X),b)),b);
	% new inverse of deterministic correlation matrix
      W=W+ep*SD*X;		
      W=qround(W,b);	% new coefficient vector
      MSE(k)=MSE(k)+ep^2;	% accumulate MSE*I
      MSNDW(k)=MSNDW(k)+norm((Wo-W),2)^2; 
                        % accumulate MSNDW*I
   end
end

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