function [e,w] = per_rls(x,d,N,S,w0,lambda,epsilon)
%per_rls(x,d,N,w0,lambda,epsilon) implements periodic-partial-update RLS algorithm
%   ----------------
%   input parameters
%   ----------------
%   x : Lx1 input signal
%   d : Lx1 desired response 
%   N : filter length
%   S : period of updates (coefficients updated every Sth iteration)
%   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+S,N);
e = zeros(L,1);

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

for i = 1:L-1
    xvec = [x(i);xvec(1:N-1)];
    e(i) = d(i)-w(i,:)*xvec;    %error
    if mod(i-1,S) == 0
        Px = invlambda*P*xvec;      %P(k) computation
        P = invlambda*P - Px*Px'/(1+Px'*xvec);
        w(i+S,:) = w(i,:) + (P*xvec*e(i))';
        w(i:i+S-1,:) = ones(S,1)*w(i,:);
    end
end
xvec = [x(L);xvec(1:N-1)];
e(L) = d(L)-w(L,:)*xvec;
w = w(1:L,:);