function [e,w] = seq_rls(x,d,N,M,w0,lambda,epsilon)
%seq_rls(x,d,N,M,w0,lambda,epsilon) implements sequential-partial-update
%RLS algorithm 
%   ----------------
%   input parameters
%   ----------------
%   x : Lx1 input signal
%   d : Lx1 desired response 
%   N : filter length
%   M : number of coefficients to be updated
%   w0 : Nx1 initialization
%   lambda : exponential forgetting factor
%   epsilon: regularization parameter 
%   ----------------
%   function outputs
%   ----------------
%   e : Lx1 output error vector
%   w : LxN coefficients vectors

x = x(:);
d = d(:);
w0 = w0(:);

L = length(x);
w = zeros(L,N);
e = zeros(L,1);

w(1,:) = w0';
xvec = zeros(N,1);
P = (1/epsilon) * eye(N);
invlambda = 1/lambda;

B = N/M;
if B == floor(B),
    for i = 1:B
        J(i,1:M) = (i-1)*M+1:i*M;
    end
else
    B = ceil(B);
    for i = 1:B-1
        J(i,1:M) = (i-1)*M+1:i*M;
    end
    J(B,1:N-(B-1)*M) = (B-1)*M+1:N;
end

for i = 1:L-1
    xvec = [x(i);xvec(1:N-1)];
    e(i) = d(i)-w(i,:)*xvec;    %error
    
    ix = mod(i-1,B)+1;
    de = zeros(1,N);
    de(J(ix,:)) = 1;
    IM = diag(de);
    
    %PM(k) computation
    aa = 1/(1+invlambda*xvec'*P*IM*xvec);
    P = invlambda*P - (invlambda^2*P*IM*xvec*xvec'*P)*aa;
    
    upd = P*IM*xvec*e(i);
    w(i+1,:) = w(i,:) + upd';  
end
xvec = [x(L);xvec(1:N-1)];
e(L) = d(L)-w(L,:)*xvec;