Source

ITE / code / shared / embedded / CFastICA / cfastica_public.m

Full commit
function [e_hat,W_hat] = cfastica_public(x,defl)
% Complex FastICA
% Ella Bingham 1999
% Neural Networks Research Centre, Helsinki University of Technology

% This is simple Matlab code for computing FastICA on complex valued signals.
% The algorithm is reported in:
% Ella Bingham and Aapo Hyvärinen, "A fast fixed-point algorithm for 
% independent component analysis of complex valued signals", International 
% Journal of Neural Systems, Vol. 10, No. 1 (February, 2000) 1-8.

% When using the code, please refer to the abovementioned publication.

% Nonlinearity G(y) = log(eps+y) is used in this code; this corresponds to 
% G_2 in the above paper with eps=a_2.

% Some bugs corrected on Oct 2003, thanks to Ioannis Andrianakis for pointing them out

% Compatible with GPLv2 or later. http://www.gnu.org/licenses/gpl-2.0.txt

% Copyright (C) 1999  Ella Bingham
%
% 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 in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA.
%
%Modification (Zoltan Szabo, Sep. 2012): whitening is carried out as a preprocessing step, separately; source generation, figures, plots: deleted.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eps = 0.1; % epsilon in G
%defl = 1; % components are estimated one by one in a deflationary manner; set this to 0 if you want them all estimated simultaneously
n = size(x,1);

% FIXED POINT ALGORITHM
if defl % Components estimated one by one
  W = zeros(n,n);
  maxcounter = 40;
  for k = 1:n
    w = rand(n,1) + i*rand(n,1);
    counter = 0;
    wold = zeros(n,1);
    while min(sum(abs(abs(wold) - abs(w))), maxcounter - counter) > 0.001;
      wold = w;
      g = 1./(eps + abs(w'*x).^2);
      dg = -1./(eps + abs(w'*x).^2).^2;
      w = mean(x .* (ones(n,1)*conj(w'*x)) .* (ones(n,1)*g), 2) - ...
	  mean(g + abs(w'*x).^2 .* dg) * w;
      w = w / norm(w);
      % Decorrelation:
      w = w - W*W'*w;
      w = w / norm(w);
      counter = counter + 1;
      G = log(eps + abs(w'*x).^2);
    end
    W(:,k) = w;
    counters(k) = counter;
  end; 

else %symmetric approach, all components estimated simultaneously
  maxcounter = 10;
  counter = 0;
  W = randn(n,n) + i*randn(n,n);
  while counter < maxcounter;
    for j = 1:n
      gWx(j,:) = 1./(eps + abs(W(:,j)'*x).^2);
      dgWx(j,:) = -1./(eps + abs(W(:,j)'*x).^2).^2;
      W(:,j) = mean(x .* (ones(n,1)*conj(W(:,j)'*x)) .* (ones(n,1)*gWx(j,:)),2) - mean(gWx(j,:) + abs(W(:,j)'*x).^2 .* dgWx(j,:)) * W(:,j);
    end;
    % Symmetric decorrelation:
    %W = W * sqrtm(inv(W'*W));
    [E,D] = eig(W'*W);
    W = W * E * inv(sqrt(D)) * E';
    counter = counter + 1;
    GWx = log(eps + abs(W'*x).^2);
  end
  
end

W_hat = W';
e_hat = W_hat * x;