%***********************************************************
%   This mfile implements the exponentially-weighted 
%   RLS algorithm (Hayes 1996)
%   AUTHOR  :  Santhanam  ,   Balu
%   DATE    :  11/06/97
%
%   USAGE
% 
%   [W,y,e] = rls(x,d,lambda,ord,delta);
%   SYNOPSIS: 
%   W      : (Nxord) (Adaptive Filter Coefficients)
%   y      : Filtered output signal
%   x      : observation signal
%   d      : desired signal
%   lambda : Memory factor (0 - 1) 
%   delta  : Positive constant for P[n] initial value
%************************************************************
function [W,y,e] = rls(x,d,lambda,ord,delta)
format long
x = x(:); d = d(:);
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
for k = 2:1:M-ord+1
    z = P*X(k,:)';
%   
%   RLS  gain  vector computation: 
    g = z/(lambda + X(k,:)*z);
%  
%   Innovations process: 
    alpha = d(k) - X(k,:)*W(k-1,:).';
%
%   Tap weight update: 
    W(k,:) =  W(k-1,:) + (alpha*g.');
%   
%   Filtered output signal
    y(k) = X(k,:)*W(k,:).';
%   
%   Estimation error:
    e(k) = d(k) - y(k);
%
%   Inverse matrix P[n] update: 
    P = (P - g*z.') / lambda;
end
%********************************************************  

