%System identification with adaptive filter using full and partial-update Affine Projection Algorithm 



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
nvar = .001;    %output noise variance   
w0 = zeros(N,1);    %initialization for adaptive filters
mu = 0.1;       %step-size parameter
epsilon = 1e-5; %regularization parameter
gamma1 = 0.08;  %error bound
gamma2 = 0.03;
MCruns = 100;   %Monte Carlo simulation runs

h = [1, -1.1, 0.3, 0.2, 0.1, -0.2, 0.1, 0.05, -0.01, 0.02];      %system to be identified

if mm == 1,         %input is iid
    a = 0;
elseif mm == 2,     %input is AR(1)
    a = -0.5;
elseif mm == 3,     %input is AR(1)
    a = -0.9;
end

mse_apa = zeros(L,1);  %MSE vectors for Monte Carlo simulations
mse_seqapa = zeros(L,1);
mse_Mmaxapa1 = zeros(L,1);
mse_Mmaxapa2 = zeros(L,1);
mse_spuapa = zeros(L,1);
mse_smapa1 = zeros(L,1);
mse_smapa2 = zeros(L,1);

for I = 1:MCruns
    x = filter(1,[1,a],randn(L,1));    %Gaussian input
    n = sqrt(nvar)*randn(L,1);          %output noise
    d = filter(h,1,x) + n;              %desired response

    %
    %Full-update APA
    %
    P = 5;      %APA order
    eAPA = full_apa(x,d,N,P,w0,mu,epsilon);
    mse_apa = mse_apa + eAPA.^2;
    
    %
    %Sequential-partial-update APA
    %
    P = 5;      %APA order
    M = 1;      %no of coefficients to be updated at each iteration
    eSAPA = seq_apa(x,d,N,M,P,w0,mu,epsilon);
    mse_seqapa = mse_seqapa + eSAPA.^2;
    
    %
    %M-max APA (simplified)
    %    
    P = 5;      %APA order
    M = 1;      %no of coefficients to be updated at each iteration
    eMmaxapa1 = Mmax_apa(x,d,N,M,P,w0,mu,epsilon,0);
    mse_Mmaxapa1 = mse_Mmaxapa1 + eMmaxapa1.^2;
 
    %
    %M-max APA (optimal)
    %    
    P = 5;      %APA order
    M=1;        %no of coefficients to be updated at each iteration
    eMmaxapa2 = Mmax_apa(x,d,N,M,P,w0,mu,epsilon,1);
    mse_Mmaxapa2 = mse_Mmaxapa2 + eMmaxapa2.^2;
     
    %
    %Selective-partial-update APA (NB P<=M)
    %
    P = 2;
    M = 8;
    espuAPA = spu_apa(x,d,N,M,P,w0,mu,epsilon);
    mse_spuapa = mse_spuapa + espuAPA.^2;
      
    %
    %Set-membership partial-update APA-1
    %
    P = 2;
    M = 8;
    esmAPA1 = setmem_apa1(x,d,N,M,P,w0,gamma1,epsilon);
    mse_smapa1 = mse_smapa1 + esmAPA1.^2;
    
    %
    %Set-membership partial-update APA-2
    %
    P = 2;
    M = 8;
    esmAPA2 = setmem_apa2(x,d,N,M,P,w0,gamma2,epsilon);
    mse_smapa2 = mse_smapa2 + esmAPA2.^2;
end


mse_apa = mse_apa/MCruns;
mse_seqapa = mse_seqapa/MCruns;
mse_Mmaxapa1 = mse_Mmaxapa1/MCruns;
mse_Mmaxapa2 = mse_Mmaxapa2/MCruns;
mse_spuapa = mse_spuapa/MCruns;
mse_smapa1 = mse_smapa1/MCruns;
mse_smapa2 = mse_smapa2/MCruns;

%MSE plot
figure
semilogy(1:L,mse_apa,'k-',...
    1:L,mse_seqapa,'b:',...
    1:L,mse_Mmaxapa1,'r-.',...
    1:L,mse_Mmaxapa2,'g-',...
    1:L,mse_spuapa,'c-.',...
    1:L,mse_smapa1,'m--',...
    1:L,mse_smapa2,'y:');
ylim([1e-4,10]);
xlabel('k'), ylabel('MSE')
legend('APA','Sequential-partial-update APA','M-max APA (simplified)',...
    'M-max APA (optimal)','Sel-par-upd APA',...
    'Set-mem par-upd APA-1','Set-mem par-upd APA-2','Location','NorthEast');