function [ica,Ae,cerr]=coshFpDeIca(X,epsilon,maxNumIter,debug)
% -------------------------------------------------------------------------------------------
% Fixed point algorithm (FastICA, Hyvarinen) using hyperbolic cosine (log(cosh())
%
% function [ica,Ae,cerr]=coshFpDeIca(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;
a1=1;
fprintf('FastIca, cosh , 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
% minNormDiff=min(norm(w - wOld), norm(w + wOld));
% if minNormDiff < epsilon, break, end
% fprintf('Iteration= %3d, %12.6f, eps=%12.6f \n',iter, minNormDiff, epsilon);
wOld=w;
% Fixed point equation:
% w <-- E[x g(w'x)] - E[g'(w'x)]w, where g' is the derivative of g
% chosing G(y)= 1/a log(cosh(ay)) we have: g=tanh(ay)) and g'=a(1 -tanh^2(ay))
% w <-- E[x tanh(a w'x)] - E[a (1-tanh^2(a w'x))] w
Tanh = tanh(a1 * Z' * w);
w = (Z * Tanh) / numSamples - a1 * sum(1 - Tanh .^ 2)'/ numSamples * 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('\n'); 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 ================================================================
%