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: fpica.m,v 1.6 2004/07/27 11:37:17 jarmo Exp $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global variable for stopping the ICA calculations from the GUIglobal g_FastICA_interrupt;if isempty(g_FastICA_interrupt)  clear global g_FastICA_interrupt;  interruptible = 0;else  interruptible = 1;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Default valuesif nargin < 3, error('Not enough arguments!'); end[vectorSize, numSamples] = size(X);if nargin < 20, s_verbose = 'on'; endif nargin < 19, displayInterval = 1; endif nargin < 18, displayMode = 'on'; endif nargin < 17, sampleSize = 1; endif nargin < 16, guess = 1; endif nargin < 15, initState = 'rand'; endif nargin < 14, maxFinetune = 100; endif nargin < 13, maxNumIterations = 1000; endif nargin < 12, epsilon = 0.0001; endif nargin < 11, stabilization = 'on'; endif nargin < 10, myy = 1; endif nargin < 9, a2 = 1; endif nargin < 8, a1 = 1; endif nargin < 7, finetune = 'off'; endif nargin < 6, g = 'pow3'; endif nargin < 5, numOfIC = vectorSize; end     % vectorSize = Dimif nargin < 4, approach = 'defl'; end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Checking the dataif ~isreal(X)  error('Input has an imaginary part.');end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Checking the value for verboseswitch 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 approachswitch lower(approach) case 'symm'  approachMode = 1; case 'defl'  approachMode = 2; otherwise  error(sprintf('Illegal value [ %s ] for parameter: ''approach''\n', approach));endif b_verbose, fprintf('Used approach [ %s ].\n', approach); end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Checking the value for numOfICif vectorSize < numOfIC  error('Must have numOfIC <= Dimension!');end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Checking the sampleSizeif 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    endendif b_verbose  if  b_verbose & (sampleSize < 1)    fprintf('Using about %0.0f%% of the samples in random order in every step.\n',sampleSize*100);  endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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));endif sampleSize ~= 1  gOrig = gOrig + 2;endif myy ~= 1  gOrig = gOrig + 1;endif b_verbose,  fprintf('Used nonlinearity [ %s ].\n', g);endfinetuningEnabled = 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));endif b_verbose & finetuningEnabled  fprintf('Finetuning enabled (nonlinearity: [ %s ]).\n', finetune);endswitch 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)); endif b_verbose & stabilizationEnabled  fprintf('Using stabilized algorithm.\n');end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Some other parametersmyyOrig = myy;% When we start fine-tuning we'll set myy = myyK * myymyyK = 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 APPROACHif 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(rand(vectorSize, numOfIC) - .5);  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)	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鋗