% ==========================================================================
%
% This function performs ICA using fixed-point algorithm and the deflation
% scheme.
%
% Usage:  
%         >> [icasig,W] = Fasticadef(x);
% Inputs: 
%        x - the data  
%
% Outputs:
%   icasig - the independent components (estimated sources) 
%        W - the separating matrix
%
% For details, please see the following paper:
% A. Hyvarinen and E. Oja. Independent component
% analysis: Algorithms and applications. Neural
% Networks, 13(4-5):411–430, 2000.
%
%  The purpose of this software is the dissemination of
%  scientific work for scientific use. The commercial
%  distribution or use of this source code is prohibited. 
%
% by Zhenwei Shi
% shizhenwei@buaa.edu.cn
% 8-23-04
% ==========================================================================

function [icasig,W] = Fasticadef(x);

[N,P] = size(x);

% ------ sphering (whitening)
fprintf('\n sphering ... ');
x = x - mean(x')'*ones(1,P);               % subtract means from mixtures
mixedmean = mean(x')'; 
mixedsig = x;
whiteningMatrix = inv(sqrtm(cov(x')));              % get decorrelating matrix
dewhiteningMatrix = sqrtm(cov(x'));
whitesig = whiteningMatrix * mixedsig; 
X = whitesig;                                % decorrelate mixtures
fprintf(' done \n');

a1=1;
maxmaxNumIterations = 100;
[vectorSize,numSamples] = size(X);
B = zeros(vectorSize);
numOFIC = vectorSize;
epsilon = 0.00001;


for r = 1:numOFIC,
    i = 0;
    w = rand (vectorSize,1) - .5;
    w = w/norm(w); wOld = zeros(size(w));
    while i<=maxmaxNumIterations,
         
         if norm(w - wOld) < epsilon | norm(w + wOld) < epsilon,
            break;
         end   
            wOld = w;
            w = w - B * B' * w;
            w = w / norm(w);
            hypTan = tanh(a1 * X' * w);
            w = (X * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / numSamples;
           % Normalize the new w.
            w = w / norm(w);
            i = i + 1;
           
    end;
    i
    B(:,r) = w;
    W(r,:) = w' * whiteningMatrix;
	A(:,r) = dewhiteningMatrix * w;
end;
icasig = W * mixedsig + (W * mixedmean) * ones(1,numSamples);