rungsdlms.m

System identification with adaptive filter using full and partial-update Generalised-Sideband-Decomp

%System identification with adaptive filter using full and partial-update Generalised-Sideband-Decomposition Least-Mean-Squares 

clear all
close all

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

Ldat = 3000;    %input data length
Ns = 5;         %transform size
K = 2;          %adaptive non-zero coefficients in each subfilter
L = Ns;         %sparsity factor
M = 1;          %number of coefficients to be updated
nvar = .0001;   %output noise variance 
w0 = zeros(K*Ns,1);     %initialization for adaptive filters
mu = 2;         %step-size parameter
lambda = 0.99;  %exponential forgetting factor
epsilon = 50;   %initialization for power estimate
MCruns = 50;    %Monte Carlo simulation runs

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
    
h = [1, -0.8, 0.3, 0.2, 0.1, -0.2, 0.1, 0.4, -0.2, 0.1];    %system to be identified

mse_full = zeros(Ldat,1);      %MSE vectors for Monte Carlo simulations
mse_seq = zeros(Ldat,1);
mse_Mmax = zeros(Ldat,1);

for I = 1:MCruns
    x = filter(1,[1,a],randn(Ldat,1));     %Gaussian input
    n = sqrt(nvar)*randn(Ldat,1);          %output noise
    d = filter(h,1,x) + n;              %desired response
    
    %
    %Full-update GSD-LMS
    %
    e_full = full_gsdlms(x,d,Ns,K,L,w0,mu,lambda,epsilon);
    mse_full = mse_full + e_full.^2;
    
    %
    %Sequential-partial-update GSD-LMS
    %
    e_seq = seq_gsdlms(x,d,Ns,K,L,M,w0,mu,lambda,epsilon);
    mse_seq = mse_seq + e_seq.^2;
    
    %
    %M-max GSD-LMS
    %
    e_Mmax = Mmax_gsdlms(x,d,Ns,K,L,M,w0,mu,lambda,epsilon);
    mse_Mmax = mse_Mmax + e_Mmax.^2;
end

mse_full = mse_full/MCruns;
mse_seq = mse_seq/MCruns;
mse_Mmax = mse_Mmax/MCruns;


%MSE plot
figure
semilogy(1:Ldat,mse_full,'k-',...
    1:Ldat,mse_seq,'b:',...
    1:Ldat,mse_Mmax,'r-.');
%ylim([1e-4,10]);
xlabel('k'), ylabel('MSE')
legend('GSD-LMS','Sequential-partial-update GSD-LMS',...
    'M-max GSD-LMS','Location','NorthEast');