%SFRLS1 Problem 1.1.1.1.2.7
%
%   '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
%      kappa1, kappa2, kappa3 - auxiliary parameters
%      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
L1=L+1;			% auxiliary constant
MSE=zeros(K,1);		% prepare to accumulate MSE*I
MSNDW=zeros(K,1);	% prepare to accumulate MNSDW*I

for i=1:I,		% ensemble
   Xe=zeros(L1,1);	% extended initial memory
   x=randn(K,1)*sigmax;		% input 
   n=randn(K,1)*sigman;		% measurement noise 
   %
   % initial coefficient vectors
   Wb=zeros(L,1);	
   Wf=Wb;	
   W=Wb;	
   %
   % initialization
   phih=zeros(L,1);
   gamma=1;
   xibmin=epsilon;
   ixifmin=1/epsilon;
   for k=1:K,		% iterations
      Xe=[x(k)
          Xe(1:L)];	% new extended input vector
      %
      % predictions
      efp=Xe'*[1 
               -Wf];
      efp=qround(efp,b);
      ef=efp*gamma;
      phihe=[0
             phih]+ixifmin/lambda*[1
                                  -Wf]*efp;
      phihe=qround(phihe,b);
      gamma=1/qround((1/gamma+phihe(1)*efp),b);	
			% gamma1
      ixifmin=ixifmin/lambda-gamma*(phihe(1))^2;
      ixifmin=qround(ixifmin,b);
      Wf=Wf+phih*ef;
      Wf=qround(Wf,b);
      ebp1=lambda*xibmin*phihe(L1);    
      ebp2=[-Wb' 1]*Xe;     
      ebp2=qround(ebp2,b);
      ebp=[1-kappa1 kappa1
           1-kappa2 kappa2
           1-kappa3 kappa3]*[ebp1 
                             ebp2];
      ebp=qround(ebp,b);
      gamma=1/qround((1/gamma-phihe(L1)*ebp(3)),b);
			% gamma2
      eb=gamma*ebp;
      xibmin=lambda*xibmin+eb(2)*ebp(2);
      xibmin=qround(xibmin,b);
      phih=phihe(1:L)+phihe(L1)*Wb;
      phih=qround(phih,b);
      Wb=Wb+phih*eb(1);
      Wb=qround(Wb,b);
      gamma=1/qround((1+phih'*Xe(1:L)),b);
			% gamma3
      %
      % joint-process estimation
      d=Wo'*Xe(1:L);		% desired signal sample
      ep=d+n(k)-W'*Xe(1:L);
      ep=qround(ep,b);
      e=ep*gamma;		% error sample
      W=W+phih*e;		
      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