 =========================================================================
%                         CR-RLS Algorithm realization
% =========================================================================
%  For details please see the paper  
% “RLS Algorithm with Convex Regularization”
% IEEE Signal Process. Lett., vol. 18, pp. 470-473, Aug. 2011,
% written by E. M. Eksioglu and A. Korhan Tanc, 
% Realizes the l1-RLS (Auto) curve for Fig.1
% Copyright 2011 Ender M. Eksioglu & A. Korhan Tanc
% =========================================================================

clear all;
close all;

L=64;
S=4;
SNR=10^(20/10);

ff=0.995;
rf=1;

NR=2000;
N=1000+L-1;

MSD=zeros(N,1);

for Nr=1:NR
    
    x=randn(N,1);
    
    h1=sqrt(1/S)*[randn(S,1); zeros(L-S,1)];
    
    nk=randperm(L);

    h=[];

    for n=1:L
    
        h(n)=h1(nk(n));
    
    end
    
    h=transpose(h);
    
    y=conv(x,h)+sqrt(1/SNR)*randn(N+L-1,1);
    
    W=zeros(L,1);
    
    WP=zeros(L,1);

    K=[];

    P=1/rf*eye(L,L);

    for n=L:N
    
        X=x(n:-1:n-L+1);
        
        ye(n)=W'*X;
        
        yp(n)=WP'*X;
        
        e(n)=y(n)-ye(n);
        
        ep(n)=y(n)-yp(n);
        
        K=P*X/(ff+X'*P*X);
        
        gm(n)=max(0,2*( trace(P)/L*(norm(W,1)-norm(h,1)) + sign(W')*P*(WP-W) ) / (sign(W')*P*P*sign(W)) );
                    
        P=(1/ff*P)-(1/ff*K*X'*P);
                   
        W=W+K*conj(e(n))+gm(n)*(ff-1)*P*sign(W);
        
        WP=WP+K*conj(ep(n));
        
        msd(n)=(h-W)'*(h-W);
        
    end
    
    MSD=MSD+transpose(msd);
    
    [Nr, MSD(N)/Nr]
    
end

semilogy(MSD(L:N)/NR);
hold
semilogy(1:100:1000,MSD(L:100:N)/NR,'.');

MSD(N)/NR
