ica_6.m

盲信号处理牛人cardoso

echo off; 
clear all;  
close all;
format long;

NT=2;
NG=20000;
L=20;
NN=NG+NT;
N=6;     
M=6;     
T=0.0001; 
NNN=T*(1:NN);
snr1=sqrt(0.);
arfa=0.995;
  
% five deterministic sub-Gaussian signals
s1=sign(cos(2*pi*155*NNN));  
s2=sin(2*pi*800*NNN);
s3=sin(2*pi*90*NNN);
s4=sin(2*pi*10*NNN).*sin(2*pi*300*NNN);
s5=sin(2*pi*300*NNN+6*cos(2*pi*60*NNN)); 
% one random over-Gaussian signal 
ar1=5/20;
hh=zeros(L,1);
for n1=1:L
    a1=n1*ar1;
    hh(n1)=randn*exp(-a1);
end    
NN1=NN+L-1
nh=1-2*rand(1,NN1); 
for n1=1:NN1       
    if nh(n1)>0
        nh(n1)=-log(nh(n1));
    end
    if nh(n1)<0        
        nh(n1)=log(-nh(n1));
    end
end
for n1=1:NN
    nj(n1)=nh(n1:n1+L-1)*hh;
end
nj=1-2*rand(1,NN); 
S=[s1;s2;s3;s4;s5;nj]; 
clear s1; clear s2; clear s3; clear s4; clear s5; clear n;
% 100 independent runs are executed
for kkk=1:1;   kkk,    
    % A is randomly generated in each run
    A=(1-2*rand(M,N));       
    X=A*S;
    Xn=randn(M,NN);
    aa1=snr1*norm(X,'fro')/norm(Xn,'fro');
    Xn=aa1*Xn;
    X=X+Xn;        
    W1=eye(M);
    W2=W1;
    W3=W1;
    W4=W1;
    Sest1=zeros(M,NN);
    R=zeros(M,M);
    arpha=0.99;

    for i1=1:NG          
        % the Yong algorithm A' 
        if i1<=600   eta1=120*T; 
        else eta1=120*T*exp(-15*T*(i1-600));
        end 
        eta1=0.005;
        if i1<=600   MU=60; 
        else MU=60*exp(-15*T*(i1-600));
        end 
        MU=60;
        %------------ Natural Gradient -------------------- 
        eta6=0.001;
        x1=X(:,i1); 
        x2=X(:,i1+NT);      
        Y1=W1*x1;         
        Sest1(:,i1)=Y1;
        if  i1==1;  
            K3=Y1.^3;
            K4=Y1.^4-3;
        else
            K3=K3-MU*T*(K3-Y1.^3);
            K4=K4-MU*T*(K4-Y1.^4+3);
        end   
        f=-0.5*K3+2.25*K3.*K4;
        g=-1/6*K4+1.5*K3.^2+0.75*K4.^2;       
      
        mat1=Y1.^3*Y1';              
        W1=W1+3*eta6*(eye(M)-mat1)*W1;
        P1EFF(kkk,i1)=Eff(W1*A);
        %---------------- Relative Gradient-------------------
        eta2=eta1/1.5;
        Y1=W2*x1; 
        Y2=W2*x2;        
        mat1=Y1.^3*Y1'; 
        %mat1=Y2*Y1';                             
        W2=W2+((eye(M)-Y1*Y1')-mat1+mat1')*W2*(eta6);            
        P2EFF(kkk,i1)=Eff(W2*A); 
        %---------------- New 1 -----------------------------
        if i1>1
            a1=1/i1;          
            R=R+(X(:,i1)*X(:,i1)'-R)*a1;
        else
            R=X(:,i1)*X(:,i1)';
        end
        B1=R^(-0.5);         
        z1=B1*x1;
        z2=B1*x2;        
      
        eta3=0.001;        
        Aj=W3*W3';
        %Aj=zeros(N,N);
        Y1=W3*z1; 
        Y2=W3*z2;          
        mat1=Y1.^3*Y1';
        %mat1=Y2*Y1';                           
        W3=W3+( -mat1+mat1')*W3*(eta3);  %   (eye(M)-Aj)
        P3EFF(kkk,i1)=Eff(W3*B1*A);       
        tx(kkk,i1)=i1;
        %---------------- New N -----------------------------
        Y1=W4*z1;
        mat1=Y1.^3*Y1';              
        W4=W4+3*eta6*(eye(M)-mat1)*W4;
        P4EFF(kkk,i1)=Eff(W4*B1*A);
        
   end    

   figure;  grid off;    hold on;    
   plot(tx,10*log10(P1EFF(kkk,:)),'b',tx,10*log10(P2EFF(kkk,:)),'r',tx,10*log10(P3EFF(kkk,:)),'g',tx,10*log10(P4EFF(kkk,:)),'y');
   
end