www.pudn.com > FastICA_2.1.zip > fpica.m
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
% 15.1.2001
% Hugo Gävert
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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(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,
if b_verbose
fprintf('No convergence after %d steps\n', maxNumIterations);
fprintf('Note that the plots are probably wrong.\n');
end
A=[];
W=[];
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: %.3f \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)');
end
case 2
if rem(round, displayInterval) == 0,
% ... and now there are :-)
% 1D basis
A = dewhiteningMatrix * B;
icaplot('dispsig',A');
end
case 3
if rem(round, displayInterval) == 0,
% ... and now there are :-)
% 1D filters
W = B' * whiteningMatrix;
icaplot('dispsig',W);
end
otherwise
end
drawnow;
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 = rand(vectorSize, 1) - .5;
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
drawnow;
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
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 = (Xub * ((Xub' * 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);