%System identification with adaptive filter using partial-update Recursive-Least-Squares 


clear all
close all

mm = menu('Choose an input',...
    'iid Gaussian',...
    'AR(1) Gaussian a=-0.5',...
    'AR(1) Gaussian a=-0.9');

L = 1000;       %input data size
N = 10;         %filter length
M = 1;          %number of coefficients to be updated
nvar = .0001;   %output noise variance   
w0 = zeros(N,1);    %initialization for adaptive filters
lambda = 0.999;     %exponential forgetting factor
epsilon = 1e-4;     %regularization
gamma = 0.02;   %error bound

u = randn(L,1);     %white Gaussian input

if mm == 1,         %input is iid
    x = u;
elseif mm == 2,     %input is AR(1)
    x = filter(1,[1,-.5],u);
elseif mm == 3,     %input is AR(1)
    x = filter(1,[1,-.9],u);
end
    
h = [1,-.8,.3,.2,.1,-.2,.1,.4,-.2,.1];      %system to be identified
n = sqrt(nvar)*randn(L,1);          %output noise
d = filter(h,1,x) + n;              %desired response

%
%Full-update RLS
%
[e,wRLS] = full_rls(x,d,N,w0,lambda,epsilon);
maRLS = sum((wRLS - ones(L,1)*h).^2,2)./norm(h)^2;  %misalignment

%
%Periodic-partial-update RLS
%
S = N;  %period of updates
[e,wPRLS] = per_rls(x,d,N,S,w0,lambda,epsilon); 
maPRLS = sum((wPRLS - ones(L,1)*h).^2,2)./norm(h)^2;  

%
%Sequential-partial-update RLS
%
[e,wSRLS] = seq_rls(x,d,N,M,w0,lambda,epsilon);
maSRLS = sum((wSRLS - ones(L,1)*h).^2,2)./norm(h)^2; 

%
%Stochastic-partial-update RLS
%
[e,wStochRLS] = stoch_rls(x,d,N,M,w0,lambda,epsilon);
maStochRLS = sum((wStochRLS - ones(L,1)*h).^2,2)./norm(h)^2; 

%
%Selective-partial-update RLS
%
[e,wSpuRLS] = spu_rls(x,d,N,M,w0,lambda,epsilon);
maSpuRLS = sum((wSpuRLS - ones(L,1)*h).^2,2)./norm(h)^2; 

%
%Set-membership partial-update RLS
%
[e,wSMRLS,upd] = setmem_rls(x,d,N,M,w0,lambda,epsilon,gamma);
maSMRLS = sum((wSMRLS - ones(L,1)*h).^2,2)./norm(h)^2;

figure(1)
k = 0:L-1;
semilogy(k,maRLS,'k-',k,maPRLS,'b-',k,maSRLS,'r-',k,maStochRLS,'k:',...
    k,maSpuRLS,'b:',k,maSMRLS,'r:');
xlabel('k'); ylabel('Misalignment');
legend('RLS','Per-RLS','Seq-RLS','Stoch-RLS','SPU-RLS','Set-mem par-upd RLS');

figure(2)
stem(k,upd,'filled')
xlabel('k');
ylabel('Set-membership updates');
ylim([0,1.5]);