%***********************************************************
%   This mfile implements IF tracking using 
%   the hybrid LMS/RLS algorithm that uses 
%   a steepest descent update for the memory factor 
%   AUTHOR  :  Santhanam  ,   Balu
%   DATE    :  11/06/97
%
%   USAGE
% 
%   [W,y,MSE,lambda] = grls(sig,lambda_in,ord,delta,rho);
%   SYNOPSIS: 
%
%   lambda_in      : Initial memory factor 
%   ord            : order of the predictor
%   delta          : Positive constant for P[n] initial value
%   rho            : Learning rate for the memory factor
%************************************************************
function  [W,y,MSE,lambda] = grls(sig,lambda_in,ord,delta,rho)
format long
x = rshft(sig(:),1); 
d = sig(:);
%-----------Initialization--------------------------------------
X = convmtx(x,ord); % Generate convolution matrix
[M,N] = size(X); % 
P = eye(N)/delta; % Initialize P[n] recursion
W(1,:) = zeros(1,N); % Initialize Tap-weight vector
lambda(1) = lambda_in; % Intitialize the memory factor update
g = zeros(N,N);
z = zeros(N,N);
MSE = zeros(1,N);
%----------------------------------------------------------------
for k = 2:1:M-ord+1
    %   Innovations process: 
    alpha(k) = d(k) - X(k,:)*W(k-1,:).';
%   
%----------------------Memory factor update----------------------
    num = alpha(k)*alpha(k-1)*X(k,:)*z(:,k-1);
    den = (lambda(k-1) + X(k-1,:)*z(:,k-1))^2;
    lambda(k) = lambda(k-1) - rho*num/den;
%   Truncate the forgetting factor update
    if lambda(k) > 1 
       lambda(k) = 0.999;
   elseif lambda(k) < 0
       lambda(k) = 0.001;
   end
%----------------------------------------------------------------
%  Least-squares solution vector
    z(:,k) = P*X(k,:)';
%  

%   RLS  gain  vector computation: 
    g(:,k) = z(:,k)/(lambda(k) + X(k,:)*z(:,k));
%
%   Tap weight update: 
    W(k,:) =  W(k-1,:) + (alpha(k)*g(:,k).');
%   
%   Filtered output signal
    y(k) = X(k,:)*W(k,:).';
%   
%   Inverse P[n] update update: 
    P = (P - P*g(:,k)*z(:,k).') / lambda(k);
%   
%   MSE update 
    MSE(k) = lambda(k)*MSE(k-1) + (d(k) - y(k))*alpha(k);
end
clear z g P num den alpha y X
%-------------------------------------------------------------------