efica.m

ICALAB for Signal Processing Toolbox for BSS, ICA, cICA, ICA-R, SCA and MCA

function [Wefica, ISRef, Hef, Wsymm, ISRsymm, status]= efica(X, ini, FASTICApars, Saddlepars,EFICApars)
% function [Wefica, ISRef, Hef, Wsymm, ISRsymm, status]=efica(X, ini, g, SaddleTest)

%EFICA: [Wefica, ISRef, Hef, Wsymm, SIRsymm, status]=efica(X, ini, g, SaddleTest)
%
% version: 1.9  release: 24.11.2006
%
%          %improved speed
%
%Input data:
% X ... mixed data dim x N, where dim is number of signals and N is number of
%       samples
% ini ... starting point of the iterations for W matrix
% g ... nonlinearity used in the symmetric algorithm ('tanh'-'rati'-'gaus'-'pow3')
% SaddleTest ... if true (default) the test of saddle points on the demixed signals is done
%
%
%Output data:
% W_efica - demixing matrix produced by EFICA
% ISRef - ISR matrix estimator of EFICA components
% Hef - matrix for EFICA bias computation (in the noisy environment)
% Wsymm - demixing matrix produced by symmetric Fast-ICA
% ISRsymm - ISR matrix estimator of Fast-ICA components
% status - one if test of saddle points was positive, zero otherwise
%
% Examples of usage:
%
%  [Wefica, ISRef, Hef, Wsymm, SIRsymm]=efica(X);
%
%    use default settings (random initialization, 'tanh' nonlinearity,
%    the test of saddle points will be done)
%
%
%  [Wefica, ISRef, Hef, Wsymm, SIRsymm]=efica(X, randn(dim), 'rati', true);
%
%    random initialization, 'rati' nonlinearity used in symmetric FastICA,
%    the test of saddle points will be done
%
%
%  [Wefica, ISRef, Hef, Wsymm, SIRsymm]=efica(X, eye(dim), 'tanh', false);
%
%    random initialization, 'tanh' nonlinearity used in symmetric FastICA,
%    without the test of saddle points
%
%
%
%
%

[dim N]=size(X);

%Default values of parameters
if nargin<5 %EFICApars 
%     EFICApars = true;
    EFICAMaxIt=50; %%Maximum number of improving iterations
    eficaGaussianNL=3; %'ggda';  %Nonlinearity used for eficagaussian signals 'gaus'-'ggda'-'npnl'
else
    EFICAMaxIt=EFICApars.EFICAMaxIt ; %%Maximum number of improving iterations
    eficaGaussianNL=EFICApars.eficaGaussianNL;  %Nonlinearity used for eficagaussian signals 'gaus'-'ggda'-'npnl'
end

if nargin<4     %Saddlepars
    SaddleTest = 1;
    test_of_saddle_points_nonln = 4 ; %'rtsp';
    MaxItAfterSaddleTest=30; %%Maximum number of iterations after a saddle point was indicated
else
    SaddleTest = Saddlepars.SaddleTest ;
    test_of_saddle_points_nonln=Saddlepars.test_of_saddle_points_nonln;
    MaxItAfterSaddleTest=Saddlepars.MaxItAfterSaddleTest; %%Maximum number of iterations after a saddle point was indicated

end
if nargin<3 %FASTICApars
    g= 3; %'rati';
    MaxIt=100; %% Maximum number of FastICA iterations
else
    g=FASTICApars.g;
    MaxIt=FASTICApars.MaxIt; %% Maximum number of FastICA iterations

end

if nargin<2
    
    while 1
        ini=randn(dim);
        if cond(ini) < 10^3
            break
        end
    end
end

epsilon=0.0001; %Stop criterion
fineepsilon=1e-5; %Stop criterion for post-estimation
repeat=1;
rot2d=[1/sqrt(2) 1/sqrt(2);-1/sqrt(2) 1/sqrt(2)];
status=0;
min_correlation=0.8; %additive noise...0.75, noise free... 0.95, turn off (unstable)...0

%removing mean
X = X-mean(X,2)*ones(1,N);


%fprintf('EFICA...\n');

%preprocessing
C = cov(X');
CC = C^(-1/2);
Z = CC*X;

%FastICA
W=ini;
W_before_decor=W;
W=W*real(inv(W'*W)^(1/2));
%Wold=zeros(dim,dim);
crit=zeros(1,dim);
NumIt=0;
while repeat
    %%% Symmetric approach
    while (1-min(crit)>epsilon && NumIt<MaxIt)% && sum(double(changed>10))<2)
        Wold=W;
        switch g
            case 1 %'tanh'
                hypTan = tanh(Z'*W);
                W=Z*hypTan/N-ones(dim,1)*sum(1-hypTan.^2).*W/N;
            case 2 %'pow3'
                W=(Z*((Z'*W).^ 3))/N-3*W;
            case 3 %'rati'
                U=Z'*W;
                Usquared=U.^2;
                RR=4./(4+Usquared);
                Rati=U.*RR;
                Rati2=Rati.^2;
                dRati=RR-Rati2/2;
                nu=mean(dRati);
                %    beta=mean(Rati2);
                hlp=Z*Rati/N;
                %     mu=diag(W'*hlp)';

                %      J=ones(1,dim); tau=abs(nu-mu); gam=(beta-mu.^2);
                %      ISR=(gam'*J+J'*gam+J'*tau.^2)./(tau'*J+J'*tau).^2;
                %      ISR=ISR-diag(diag(ISR));

                W=hlp-ones(dim,1)*nu.*W;

            case 4%'gaus'
                U=Z'*W;
                Usquared=U.^2;
                ex=exp(-Usquared/2);
                gauss=U.*ex;
                dGauss=(1-Usquared).*ex;
                W=Z*gauss/N-ones(dim,1)*sum(dGauss).*W/N;
        end
        %   crit=abs(sum((W*diag(sqrt(1./sum(W.^2)))).*(W_before_decor*diag(sqrt(1./sum(W_before_decor.^2))))));
        W_before_decor=W;
        W=W*real(inv(W'*W)^(1/2));
        crit=abs(sum(W.*Wold));
        NumIt=NumIt+1;
    end %while iteration
    fprintf('EFICA v1.9 NumIt: %d\n',NumIt);
    repeat=0;
    %%%The SaddleTest of the separated components
    if SaddleTest
        SaddleTest=false; %%The SaddleTest may be done only one times
        u=Z'*W;
        switch test_of_saddle_points_nonln
            case 1 %'tanh'
                table1=(mean(log(cosh(u)))-0.37456).^2;
            case 2  %'gaus'
                table1=(mean(u.*exp(-u.^2))-1/sqrt(2)).^2;
            case 3 %'rati'
                table1=(mean(2*log(4+u.^2))-3.16).^2;
            case 4 %'rtsp'
                table1=(mean(u.^2./(1+abs(u)))-0.4125).^2;
            case 5 %'pow3'
                table1=(mean((pwr(u,4)))-3).^2;
        end
        rotated=zeros(1,dim);
        checked=1:dim;
        for i=checked
            for j=checked(checked>i)
                if (~rotated(i) && ~rotated(j))
                    h=[u(:,i) u(:,j)]*rot2d;   % rotate two columns i,j
                    switch test_of_saddle_points_nonln
                        case 1 %'tanh'
                            ctrl=(mean(log(cosh(h)))-0.37456).^2;
                        case 2 %'gaus'
                            ctrl=(mean(exp(-h.^2/2)-1/sqrt(2))).^2;
                        case 3 %'rati'
                            ctrl=(mean(2*log(4+h.^2))-3.16).^2;
                        case 4 %'rtsp'
                            ctrl=(mean(h.^2./(1+abs(h)))-0.4125).^2;
                        case 5 %'pow3'
                            ctrl=(mean((pwr(h,4)))-3).^2;
                    end
                    if sum(ctrl)>table1(i)+table1(j)
                        %bad extrem indicated
                        rotated(i)=1;rotated(j)=1; %do not test the rotated signals anymore
                        W(:,[i j])=W(:,[i j])*rot2d;
                        repeat=1; %continue in iterating - the test of saddle points is positive
                        NumIt=MaxIt-MaxItAfterSaddleTest;
                        fprintf('EFICA: rotating components: %d and %d\n',i,j);
                        status=1;
                    end
                end
            end
        end
    end %if SaddleTest
    crit=zeros(1,dim);
end %while repeat

Wsymm=W'*CC;

%estimated signals
s=W'*Z;

%estimate SIRs of the symmetric approach
switch g
    case  1 %'tanh'
        mu=mean(s.*tanh(s),2);
        nu=mean(1./cosh(s).^2,2);
        beta=mean(tanh(s).^2,2);
    case 2 %'rati'
        ssquare=s.^2;
        mu=mean(ssquare./(1+ssquare/4),2);
        nu=mean((1-ssquare/4)./(ssquare/4+1).^2,2);
        beta=mean(ssquare./(1+ssquare/4).^2,2);
    case 3 %'gaus'
        aexp=exp(-s.^2/2);
        mu=mean(s.^2.*aexp,2);
        nu=mean((1-s.^2).*aexp,2);
        beta=mean((s.*aexp).^2,2);