function [ica,Ae,cerr]=kurtFpDeIca(X,epsilon,maxNumIter,debug)
% --------------------------------------------------------------------------------------------
% Fixed point algorithm (FastICA, Hyvarinen) using kurtosis
%
% function [ica,Ae,cerr]=kurtFpDeIca(X,epsilon,maxNumIter,debug)
% X : mixed signals (each row is a signal, must be centered)
% epsilon : convergence stopping criterion (default=0.0001)
% maxNumIter : maximum number of iterations (default=200)
% debug : display each iteration of algorithm
%
% ica : separated signals
% Ae : estimated mixing matrix (each row is a component)
% cerr : convergence error (1-> convergence failure)
%
% Simple version using as many variables as sources
% and deflationay orthogonalization
%
% --------------------------------------------------------------------------------------------
% Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% --------------------------------------------------------------------------------------------
% This program is free software; you can redistribute it and/or modify it under the terms
% of the GNU General Public License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.
%
% This program is distributed "as is" and "as available" in the hope that it will be useful,
% but WITHOUT ANY WARRANTY of any kind; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
% --------------------------------------------------------------------------------------------
if(nargin<2), epsilon= 0.0001; end
if(nargin<3), maxNumIter=200; end
if(nargin<4), debug=0; end
[numVars, numSamples] = size(X);
cerr = 0;
% -------------------------- Whitening matrix estimation -----------------------
Rx=(X*X')/numSamples; % covariance matrix
% Calculate the eigenvalues and eigenvectors of covariance matrix.
[E, D] = eig (Rx);
W = inv(sqrt (D)) * E'; % whitening matrix
% W = E*inv(sqrt (D))*E'; % whitening matrix (W=sqrtm(Rx))
Z=W*X; % Whitened data
numOfIC=numVars;
fprintf('FastIca, kurt , deflationary \n');
fprintf('epsilon=%12.6f \n', epsilon);
W = zeros(numVars);
% The search for an independent component is repeated numOfIC times.
for ic = 1:numOfIC,
fprintf('Component=%3d\n',ic);
% Take a random initial vector
w = rand(numVars, 1) - .5;
w = w / norm(w);
% wOld = w + epsilon;
wOld = zeros(size(w));
% Fixed-point iteration loop for a single IC.
for iter = 1 : maxNumIter + 1
% Test for termination condition.
% The algorithm converged if the direction of w and wOld is the same.
absCos = abs(w' * wOld);
if(1-absCos < epsilon), break, end
if(debug), fprintf('it=%3d,%9.6f; ',iter, 1-absCos);end
wOld=w;
% Fixed point equation:
% the Kurtosis gradient, on the right-hand side gives the new value for w
w = (Z * ((Z' * w) .^ 3)) / numSamples - 3 * w;
% Deflationary orthogonalization.
% The current w vector is projected into the space orthogonal to the space
% spanned by the previuosly found W vectors.
w = w - W * W' * w;
w = w / norm(w);
end
if(debug), fprintf('\'); end
fprintf('numIt=%3d,%9.6f\n',iter, 1-absCos);
% Save the w vector
W(:, ic) = w;
cerr = cerr | (1-absCos >= epsilon);
if(cerr), fprintf(' ==> Convergence failure!\n');end
end
% if(cerr), Ae=[]; ica=[]; return; end
Ae = E*sqrt(D)*W; % estimated mixing matrix
ica = W'*Z; % separated components
end %== function ================================================================
%