mmax_gsdlms.m

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

function [e,w] = Mmax_gsdlms(x,d,Ns,K,L,M,w0,mu,lambda,epsilon)
%[e,w] = Mmax_gsdlms(x,d,Ns,K,L,M,w0,mu,lambda,epsilon) implements the
%M-max GSD-LMS algorithm using DCT and single-division power normalization.
%   ----------------
%   input parameters
%   ----------------
%   x : Ld x 1 input signal
%   d : Ld x 1 desired response 
%   Ns : DCT size
%   K : number of non-zero adaptive coefficients in each sparse subfilter
%   L : sparsity factor for subfilters
%   M : number of coefficients to be updated
%   w0 : Nx1 initialization
%   mu : step-size parameter
%   lambda : exponential forgetting factor for transform output power
%   estimation 
%   epsilon : initialization for power estimates
%   ----------------
%   function outputs
%   ----------------
%   w : Ld x KNs evolution of coefficients
%   e : Ld x 1 output error vector

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

Ld = length(x);
w = zeros(Ld,K*Ns);
e = zeros(Ld,1);

w(1,:) = w0';
xvec = zeros(Ns,1);

rpk = (1/epsilon)*ones(Ns,1);    %initial inv power estimate for single-division power normalization

invlambda = 1/lambda;

vnorm = zeros(Ld,Ns);
vvec = zeros(Ld,Ns);
mmax = zeros(Ld,Ns);

for i = 1:Ld-1
    xvec = [x(i);xvec(1:Ns-1)];
    vvec(i,:) = dct(xvec)';
    V = vvec(i:-L:max([i-(K-1)*L,1]),:)';
    V = V(:);
    V = [V;zeros(K*Ns-length(V),1)];
    e(i) = d(i) - w(i,:)*V;
    
    %single-division power normalization
    bk = invlambda*rpk;
    ck = bk.*vvec(i,:)';
    vnorm(i,:) = ck';        %delayed normalized v(k)
    ck = ck.* vvec(i,:)';
    dk = bk.*ck;
    mmax(i,:) = ck';        %entries to be sorted
    rpk = bk - dk/(1+sum(ck));      %update reciprocal power estimate
    
    vv = vnorm(i:-L:max([i-(K-1)*L,1]),:)';  %delayed normalized v(k), v(k-L), ..., v(k-(K-1)L)
    vv = vv(:);
    vv = [vv;zeros(K*Ns-length(vv),1)];
    upd = mu * e(i) * vv;           %delayed normalized v_a(k)
   
    mm = mmax(i:-L:max([i-(K-1)*L,1]),:)'; 
    mm = mm(:);
    mm = [mm;zeros(K*Ns-length(mm),1)];
    
    [xs,ix] = sort(-mm);     %rank
    de = zeros(1,K*Ns);
    de(ix(1:M)) = 1;
    IM = diag(de);

    upd = IM*upd;           %Mmax selection
    
    w(i+1,:) = w(i,:) + upd';
end

xvec = [x(Ld);xvec(1:Ns-1)];
vvec(Ld,:) = dct(xvec);
V = vvec(Ld:-L:max([Ld-(K-1)*L,1]),:)';
V = V(:);
V = [V;zeros(K*Ns-length(V),1)];
e(Ld) = d(Ld) - w(Ld,:)*V;