515 lines
18 KiB
Matlab
515 lines
18 KiB
Matlab
function [Out1, Out2, Out3] = fastica(mixedsig, varargin)
|
||
%FASTICA - Fast Independent Component Analysis
|
||
%
|
||
% FastICA for Matlab 5.x
|
||
% Version 2.1, January 15 2001
|
||
% Copyright (c) Hugo G<>vert, Jarmo Hurri, Jaakko S<>rel<65>, and Aapo Hyv<79>rinen.
|
||
%
|
||
% FASTICA(mixedsig) estimates the independent components from given
|
||
% multidimensional signals. Each row of matrix mixedsig is one
|
||
% observed signal. FASTICA uses Hyvarinen's fixed-point algorithm,
|
||
% see http://www.cis.hut.fi/projects/ica/fastica/. Output from the
|
||
% function depends on the number output arguments:
|
||
%
|
||
% [icasig] = FASTICA (mixedsig); the rows of icasig contain the
|
||
% estimated independent components.
|
||
%
|
||
% [icasig, A, W] = FASTICA (mixedsig); outputs the estimated separating
|
||
% matrix W and the corresponding mixing matrix A.
|
||
%
|
||
% [A, W] = FASTICA (mixedsig); gives only the estimated mixing matrix
|
||
% A and the separating matrix W.
|
||
%
|
||
% Some optional arguments induce other output formats, see below.
|
||
%
|
||
% A graphical user interface for FASTICA can be launched by the
|
||
% command FASTICAG
|
||
%
|
||
% FASTICA can be called with numerous optional arguments. Optional
|
||
% arguments are given in parameter pairs, so that first argument is
|
||
% the name of the parameter and the next argument is the value for
|
||
% that parameter. Optional parameter pairs can be given in any order.
|
||
%
|
||
% OPTIONAL PARAMETERS:
|
||
%
|
||
% Parameter name Values and description
|
||
%
|
||
%======================================================================
|
||
% --Basic parameters in fixed-point algorithm:
|
||
%
|
||
% 'approach' (string) The decorrelation approach used. Can be
|
||
% symmetric ('symm'), i.e. estimate all the
|
||
% independent component in parallel, or
|
||
% deflation ('defl'), i.e. estimate independent
|
||
% component one-by-one like in projection pursuit.
|
||
% Default is 'defl'.
|
||
%
|
||
% 'numOfIC' (integer) Number of independent components to
|
||
% be estimated. Default equals the dimension of data.
|
||
%
|
||
%======================================================================
|
||
% --Choosing the nonlinearity:
|
||
%
|
||
% 'g' (string) Chooses the nonlinearity g used in
|
||
% the fixed-point algorithm. Possible values:
|
||
%
|
||
% Value of 'g': Nonlinearity used:
|
||
% 'pow3' (default) g(u)=u^3
|
||
% 'tanh' g(u)=tanh(a1*u)
|
||
% 'gauss g(u)=u*exp(-a2*u^2/2)
|
||
% 'skew' g(u)=u^2
|
||
%
|
||
% 'finetune' (string) Chooses the nonlinearity g used when
|
||
% fine-tuning. In addition to same values
|
||
% as for 'g', the possible value 'finetune' is:
|
||
% 'off' fine-tuning is disabled.
|
||
%
|
||
% 'a1' (number) Parameter a1 used when g='tanh'.
|
||
% Default is 1.
|
||
% 'a2' (number) Parameter a2 used when g='gaus'.
|
||
% Default is 1.
|
||
%
|
||
% 'mu' (number) Step size. Default is 1.
|
||
% If the value of mu is other than 1, then the
|
||
% program will use the stabilized version of the
|
||
% algorithm (see also parameter 'stabilization').
|
||
%
|
||
%
|
||
% 'stabilization' (string) Values 'on' or 'off'. Default 'off'.
|
||
% This parameter controls wether the program uses
|
||
% the stabilized version of the algorithm or
|
||
% not. If the stabilization is on, then the value
|
||
% of mu can momentarily be halved if the program
|
||
% senses that the algorithm is stuck between two
|
||
% points (this is called a stroke). Also if there
|
||
% is no convergence before half of the maximum
|
||
% number of iterations has been reached then mu
|
||
% will be halved for the rest of the rounds.
|
||
%
|
||
%======================================================================
|
||
% --Controlling convergence:
|
||
%
|
||
% 'epsilon' (number) Stopping criterion. Default is 0.0001.
|
||
%
|
||
% 'maxNumIterations' (integer) Maximum number of iterations.
|
||
% Default is 1000.
|
||
%
|
||
% 'maxFinetune' (integer) Maximum number of iterations in
|
||
% fine-tuning. Default 100.
|
||
%
|
||
% 'sampleSize' (number) [0 - 1] Percentage of samples used in
|
||
% one iteration. Samples are chosen in random.
|
||
% Default is 1 (all samples).
|
||
%
|
||
% 'initGuess' (matrix) Initial guess for A. Default is random.
|
||
% You can now do a "one more" like this:
|
||
% [ica, A, W] = fastica(mix, 'numOfIC',3);
|
||
% [ica2, A2, W2] = fastica(mix, 'initGuess', A, 'numOfIC', 4);
|
||
%
|
||
%======================================================================
|
||
% --Graphics and text output:
|
||
%
|
||
% 'verbose' (string) Either 'on' or 'off'. Default is
|
||
% 'on': report progress of algorithm in text format.
|
||
%
|
||
% 'displayMode' (string) Plot running estimates of independent
|
||
% components: 'signals', 'basis', 'filters' or
|
||
% 'off'. Default is 'signals'.
|
||
%
|
||
% 'displayInterval' Number of iterations between plots.
|
||
% Default is 1 (plot after every iteration).
|
||
%
|
||
%======================================================================
|
||
% --Controlling reduction of dimension and whitening:
|
||
%
|
||
% Reduction of dimension is controlled by 'firstEig' and 'lastEig', or
|
||
% alternatively by 'interactivePCA'.
|
||
%
|
||
% 'firstEig' (integer) This and 'lastEig' specify the range for
|
||
% eigenvalues that are retained, 'firstEig' is
|
||
% the index of largest eigenvalue to be
|
||
% retained. Default is 1.
|
||
%
|
||
% 'lastEig' (integer) This is the index of the last (smallest)
|
||
% eigenvalue to be retained. Default equals the
|
||
% dimension of data.
|
||
%
|
||
% 'interactivePCA' (string) Either 'on' or 'off'. When set 'on', the
|
||
% eigenvalues are shown to the user and the
|
||
% range can be specified interactively. Default
|
||
% is 'off'. Can also be set to 'gui'. Then the user
|
||
% can use the same GUI that's in FASTICAG.
|
||
%
|
||
% If you already know the eigenvalue decomposition of the covariance
|
||
% matrix, you can avoid computing it again by giving it with the
|
||
% following options:
|
||
%
|
||
% 'pcaE' (matrix) Eigenvectors
|
||
% 'pcaD' (matrix) Eigenvalues
|
||
%
|
||
% If you already know the whitened data, you can give it directly to
|
||
% the algorithm using the following options:
|
||
%
|
||
% 'whiteSig' (matrix) Whitened signal
|
||
% 'whiteMat' (matrix) Whitening matrix
|
||
% 'dewhiteMat' (matrix) dewhitening matrix
|
||
%
|
||
% If values for all the 'whiteSig', 'whiteSig' and 'dewhiteMat' are
|
||
% suplied, they will be used in computing the ICA. PCA and whitening
|
||
% are not performed. Though 'mixedsig' is not used in the main
|
||
% algorithm it still must be entered - some values are still
|
||
% calculated from it.
|
||
%
|
||
% Performing preprocessing only is possible by the option:
|
||
%
|
||
% 'only' (string) Compute only PCA i.e. reduction of
|
||
% dimension ('pca') or only PCA plus whitening
|
||
% ('white'). Default is 'all': do ICA estimation
|
||
% as well. This option changes the output
|
||
% format accordingly. For example:
|
||
%
|
||
% [whitesig, WM, DWM] = FASTICA(mixedsig,
|
||
% 'only', 'white')
|
||
% returns the whitened signals, the whitening matrix
|
||
% (WM) and the dewhitening matrix (DWM). (See also
|
||
% WHITENV.) In FastICA the whitening matrix performs
|
||
% whitening and the reduction of dimension. Dewhitening
|
||
% matrix is the pseudoinverse of whitening matrix.
|
||
%
|
||
% [E, D] = FASTICA(mixedsig, 'only', 'pca')
|
||
% returns the eigenvector (E) and diagonal
|
||
% eigenvalue (D) matrices containing the
|
||
% selected subspaces.
|
||
%
|
||
%======================================================================
|
||
% EXAMPLES
|
||
%
|
||
% [icasig] = FASTICA (mixedsig, 'approach', 'symm', 'g', 'tanh');
|
||
% Do ICA with tanh nonlinearity and in parallel (like
|
||
% maximum likelihood estimation for supergaussian data).
|
||
%
|
||
% [icasig] = FASTICA (mixedsig, 'lastEig', 10, 'numOfIC', 3);
|
||
% Reduce dimension to 10, and estimate only 3
|
||
% independent components.
|
||
%
|
||
% [icasig] = FASTICA (mixedsig, 'verbose', 'off', 'displayMode', 'off');
|
||
% Don't output convergence reports and don't plot
|
||
% independent components.
|
||
%
|
||
%
|
||
% A graphical user interface for FASTICA can be launched by the
|
||
% command FASTICAG
|
||
%
|
||
% See also FASTICAG
|
||
|
||
% 15.1.2001
|
||
% Hugo G<>vert
|
||
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Remove the mean and check the data
|
||
|
||
[mixedsig, mixedmean] = remmean(mixedsig);
|
||
|
||
[Dim, NumOfSampl] = size(mixedsig);
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Default values for optional parameters
|
||
|
||
% All
|
||
verbose = 'on';
|
||
|
||
% Default values for 'pcamat' parameters
|
||
firstEig = 1;
|
||
lastEig = Dim;
|
||
interactivePCA = 'off';
|
||
|
||
% Default values for 'fpica' parameters
|
||
approach = 'defl';
|
||
numOfIC = Dim;
|
||
g = 'pow3';
|
||
finetune = 'off';
|
||
a1 = 1;
|
||
a2 = 1;
|
||
myy = 1;
|
||
stabilization = 'off';
|
||
epsilon = 0.0001;
|
||
maxNumIterations = 1000;
|
||
maxFinetune = 5;
|
||
initState = 'rand';
|
||
guess = 0;
|
||
sampleSize = 1;
|
||
displayMode = 'signals';
|
||
displayInterval = 1;
|
||
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Parameters for fastICA - i.e. this file
|
||
|
||
b_verbose = 1;
|
||
jumpPCA = 0;
|
||
jumpWhitening = 0;
|
||
only = 3;
|
||
userNumOfIC = 0;
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Parameters passed to icaplot
|
||
% Added by DR, 12.18.2001
|
||
|
||
plottype = 'histogram';
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Read the optional parameters
|
||
|
||
args = varargin;
|
||
if(rem(nargin-1,2)==1)
|
||
error('Optional parameters should always go by pairs');
|
||
else
|
||
a = 1;
|
||
for i=1:(nargin-1)/2
|
||
% get the name and value of parameter
|
||
str_param = args{a};
|
||
val_param = args{a + 1}
|
||
a = a + 2;
|
||
|
||
% change the value of parameter
|
||
if strcmp (str_param, 'stabilization')
|
||
stabilization = val_param;
|
||
elseif strcmp (str_param, 'maxFinetune')
|
||
maxFinetune = val_param;
|
||
elseif strcmp (str_param, 'sampleSize')
|
||
sampleSize = val_param;
|
||
elseif strcmp (str_param, 'verbose')
|
||
verbose = val_param;
|
||
% silence this program also
|
||
if strcmp (verbose, 'off'), b_verbose = 0; end
|
||
elseif strcmp (str_param, 'firstEig')
|
||
firstEig = val_param;
|
||
elseif strcmp (str_param, 'lastEig')
|
||
lastEig = val_param;
|
||
elseif strcmp (str_param, 'interactivePCA')
|
||
interactivePCA = val_param;
|
||
elseif strcmp (str_param, 'approach')
|
||
approach = val_param;
|
||
elseif strcmp (str_param, 'numOfIC')
|
||
numOfIC = val_param;
|
||
% User has suplied new value for numOfIC.
|
||
% We'll use this information later on...
|
||
userNumOfIC = 1;
|
||
elseif strcmp (str_param, 'g')
|
||
g = val_param;
|
||
elseif strcmp (str_param, 'finetune')
|
||
finetune = val_param;
|
||
elseif strcmp (str_param, 'a1')
|
||
a1 = val_param;
|
||
elseif strcmp (str_param, 'a2')
|
||
a2 = val_param;
|
||
elseif strcmp (str_param, 'mu') | strcmp(str_param, 'myy')
|
||
myy = val_param;
|
||
elseif strcmp (str_param, 'epsilon')
|
||
epsilon = val_param;
|
||
elseif strcmp (str_param, 'maxNumIterations')
|
||
maxNumIterations = val_param;
|
||
elseif strcmp (str_param, 'initGuess')
|
||
% no use setting 'guess' if the 'initState' is not set
|
||
initState = 'guess';
|
||
guess = val_param;
|
||
elseif strcmp (str_param, 'displayMode')
|
||
displayMode = val_param;
|
||
elseif strcmp (str_param, 'displayInterval')
|
||
displayInterval = val_param;
|
||
elseif strcmp (str_param, 'pcaE')
|
||
% calculate if there are enought parameters to skip PCA
|
||
jumpPCA = jumpPCA + 1;
|
||
E = val_param;
|
||
elseif strcmp (str_param, 'pcaD')
|
||
% calculate if there are enought parameters to skip PCA
|
||
jumpPCA = jumpPCA + 1;
|
||
D = val_param;
|
||
elseif strcmp (str_param, 'whiteSig')
|
||
% calculate if there are enought parameters to skip PCA and whitening
|
||
jumpWhitening = jumpWhitening + 1;
|
||
whitesig = val_param;
|
||
elseif strcmp (str_param, 'whiteMat')
|
||
% calculate if there are enought parameters to skip PCA and whitening
|
||
jumpWhitening = jumpWhitening + 1;
|
||
whiteningMatrix = val_param;
|
||
elseif strcmp (str_param, 'dewhiteMat')
|
||
% calculate if there are enought parameters to skip PCA and whitening
|
||
jumpWhitening = jumpWhitening + 1;
|
||
dewhiteningMatrix = val_param;
|
||
elseif strcmp (str_param, 'only')
|
||
% if the user only wants to calculate PCA or...
|
||
if strcmp (val_param, 'pca')
|
||
only = 1;
|
||
elseif strcmp(val_param, 'white')
|
||
only = 2;
|
||
elseif strcmp(val_param, 'all')
|
||
only = 3;
|
||
end
|
||
% Added by DR, define the type of plot
|
||
elseif strcmp (str_param, 'plottype')
|
||
plottype = val_param;
|
||
|
||
else
|
||
% Hmmm, something wrong with the parameter string
|
||
error(['Unrecognized parameter: ''' str_param '''']);
|
||
end;
|
||
end;
|
||
end
|
||
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% print information about data
|
||
if b_verbose
|
||
fprintf('Number of signals: %d\n', Dim);
|
||
fprintf('Number of samples: %d\n', NumOfSampl);
|
||
end
|
||
|
||
% Check if the data has been entered the wrong way,
|
||
% but warn only... it may be on purpose
|
||
|
||
if Dim > NumOfSampl
|
||
if b_verbose
|
||
fprintf('Warning: ');
|
||
fprintf('The signal matrix may be oriented in the wrong way.\n');
|
||
fprintf('In that case transpose the matrix.\n\n');
|
||
end
|
||
end
|
||
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Calculating PCA
|
||
|
||
% We need the results of PCA for whitening, but if we don't
|
||
% need to do whitening... then we dont need PCA...
|
||
if jumpWhitening == 3
|
||
if b_verbose,
|
||
fprintf ('Whitened signal and corresponding matrices supplied.\n');
|
||
fprintf ('PCA calculations not needed.\n');
|
||
end;
|
||
else
|
||
|
||
% OK, so first we need to calculate PCA
|
||
% Check to see if we already have the PCA data
|
||
if jumpPCA == 2,
|
||
if b_verbose,
|
||
fprintf ('Values for PCA calculations suplied.\n');
|
||
fprintf ('PCA calculations not needed.\n');
|
||
end;
|
||
else
|
||
% display notice if the user entered one, but not both, of E and D.
|
||
if (jumpPCA > 0) & (b_verbose),
|
||
fprintf ('You must supply all of these in order to jump PCA:\n');
|
||
fprintf ('''pcaE'', ''pcaD''.\n');
|
||
end;
|
||
|
||
% Calculate PCA
|
||
[E, D]=pcamat(mixedsig, firstEig, lastEig, interactivePCA, verbose);
|
||
end
|
||
end
|
||
|
||
% skip the rest if user only wanted PCA
|
||
if only > 1
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Whitening the data
|
||
|
||
% Check to see if the whitening is needed...
|
||
if jumpWhitening == 3,
|
||
if b_verbose,
|
||
fprintf ('Whitening not needed.\n');
|
||
end;
|
||
else
|
||
|
||
% Whitening is needed
|
||
% display notice if the user entered some of the whitening info, but not all.
|
||
if (jumpWhitening > 0) & (b_verbose),
|
||
fprintf ('You must suply all of these in order to jump whitening:\n');
|
||
fprintf ('''whiteSig'', ''whiteMat'', ''dewhiteMat''.\n');
|
||
end;
|
||
|
||
% Calculate the whitening
|
||
[whitesig, whiteningMatrix, dewhiteningMatrix] = whitenv ...
|
||
(mixedsig, E, D, verbose);
|
||
end
|
||
|
||
end % if only > 1
|
||
|
||
% skip the rest if user only wanted PCA and whitening
|
||
if only > 2
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% Calculating the ICA
|
||
|
||
% Check some parameters
|
||
% The dimension of the data may have been reduced during PCA calculations.
|
||
% The original dimension is calculated from the data by default, and the
|
||
% number of IC is by default set to equal that dimension.
|
||
|
||
Dim = size(whitesig, 1);
|
||
|
||
% The number of IC's must be less or equal to the dimension of data
|
||
if numOfIC > Dim
|
||
numOfIC = Dim;
|
||
% Show warning only if verbose = 'on' and user suplied a value for 'numOfIC'
|
||
if (b_verbose & userNumOfIC)
|
||
fprintf('Warning: estimating only %d independent components\n', numOfIC);
|
||
fprintf('(Can''t estimate more independent components than dimension of data)\n');
|
||
end
|
||
end
|
||
|
||
% Calculate the ICA with fixed point algorithm.
|
||
[A, W] = fpica (whitesig, whiteningMatrix, dewhiteningMatrix, approach, ...
|
||
numOfIC, g, finetune, a1, a2, myy, stabilization, epsilon, ...
|
||
maxNumIterations, maxFinetune, initState, guess, sampleSize, ...
|
||
displayMode, displayInterval, verbose, plottype);
|
||
|
||
% Check for valid return
|
||
if ~isempty(W)
|
||
% Add the mean back in.
|
||
if b_verbose
|
||
fprintf('Adding the mean back to the data.\n');
|
||
end
|
||
%icasig = W * mixedsig + (W * mixedmean) * ones(1, NumOfSampl);
|
||
icasig = W * mixedsig;
|
||
if b_verbose & ...
|
||
(max(abs(W * mixedmean)) > 1e-9) & ...
|
||
(strcmp(displayMode,'signals') | strcmp(displayMode,'on'))
|
||
fprintf('Note that the plots don''t have the mean added.\n');
|
||
end
|
||
else
|
||
icasig = [];
|
||
end
|
||
|
||
end % if only > 2
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
% The output depends on the number of output parameters
|
||
% and the 'only' parameter.
|
||
|
||
if only == 1 % only PCA
|
||
Out1 = E;
|
||
Out2 = D;
|
||
elseif only == 2 % only PCA & whitening
|
||
if nargout == 2
|
||
Out1 = whiteningMatrix;
|
||
Out2 = dewhiteningMatrix;
|
||
else
|
||
Out1 = whitesig;
|
||
Out2 = whiteningMatrix;
|
||
Out3 = dewhiteningMatrix;
|
||
end
|
||
else % ICA
|
||
if nargout == 2
|
||
Out1 = A;
|
||
Out2 = W;
|
||
else
|
||
Out1 = icasig;
|
||
Out2 = A;
|
||
Out3 = W;
|
||
end
|
||
end
|
||
|
||
endfunction
|