--- a
+++ b/code/preprocessing/EEG/fastica/fpica.m
@@ -0,0 +1,910 @@
+function [A, W] = fpica(X, whiteningMatrix, dewhiteningMatrix, approach, ...
+			numOfIC, g, finetune, a1, a2, myy, stabilization, ...
+			epsilon, maxNumIterations, maxFinetune, initState, ...
+			guess, sampleSize, displayMode, displayInterval, ...
+			s_verbose);
+%FPICA - Fixed point ICA. Main algorithm of FASTICA.
+%
+% [A, W] = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach,
+%        numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, 
+%        maxNumIterations, maxFinetune, initState, guess, sampleSize,
+%        displayMode, displayInterval, verbose);
+% 
+% Perform independent component analysis using Hyvarinen's fixed point
+% algorithm. Outputs an estimate of the mixing matrix A and its inverse W.
+%
+% whitesig                              :the whitened data as row vectors
+% whiteningMatrix                       :can be obtained with function whitenv
+% dewhiteningMatrix                     :can be obtained with function whitenv
+% approach      [ 'symm' | 'defl' ]     :the approach used (deflation or symmetric)
+% numOfIC       [ 0 - Dim of whitesig ] :number of independent components estimated
+% g             [ 'pow3' | 'tanh' |     :the nonlinearity used
+%                 'gaus' | 'skew' ]     
+% finetune      [same as g + 'off']     :the nonlinearity used in finetuning.
+% a1                                    :parameter for tuning 'tanh'
+% a2                                    :parameter for tuning 'gaus'
+% mu                                    :step size in stabilized algorithm
+% stabilization [ 'on' | 'off' ]        :if mu < 1 then automatically on
+% epsilon                               :stopping criterion
+% maxNumIterations                      :maximum number of iterations 
+% maxFinetune                           :maximum number of iteretions for finetuning
+% initState     [ 'rand' | 'guess' ]    :initial guess or random initial state. See below
+% guess                                 :initial guess for A. Ignored if initState = 'rand'
+% sampleSize    [ 0 - 1 ]               :percentage of the samples used in one iteration
+% displayMode   [ 'signals' | 'basis' | :plot running estimate
+%                 'filters' | 'off' ]
+% displayInterval                       :number of iterations we take between plots
+% verbose       [ 'on' | 'off' ]        :report progress in text format
+%
+% EXAMPLE
+%       [E, D] = pcamat(vectors);
+%       [nv, wm, dwm] = whitenv(vectors, E, D);
+%       [A, W] = fpica(nv, wm, dwm);
+%
+%
+% This function is needed by FASTICA and FASTICAG
+%
+%   See also FASTICA, FASTICAG, WHITENV, PCAMAT
+
+% @(#)$Id$
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Global variable for stopping the ICA calculations from the GUI
+global g_FastICA_interrupt;
+if isempty(g_FastICA_interrupt)
+  clear global g_FastICA_interrupt;
+  interruptible = 0;
+else
+  interruptible = 1;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Default values
+
+if nargin < 3, error('Not enough arguments!'); end
+[vectorSize, numSamples] = size(X);
+if nargin < 20, s_verbose = 'on'; end
+if nargin < 19, displayInterval = 1; end
+if nargin < 18, displayMode = 'on'; end
+if nargin < 17, sampleSize = 1; end
+if nargin < 16, guess = 1; end
+if nargin < 15, initState = 'rand'; end
+if nargin < 14, maxFinetune = 100; end
+if nargin < 13, maxNumIterations = 1000; end
+if nargin < 12, epsilon = 0.0001; end
+if nargin < 11, stabilization = 'on'; end
+if nargin < 10, myy = 1; end
+if nargin < 9, a2 = 1; end
+if nargin < 8, a1 = 1; end
+if nargin < 7, finetune = 'off'; end
+if nargin < 6, g = 'pow3'; end
+if nargin < 5, numOfIC = vectorSize; end     % vectorSize = Dim
+if nargin < 4, approach = 'defl'; end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the data
+
+if ~isreal(X)
+  error('Input has an imaginary part.');
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the value for verbose
+
+switch lower(s_verbose)
+ case 'on'
+  b_verbose = 1;
+ case 'off'
+  b_verbose = 0;
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the value for approach
+
+switch lower(approach)
+ case 'symm'
+  approachMode = 1;
+ case 'defl'
+  approachMode = 2;
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''approach''\n', approach));
+end
+if b_verbose, fprintf('Used approach [ %s ].\n', approach); end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the value for numOfIC
+
+if vectorSize < numOfIC
+  error('Must have numOfIC <= Dimension!');
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the sampleSize
+if sampleSize > 1
+  sampleSize = 1;
+  if b_verbose
+    fprintf('Warning: Setting ''sampleSize'' to 1.\n');
+  end  
+elseif sampleSize < 1
+  if (sampleSize * numSamples) < 1000
+    sampleSize = min(1000/numSamples, 1);
+    if b_verbose
+      fprintf('Warning: Setting ''sampleSize'' to %0.3f (%d samples).\n', ...
+	      sampleSize, floor(sampleSize * numSamples));
+    end  
+  end
+end
+if b_verbose
+  if  b_verbose & (sampleSize < 1)
+    fprintf('Using about %0.0f%% of the samples in random order in every step.\n',sampleSize*100);
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the value for nonlinearity.
+
+switch lower(g)
+ case 'pow3'
+  gOrig = 10;
+ case 'tanh'
+  gOrig = 20;
+ case {'gaus', 'gauss'}
+  gOrig = 30;
+ case 'skew'
+  gOrig = 40;
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''g''\n', g));
+end
+if sampleSize ~= 1
+  gOrig = gOrig + 2;
+end
+if myy ~= 1
+  gOrig = gOrig + 1;
+end
+
+if b_verbose,
+  fprintf('Used nonlinearity [ %s ].\n', g);
+end
+
+finetuningEnabled = 1;
+switch lower(finetune)
+ case 'pow3'
+  gFine = 10 + 1;
+ case 'tanh'
+  gFine = 20 + 1;
+ case {'gaus', 'gauss'}
+  gFine = 30 + 1;
+ case 'skew'
+  gFine = 40 + 1;
+ case 'off'
+  if myy ~= 1
+    gFine = gOrig;
+  else 
+    gFine = gOrig + 1;
+  end
+  finetuningEnabled = 0;
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''finetune''\n', ...
+		finetune));
+end
+
+if b_verbose & finetuningEnabled
+  fprintf('Finetuning enabled (nonlinearity: [ %s ]).\n', finetune);
+end
+
+switch lower(stabilization)
+ case 'on'
+  stabilizationEnabled = 1;
+ case 'off'
+  if myy ~= 1
+    stabilizationEnabled = 1;
+  else
+    stabilizationEnabled = 0;
+  end
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''stabilization''\n', ...
+		stabilization)); 
+end
+
+if b_verbose & stabilizationEnabled
+  fprintf('Using stabilized algorithm.\n');
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Some other parameters
+myyOrig = myy;
+% When we start fine-tuning we'll set myy = myyK * myy
+myyK = 0.01;
+% How many times do we try for convergence until we give up.
+failureLimit = 5;
+
+
+usedNlinearity = gOrig;
+stroke = 0;
+notFine = 1;
+long = 0;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the value for initial state.
+
+switch lower(initState)
+ case 'rand'
+  initialStateMode = 0;
+ case 'guess'
+  if size(guess,1) ~= size(whiteningMatrix,2)
+    initialStateMode = 0;
+    if b_verbose
+      fprintf('Warning: size of initial guess is incorrect. Using random initial guess.\n');
+    end
+  else
+    initialStateMode = 1;
+    if size(guess,2) < numOfIC
+      if b_verbose
+	fprintf('Warning: initial guess only for first %d components. Using random initial guess for others.\n', size(guess,2)); 
+      end
+      guess(:, size(guess, 2) + 1:numOfIC) = ...
+					     rand(vectorSize,numOfIC-size(guess,2))-.5;
+    elseif size(guess,2)>numOfIC
+      guess=guess(:,1:numOfIC);
+      fprintf('Warning: Initial guess too large. The excess column are dropped.\n');
+    end
+    if b_verbose, fprintf('Using initial guess.\n'); end
+  end
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''initState''\n', initState));
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Checking the value for display mode.
+
+switch lower(displayMode)
+ case {'off', 'none'}
+  usedDisplay = 0;
+ case {'on', 'signals'}
+  usedDisplay = 1;
+  if (b_verbose & (numSamples > 10000))
+    fprintf('Warning: Data vectors are very long. Plotting may take long time.\n');
+  end
+  if (b_verbose & (numOfIC > 25))
+    fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
+  end
+ case 'basis'
+  usedDisplay = 2;
+  if (b_verbose & (numOfIC > 25))
+    fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
+  end
+ case 'filters'
+  usedDisplay = 3;
+  if (b_verbose & (vectorSize > 25))
+    fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
+  end
+ otherwise
+  error(sprintf('Illegal value [ %s ] for parameter: ''displayMode''\n', displayMode));
+end
+
+% The displayInterval can't be less than 1...
+if displayInterval < 1
+  displayInterval = 1;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if b_verbose, fprintf('Starting ICA calculation...\n'); end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SYMMETRIC APPROACH
+if approachMode == 1,
+
+  % set some parameters more...
+  usedNlinearity = gOrig;
+  stroke = 0;
+  notFine = 1;
+  long = 0;
+  
+  A = zeros(vectorSize, numOfIC);  % Dewhitened basis vectors.
+  if initialStateMode == 0
+    % Take random orthonormal initial vectors.
+    B = orth (randn (vectorSize, numOfIC));
+  elseif initialStateMode == 1
+    % Use the given initial vector as the initial state
+    B = whiteningMatrix * guess;
+  end
+  
+  BOld = zeros(size(B));
+  BOld2 = zeros(size(B));
+  
+  % This is the actual fixed-point iteration loop.
+  for round = 1:maxNumIterations + 1,
+    if round == maxNumIterations + 1,
+      fprintf('No convergence after %d steps\n', maxNumIterations);
+      fprintf('Note that the plots are probably wrong.\n');
+      if ~isempty(B)
+	% Symmetric orthogonalization.
+	B = B * real(inv(B' * B)^(1/2));
+
+	W = B' * whiteningMatrix;
+	A = dewhiteningMatrix * B;
+      else
+	W = [];
+	A = [];
+      end
+      return;
+    end
+    
+    if (interruptible & g_FastICA_interrupt)
+      if b_verbose 
+        fprintf('\n\nCalculation interrupted by the user\n');
+      end
+      if ~isempty(B)
+	W = B' * whiteningMatrix;
+	A = dewhiteningMatrix * B;
+      else
+	W = [];
+	A = [];
+      end
+      return;
+    end
+    
+    
+    % Symmetric orthogonalization.
+    B = B * real(inv(B' * B)^(1/2));
+    
+    % Test for termination condition. Note that we consider opposite
+    % directions here as well.
+    minAbsCos = min(abs(diag(B' * BOld)));
+    minAbsCos2 = min(abs(diag(B' * BOld2)));
+    
+    if (1 - minAbsCos < epsilon)
+      if finetuningEnabled & notFine
+        if b_verbose, fprintf('Initial convergence, fine-tuning: \n'); end;
+        notFine = 0;
+        usedNlinearity = gFine;
+        myy = myyK * myyOrig;
+        BOld = zeros(size(B));
+        BOld2 = zeros(size(B));
+	
+      else
+        if b_verbose, fprintf('Convergence after %d steps\n', round); end
+	
+        % Calculate the de-whitened vectors.
+        A = dewhiteningMatrix * B;
+        break;
+      end
+    elseif stabilizationEnabled
+      if (~stroke) & (1 - minAbsCos2 < epsilon)
+	if b_verbose, fprintf('Stroke!\n'); end;
+	stroke = myy;
+	myy = .5*myy;
+	if mod(usedNlinearity,2) == 0
+	  usedNlinearity = usedNlinearity + 1;
+	end
+      elseif stroke
+	myy = stroke;
+	stroke = 0;
+	if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
+	  usedNlinearity = usedNlinearity - 1;
+	end
+      elseif (~long) & (round>maxNumIterations/2)
+	if b_verbose, fprintf('Taking long (reducing step size)\n'); end;
+	long = 1;
+	myy = .5*myy;
+	if mod(usedNlinearity,2) == 0
+	  usedNlinearity = usedNlinearity + 1;
+	end
+      end
+    end
+    
+    BOld2 = BOld;
+    BOld = B;
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Show the progress...
+    if b_verbose
+      if round == 1
+        fprintf('Step no. %d\n', round);
+      else
+        fprintf('Step no. %d, change in value of estimate: %.3g \n', round, 1 - minAbsCos);
+      end
+    end
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Also plot the current state...
+    switch usedDisplay
+     case 1
+      if rem(round, displayInterval) == 0,
+	% There was and may still be other displaymodes...
+	% 1D signals
+	icaplot('dispsig',(X'*B)');
+	drawnow;
+      end
+     case 2
+      if rem(round, displayInterval) == 0,
+	% ... and now there are :-)
+	% 1D basis
+	A = dewhiteningMatrix * B;
+	icaplot('dispsig',A');
+	drawnow;
+      end
+     case 3
+      if rem(round, displayInterval) == 0,
+	% ... and now there are :-)
+	% 1D filters
+	W = B' * whiteningMatrix;
+	icaplot('dispsig',W);
+	drawnow;
+      end
+     otherwise
+    end
+    
+    switch usedNlinearity
+      % pow3
+     case 10
+      B = (X * (( X' * B) .^ 3)) / numSamples - 3 * B;
+     case 11
+      % optimoitu - epsilonin kokoisia eroja
+      % tämä on optimoitu koodi, katso vanha koodi esim.
+      % aikaisemmista versioista kuten 2.0 beta3
+      Y = X' * B;
+      Gpow3 = Y .^ 3;
+      Beta = sum(Y .* Gpow3);
+      D = diag(1 ./ (Beta - 3 * numSamples));
+      B = B + myy * B * (Y' * Gpow3 - diag(Beta)) * D;
+     case 12
+      Xsub=X(:, getSamples(numSamples, sampleSize));
+      B = (Xsub * (( Xsub' * B) .^ 3)) / size(Xsub,2) - 3 * B;
+     case 13
+      % Optimoitu
+      Ysub=X(:, getSamples(numSamples, sampleSize))' * B;
+      Gpow3 = Ysub .^ 3;
+      Beta = sum(Ysub .* Gpow3);
+      D = diag(1 ./ (Beta - 3 * size(Ysub', 2)));
+      B = B + myy * B * (Ysub' * Gpow3 - diag(Beta)) * D;
+      
+      % tanh
+     case 20
+      hypTan = tanh(a1 * X' * B);
+      B = X * hypTan / numSamples - ...
+	  ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / numSamples * ...
+	  a1;
+     case 21
+      % optimoitu - epsilonin kokoisia 
+      Y = X' * B;
+      hypTan = tanh(a1 * Y);
+      Beta = sum(Y .* hypTan);
+      D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
+      B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
+     case 22
+      Xsub=X(:, getSamples(numSamples, sampleSize));
+      hypTan = tanh(a1 * Xsub' * B);
+      B = Xsub * hypTan / size(Xsub, 2) - ...
+	  ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / size(Xsub, 2) * a1;
+     case 23
+      % Optimoitu
+      Y = X(:, getSamples(numSamples, sampleSize))' * B;
+      hypTan = tanh(a1 * Y);
+      Beta = sum(Y .* hypTan);
+      D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
+      B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
+      
+      % gauss
+     case 30
+      U = X' * B;
+      Usquared=U .^ 2;
+      ex = exp(-a2 * Usquared / 2);
+      gauss =  U .* ex;
+      dGauss = (1 - a2 * Usquared) .*ex;
+      B = X * gauss / numSamples - ...
+	  ones(size(B,1),1) * sum(dGauss)...
+	  .* B / numSamples ;
+     case 31
+      % optimoitu
+      Y = X' * B;
+      ex = exp(-a2 * (Y .^ 2) / 2);
+      gauss = Y .* ex;
+      Beta = sum(Y .* gauss);
+      D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
+      B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
+     case 32
+      Xsub=X(:, getSamples(numSamples, sampleSize));
+      U = Xsub' * B;
+      Usquared=U .^ 2;
+      ex = exp(-a2 * Usquared / 2);
+      gauss =  U .* ex;
+      dGauss = (1 - a2 * Usquared) .*ex;
+      B = Xsub * gauss / size(Xsub,2) - ...
+	  ones(size(B,1),1) * sum(dGauss)...
+	  .* B / size(Xsub,2) ;
+     case 33
+      % Optimoitu
+      Y = X(:, getSamples(numSamples, sampleSize))' * B;
+      ex = exp(-a2 * (Y .^ 2) / 2);
+      gauss = Y .* ex;
+      Beta = sum(Y .* gauss);
+      D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
+      B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
+      
+      % skew
+     case 40
+      B = (X * ((X' * B) .^ 2)) / numSamples;
+     case 41
+      % Optimoitu
+      Y = X' * B;
+      Gskew = Y .^ 2;
+      Beta = sum(Y .* Gskew);
+      D = diag(1 ./ (Beta));
+      B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
+     case 42
+      Xsub=X(:, getSamples(numSamples, sampleSize));
+      B = (Xsub * ((Xsub' * B) .^ 2)) / size(Xsub,2);
+     case 43
+      % Uusi optimoitu
+      Y = X(:, getSamples(numSamples, sampleSize))' * B;
+      Gskew = Y .^ 2;
+      Beta = sum(Y .* Gskew);
+      D = diag(1 ./ (Beta));
+      B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
+
+     otherwise
+      error('Code for desired nonlinearity not found!');
+    end
+  end
+
+  
+  % Calculate ICA filters.
+  W = B' * whiteningMatrix;
+  
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Also plot the last one...
+  switch usedDisplay
+   case 1 
+    % There was and may still be other displaymodes...
+    % 1D signals
+    icaplot('dispsig',(X'*B)');
+    drawnow;
+   case 2
+    % ... and now there are :-)
+    % 1D basis
+    icaplot('dispsig',A');
+    drawnow;
+   case 3
+    % ... and now there are :-)
+    % 1D filters
+    icaplot('dispsig',W);
+    drawnow;
+   otherwise
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% DEFLATION APPROACH
+if approachMode == 2
+  
+  B = zeros(vectorSize);
+  
+  % The search for a basis vector is repeated numOfIC times.
+  round = 1;
+  
+  numFailures = 0;
+  
+  while round <= numOfIC,
+    myy = myyOrig;
+    usedNlinearity = gOrig;
+    stroke = 0;
+    notFine = 1;
+    long = 0;
+    endFinetuning = 0;
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Show the progress...
+    if b_verbose, fprintf('IC %d ', round); end
+    
+    % Take a random initial vector of lenght 1 and orthogonalize it
+    % with respect to the other vectors.
+    if initialStateMode == 0
+      w = randn (vectorSize, 1);
+    elseif initialStateMode == 1
+      w=whiteningMatrix*guess(:,round);
+    end
+    w = w - B * B' * w;
+    w = w / norm(w);
+    
+    wOld = zeros(size(w));
+    wOld2 = zeros(size(w));
+    
+    % This is the actual fixed-point iteration loop.
+    %    for i = 1 : maxNumIterations + 1
+    i = 1;
+    gabba = 1;
+    while i <= maxNumIterations + gabba
+      if (usedDisplay > 0)
+	drawnow;
+      end
+      if (interruptible & g_FastICA_interrupt)
+        if b_verbose 
+          fprintf('\n\nCalculation interrupted by the user\n');
+        end
+        return;
+      end
+      
+      % Project the vector into the space orthogonal to the space
+      % spanned by the earlier found basis vectors. Note that we can do
+      % the projection with matrix B, since the zero entries do not
+      % contribute to the projection.
+      w = w - B * B' * w;
+      w = w / norm(w);
+      
+      if notFine
+	if i == maxNumIterations + 1
+	  if b_verbose
+	    fprintf('\nComponent number %d did not converge in %d iterations.\n', round, maxNumIterations);
+	  end
+	  round = round - 1;
+	  numFailures = numFailures + 1;
+	  if numFailures > failureLimit
+	    if b_verbose
+	      fprintf('Too many failures to converge (%d). Giving up.\n', numFailures);
+	    end
+	    if round == 0
+	      A=[];
+	      W=[];
+	    end
+	    return;
+	  end
+	  % numFailures > failurelimit
+	  break;
+	end
+	% i == maxNumIterations + 1
+      else
+	% if notFine
+	if i >= endFinetuning
+	  wOld = w; % So the algorithm will stop on the next test...
+	end
+      end
+      % if notFine
+      
+      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      % Show the progress...
+      if b_verbose, fprintf('.'); end;
+      
+      
+      % Test for termination condition. Note that the algorithm has
+      % converged if the direction of w and wOld is the same, this
+      % is why we test the two cases.
+      if norm(w - wOld) < epsilon | norm(w + wOld) < epsilon
+        if finetuningEnabled & notFine
+          if b_verbose, fprintf('Initial convergence, fine-tuning: '); end;
+          notFine = 0;
+	  gabba = maxFinetune;
+          wOld = zeros(size(w));
+          wOld2 = zeros(size(w));
+          usedNlinearity = gFine;
+          myy = myyK * myyOrig;
+	  
+	  endFinetuning = maxFinetune + i;
+	  
+        else
+          numFailures = 0;
+          % Save the vector
+          B(:, round) = w;
+	  
+          % Calculate the de-whitened vector.
+          A(:,round) = dewhiteningMatrix * w;
+          % Calculate ICA filter.
+          W(round,:) = w' * whiteningMatrix;
+	  
+          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+          % Show the progress...
+          if b_verbose, fprintf('computed ( %d steps ) \n', i); end
+	  
+          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+          % Also plot the current state...
+          switch usedDisplay
+	   case 1
+	    if rem(round, displayInterval) == 0,
+	      % There was and may still be other displaymodes...   
+	      % 1D signals
+	      temp = X'*B;
+	      icaplot('dispsig',temp(:,1:numOfIC)');
+	      drawnow;
+	    end
+	   case 2
+	    if rem(round, displayInterval) == 0,
+	      % ... and now there are :-) 
+	      % 1D basis
+	      icaplot('dispsig',A');
+	      drawnow;
+	    end
+	   case 3
+	    if rem(round, displayInterval) == 0,
+	      % ... and now there are :-) 
+	      % 1D filters
+	      icaplot('dispsig',W);
+	      drawnow;
+	    end
+          end
+	  % switch usedDisplay
+	  break; % IC ready - next...
+        end
+	%if finetuningEnabled & notFine
+      elseif stabilizationEnabled
+	if (~stroke) & (norm(w - wOld2) < epsilon | norm(w + wOld2) < ...
+			epsilon)
+	  stroke = myy;
+	  if b_verbose, fprintf('Stroke!'); end;
+	  myy = .5*myy;
+	  if mod(usedNlinearity,2) == 0
+	    usedNlinearity = usedNlinearity + 1;
+	  end
+	elseif stroke
+	  myy = stroke;
+	  stroke = 0;
+	  if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
+	    usedNlinearity = usedNlinearity - 1;
+	  end
+	elseif (notFine) & (~long) & (i > maxNumIterations / 2)
+	  if b_verbose, fprintf('Taking long (reducing step size) '); end;
+	  long = 1;
+	  myy = .5*myy;
+	  if mod(usedNlinearity,2) == 0
+	    usedNlinearity = usedNlinearity + 1;
+	  end
+	end
+      end
+      
+      wOld2 = wOld;
+      wOld = w;
+      
+      switch usedNlinearity
+	% pow3
+       case 10
+	w = (X * ((X' * w) .^ 3)) / numSamples - 3 * w;
+       case 11
+	EXGpow3 = (X * ((X' * w) .^ 3)) / numSamples;
+	Beta = w' * EXGpow3;
+	w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
+       case 12
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	w = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2) - 3 * w;
+       case 13
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	EXGpow3 = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2);
+	Beta = w' * EXGpow3;
+	w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
+	% tanh
+       case 20
+    % #########################################################################
+    % Riccardo Navarra "Mod" 10 Feb 2010 (*)          (r.navarra@itab.unich.it)
+    % Rewritten matrix product to speed up on multicore CPU (fpica.m)
+    % #########################################################################
+	hypTan = tanh(a1 * (X' * w));
+    % #########################################################################
+	w = (X * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / numSamples;
+       case 21
+	hypTan = tanh(a1 * X' * w);
+	Beta = w' * X * hypTan;
+	w = w - myy * ((X * hypTan - Beta * w) / ...
+		       (a1 * sum((1-hypTan .^2)') - Beta));
+       case 22
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	hypTan = tanh(a1 * Xsub' * w);
+	w = (Xsub * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / size(Xsub, 2);
+       case 23
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	hypTan = tanh(a1 * Xsub' * w);
+	Beta = w' * Xsub * hypTan;
+	w = w - myy * ((Xsub * hypTan - Beta * w) / ...
+		       (a1 * sum((1-hypTan .^2)') - Beta));
+	% gauss
+       case 30
+	% This has been split for performance reasons.
+	u = X' * w;
+	u2=u.^2;
+	ex=exp(-a2 * u2/2);
+	gauss =  u.*ex;
+	dGauss = (1 - a2 * u2) .*ex;
+	w = (X * gauss - sum(dGauss)' * w) / numSamples;
+       case 31
+	u = X' * w;
+	u2=u.^2;
+	ex=exp(-a2 * u2/2);
+	gauss =  u.*ex;
+	dGauss = (1 - a2 * u2) .*ex;
+	Beta = w' * X * gauss;
+	w = w - myy * ((X * gauss - Beta * w) / ...
+		       (sum(dGauss)' - Beta));
+       case 32
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	u = Xsub' * w;
+	u2=u.^2;
+	ex=exp(-a2 * u2/2);
+	gauss =  u.*ex;
+	dGauss = (1 - a2 * u2) .*ex;
+	w = (Xsub * gauss - sum(dGauss)' * w) / size(Xsub, 2);
+       case 33
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	u = Xsub' * w;
+	u2=u.^2;
+	ex=exp(-a2 * u2/2);
+	gauss =  u.*ex;
+	dGauss = (1 - a2 * u2) .*ex;
+	Beta = w' * Xsub * gauss;
+	w = w - myy * ((Xsub * gauss - Beta * w) / ...
+		       (sum(dGauss)' - Beta));
+	% skew
+       case 40
+	w = (X * ((X' * w) .^ 2)) / numSamples;
+       case 41
+	EXGskew = (X * ((X' * w) .^ 2)) / numSamples;
+	Beta = w' * EXGskew;
+	w = w - myy * (EXGskew - Beta*w)/(-Beta);
+       case 42
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	w = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
+       case 43
+	Xsub=X(:,getSamples(numSamples, sampleSize));
+	EXGskew = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
+	Beta = w' * EXGskew;
+	w = w - myy * (EXGskew - Beta*w)/(-Beta);
+	
+       otherwise
+	error('Code for desired nonlinearity not found!');
+      end
+      
+      % Normalize the new w.
+      w = w / norm(w);
+      i = i + 1;
+    end
+    round = round + 1;
+  end
+  if b_verbose, fprintf('Done.\n'); end
+  
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Also plot the ones that may not have been plotted.
+  if (usedDisplay > 0) & (rem(round-1, displayInterval) ~= 0)
+    switch usedDisplay
+     case 1
+      % There was and may still be other displaymodes...
+      % 1D signals
+      temp = X'*B;
+      icaplot('dispsig',temp(:,1:numOfIC)');
+      drawnow;
+     case 2
+      % ... and now there are :-)
+      % 1D basis
+      icaplot('dispsig',A');
+      drawnow;
+     case 3
+      % ... and now there are :-)
+      % 1D filters
+      icaplot('dispsig',W);
+      drawnow;
+     otherwise
+    end
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% In the end let's check the data for some security
+if ~isreal(A)
+  if b_verbose, fprintf('Warning: removing the imaginary part from the result.\n'); end
+  A = real(A);
+  W = real(W);
+end
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Subfunction
+% Calculates tanh simplier and faster than Matlab tanh.
+function y=tanh(x)
+y = 1 - 2 ./ (exp(2 * x) + 1);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function Samples = getSamples(max, percentage)
+Samples = find(rand(1, max) < percentage);