ext_infomax_runica.m

infomax的ICA算法的扩展程序源码

% runica() - Perform Independent Component Analysis (ICA) decomposition
%            of input data using the logistic infomax ICA algorithm of 
%            Bell & Sejnowski (1995) with the natural gradient feature 
%            of Amari, Cichocki & Yang, or optionally the extended-ICA 
%            algorithm of Lee, Girolami & Sejnowski, with optional PCA 
%            dimension reduction. Annealing based on weight changes is 
%            used to automate the separation process. 
% Usage:
%         >> [weights,sphere] = runica(data); % train using defaults 
%    else
%         >> [weights,sphere,compvars,bias,signs,lrates,activations] ...
%                             = runica(data,'Key1',Value1',...);
% Input:
%    data     = input data (chans,frames*epochs). 
%               Note that if data consists of multiple discontinuous epochs, 
%               each epoch should be separately baseline-zero'd using
%                  >> data = rmbase(data,frames,basevector);
%
% Optional keywords [argument]:
% 'extended'  = [N] perform tanh() "extended-ICA" with sign estimation 
%               N training blocks. If N > 0, automatically estimate the 
%               number of sub-Gaussian sources. If N < 0, fix number of 
%               sub-Gaussian comps to -N [faster than N>0] (default|0 -> off)
% 'pca'       = [N] decompose a principal component     (default -> 0=off)
%               subspace of the data. Value is the number of PCs to retain.
% 'ncomps'    = [N] number of ICA components to compute (default -> chans or 'pca' arg)
%               using rectangular ICA decomposition
% 'sphering'  = ['on'/'off'] flag sphering of data      (default -> 'on')
% 'weights'   = [W] initial weight matrix               (default -> eye())
%                            (Note: if 'sphering' 'off', default -> spher())
% 'lrate'     = [rate] initial ICA learning rate (<< 1) (default -> heuristic)
% 'block'     = [N] ICA block size (<< datalength)      (default -> heuristic)
% 'anneal'    = annealing constant (0,1] (defaults -> 0.90, or 0.98, extended)
%                         controls speed of convergence
% 'annealdeg' = [N] degrees weight change for annealing (default -> 70)
% 'stop'      = [f] stop training when weight-change < this (default -> 1e-6
%               if less than 33 channel and 1E-7 otherwise)
% 'maxsteps'  = [N] max number of ICA training steps    (default -> 512)
% 'bias'      = ['on'/'off'] perform bias adjustment    (default -> 'on')
% 'momentum'  = [0<f<1] training momentum               (default -> 0)
% 'specgram'  = [srate loHz hiHz frames winframes] decompose a complex time/frequency
%               transform of the data (Note: winframes must divide frames) 
%                            (defaults [srate 0 srate/2 size(data,2) size(data,2)])
% 'posact'    = make all component activations net-positive(default 'on'}
% 'verbose'   = give ascii messages ('on'/'off')        (default -> 'on')
%
% Outputs:    [Note: RO means output in reverse order of projected mean variance
%                    unless starting weight matrix passed ('weights' above)]
% weights     = ICA weight matrix (comps,chans)      [RO]
% sphere      = data sphering matrix (chans,chans) = spher(data)
%               Note that unmixing_matrix = weights*sphere {if sphering off -> eye(chans)}
% compvars    = back-projected component variances   [RO]
% bias        = vector of final (ncomps) online bias [RO]    (default = zeros())
% signs       = extended-ICA signs for components    [RO]    (default = ones())
%                   [ -1 = sub-Gaussian; 1 = super-Gaussian]
% lrates      = vector of learning rates used at each training step [RO]
% activations = activation time courses of the output components (ncomps,frames*epochs)
%
% Authors: Scott Makeig with contributions from Tony Bell, Te-Won Lee, 
% Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky, Delorme Arnaud,
% CNL/The Salk Institute, La Jolla, 1996-

% Uses: posact()

% Reference (please cite):
%
% Makeig, S., Bell, A.J., Jung, T-P and Sejnowski, T.J.,
% "Independent component analysis of electroencephalographic data," 
% In: D. Touretzky, M. Mozer and M. Hasselmo (Eds). Advances in Neural 
% Information Processing Systems 8:145-151, MIT Press, Cambridge, MA (1996).
%
% Toolbox Citation:
%
% Makeig, Scott et al. "EEGLAB: ICA Toolbox for Psychophysiological Research". 
% WWW Site, Swartz Center for Computational Neuroscience, Institute of Neural
% Computation, University of San Diego California
% <www.sccn.ucsd.edu/eeglab/>, 2000. [World Wide Web Publication]. 
%
% For more information:
% http://www.sccn.ucsd.edu/eeglab/icafaq.html - FAQ on ICA/EEG
% http://www.sccn.ucsd.edu/eeglab/icabib.html - mss. on ICA & biosignals
% http://www.cnl.salk.edu/~tony/ica.html - math. mss. on ICA

% Copyright (C) 1996 Scott Makeig et al, SCCN/INC/UCSD, [e-mail]
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: runica.m,v $
% Revision 1.20  2004/06/29 16:44:09  scott
% randomize data shuffling by clock state
%
% Revision 1.19  2004/05/16 01:13:19  scott
% typo
%
% Revision 1.18  2004/05/16 01:09:30  scott
% typo
%
% Revision 1.17  2004/05/16 01:07:44  scott
% disable new buggy fprintf (eval(ps)) feature - for emergency use -sm
%
% Revision 1.16  2004/05/09 22:49:45  scott
% made wchange printout hold enough decimal places when 'stop' is small
%
% Revision 1.15  2004/05/09 18:18:36  scott
% added printout of frames per weight -sm
%
% Revision 1.14  2003/12/15 23:28:34  arno
% debug nochangeupdate
%
% Revision 1.13  2003/12/11 17:51:11  arno
% stoping rule debug if more than 32 channels
%
% Revision 1.12  2003/10/23 15:48:45  arno
% indents
%
% Revision 1.11  2003/10/19 17:14:15  scott
% cosmetics, help msg
%
% Revision 1.10  2003/10/03 18:21:25  arno
% releasing constraint pca<nchans-1
%
% Revision 1.9  2003/09/19 01:42:56  arno
% documenting stop at 1E-7 for more than 32 channels
%
% Revision 1.8  2003/09/18 23:43:50  arno
% debug 'weight' input (would not sort component by variance and make some function crash
%
% Revision 1.7  2003/08/19 18:56:14  scott
% added third output arg compvar -sm
%
% Revision 1.6  2003/08/07 18:33:15  arno
% same
%
% Revision 1.5  2003/08/07 18:25:27  arno
% default lrate for more than 32 channels
%
% Revision 1.4  2003/05/23 15:31:53  arno
% adding more details for extended option of runica
%
% Revision 1.3  2003/01/15 22:08:21  arno
% typo
%
% Revision 1.2  2002/10/23 18:09:54  arno
% new interupt button
%
% Revision 1.1  2002/04/05 17:36:45  jorn
% Initial revision
%

%%%%%%%%%%%%%%%%%%%%%%%%%%% Edit history %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  runica()  - by Scott Makeig with contributions from Tony Bell, Te-Won Lee 
%              Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky et al.
%                            CNL / Salk Institute 1996-00
%  04-30-96 built from icatest.m and ~jung/.../wtwpwica.m -sm
%  07-28-97 new runica(), adds bias (default on), momentum (default off), 
%           extended-ICA (Lee & Sejnowski, 1997), cumulative angledelta 
%           (until lrate drops), keywords, signcount for speeding extended-ICA
%  10-07-97 put acos() outside verbose loop; verbose 'off' wasn't stopping -sm
%  11-11-97 adjusted help msg -sm
%  11-30-97 return eye(chans) if sphering 'off' or 'none' (undocumented option) -sm
%  02-27-98 use pinv() instead of inv() to rank order comps if ncomps < chans -sm
%  04-28-98 added 'posact' and 'pca' flags  -sm
%  07-16-98 reduced length of randperm() for kurtosis subset calc. -se & sm
%  07-19-98 fixed typo in weights def. above -tl & sm
%  12-21-99 added 'specgram' option suggested by Michael Zibulevsky, UNM -sm
%  12-22-99 fixed rand() sizing inefficiency on suggestion of Mike Spratling, UK -sm
%  01-11-00 fixed rand() sizing bug on suggestion of Jack Foucher, Strasbourg -sm
%  12-18-00 test for existence of Sig Proc Tlbx function 'specgram'; improve
%           'specgram' option arguments -sm
%  01-25-02 reformated help & license -ad 
%  01-25-02 lowered default lrate and block -ad 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [weights,sphere,meanvar,bias,signs,lrates,activations,y] = runica(data,p1,v1,p2,v2,p3,v3,p4,v4,p5,v5,p6,v6,p7,v7,p8,v8,p9,v9,p10,v10,p11,v11,p12,v12,p13,v13,p14,v14)

if nargin < 1
  help runica  
  return
end

[chans frames] = size(data); % determine the data size
urchans = chans;  % remember original data channels 
datalength = frames;
if chans<2 
   fprintf('\nrunica() - data size (%d,%d) too small.\n\n', chans,frames);
   return
end
%
%%%%%%%%%%%%%%%%%%%%%% Declare defaults used below %%%%%%%%%%%%%%%%%%%%%%%%
%
MAX_WEIGHT           = 1e8;       % guess that weights larger than this have blown up
DEFAULT_STOP         = 0.000001;  % stop training if weight changes below this
DEFAULT_ANNEALDEG    = 60;        % when angle change reaches this value,
DEFAULT_ANNEALSTEP   = 0.90;      %     anneal by multiplying lrate by this
DEFAULT_EXTANNEAL    = 0.98;      %     or this if extended-ICA
DEFAULT_MAXSTEPS     = 512;       % ]top training after this many steps 
DEFAULT_MOMENTUM     = 0.0;       % default momentum weight

DEFAULT_BLOWUP       = 1000000000.0;   % = learning rate has 'blown up'
DEFAULT_BLOWUP_FAC   = 0.8;       % when lrate 'blows up,' anneal by this fac
DEFAULT_RESTART_FAC  = 0.9;       % if weights blowup, restart with lrate
                                  % lower by this factor
MIN_LRATE            = 0.000001;  % if weight blowups make lrate < this, quit
MAX_LRATE            = 0.1;       % guard against uselessly high learning rate
DEFAULT_LRATE        = 0.00065/log(chans); 
                                  % heuristic default - may need adjustment
                                  %   for large or tiny data sets!
% DEFAULT_BLOCK        = floor(sqrt(frames/4));  % heuristic default 
DEFAULT_BLOCK        = min(floor(5*log(frames)),0.3*frames); % heuristic 
                                  % - may need adjustment!
% Extended-ICA option:
DEFAULT_EXTENDED     = 0;         % default off
DEFAULT_EXTBLOCKS    = 1;         % number of blocks per kurtosis calculation
DEFAULT_NSUB         = 1;         % initial default number of assumed sub-Gaussians
                                  % for extended-ICA
DEFAULT_EXTMOMENTUM  = 0.5;       % momentum term for computing extended-ICA kurtosis
MAX_KURTSIZE         = 6000;      % max points to use in kurtosis calculation
MIN_KURTSIZE         = 2000;      % minimum good kurtosis size (flag warning)
SIGNCOUNT_THRESHOLD  = 25;        % raise extblocks when sign vector unchanged
                                  % after this many steps
SIGNCOUNT_STEP       = 2;         % extblocks increment factor 

DEFAULT_SPHEREFLAG   = 'on';      % use the sphere matrix as the default
                                  %   starting weight matrix
DEFAULT_PCAFLAG      = 'off';     % don't use PCA reduction
DEFAULT_POSACTFLAG   = 'on';      % use posact()
DEFAULT_VERBOSE      = 1;         % write ascii info to calling screen
DEFAULT_BIASFLAG     = 1;         % default to using bias in the ICA update rule
%                                 
%%%%%%%%%%%%%%%%%%%%%%% Set up keyword default values %%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargout < 2, 
    fprintf('runica() - needs at least two output arguments.\n');
    return
end
epochs = 1;							 % do not care how many epochs in data

pcaflag    = DEFAULT_PCAFLAG;
sphering   = DEFAULT_SPHEREFLAG;     % default flags
posactflag = DEFAULT_POSACTFLAG;
verbose    = DEFAULT_VERBOSE;

block      = DEFAULT_BLOCK;          % heuristic default - may need adjustment!
lrate      = DEFAULT_LRATE;
annealdeg  = DEFAULT_ANNEALDEG;
annealstep = 0;                      % defaults declared below
nochange   = NaN;
momentum   = DEFAULT_MOMENTUM;
maxsteps   = DEFAULT_MAXSTEPS;

weights    = 0;                      % defaults defined below
ncomps     = chans;
biasflag   = DEFAULT_BIASFLAG;

extended   = DEFAULT_EXTENDED;
extblocks  = DEFAULT_EXTBLOCKS;
kurtsize   = MAX_KURTSIZE;
signsbias  = 0.02;                   % bias towards super-Gaussian components
extmomentum= DEFAULT_EXTMOMENTUM;    % exp. average the kurtosis estimates
nsub       = DEFAULT_NSUB;
wts_blowup = 0;                      % flag =1 when weights too large
wts_passed = 0;                      % flag weights passed as argument
%
%%%%%%%%%% Collect keywords and values from argument list %%%%%%%%%%%%%%%
%
   if (nargin> 1 & rem(nargin,2) == 0)
      fprintf('runica(): Even number of input arguments???')
      return
   end
   for i = 3:2:nargin % for each Keyword
      Keyword = eval(['p',int2str((i-3)/2 +1)]);
      Value = eval(['v',int2str((i-3)/2 +1)]);
      if ~isstr(Keyword)
         fprintf('runica(): keywords must be strings')
         return
      end
      Keyword = lower(Keyword); % convert upper or mixed case to lower

      if strcmp(Keyword,'weights') | strcmp(Keyword,'weight')
         if isstr(Value)
            fprintf(...
      'runica(): weights value must be a weight matrix or sphere')
            return
         else
           weights = Value;
           wts_passed =1;
         end
      elseif strcmp(Keyword,'ncomps')
         if isstr(Value)
            fprintf('runica(): ncomps value must be an integer')
            return
         end
         if ncomps < urchans & ncomps ~= Value
            fprintf('runica(): Use either PCA or ICA dimension reduction');
            return
         end
         ncomps = Value;
         if ~ncomps,
            ncomps = chans;
         end
      elseif strcmp(Keyword,'pca') 
         if ncomps < urchans & ncomps ~= Value
            fprintf('runica(): Use either PCA or ICA dimension reduction');
            return
         end
         if isstr(Value)
            fprintf(...
'runica(): pca value should be the number of principal components to retain')
            return
         end
         pcaflag = 'on';
         ncomps = Value;
         if ncomps > chans | ncomps < 1,
            fprintf('runica(): pca value must be in range [1,%d]\n',chans)
            return
         end
         chans = ncomps;
      elseif strcmp(Keyword,'posact') 
         if ~isstr(Value)
           fprintf('runica(): posact value must be on or off')
           return
         else 
           Value = lower(Value);
           if ~strcmp(Value,'on') & ~strcmp(Value,'off'),
             fprintf('runica(): posact value must be on or off')
             return
           end
           posactflag = Value;
         end
      elseif strcmp(Keyword,'lrate')
         if isstr(Value)
            fprintf('runica(): lrate value must be a number')
            return
         end
         lrate = Value;
         if lrate>MAX_LRATE | lrate <0,
           fprintf('runica(): lrate value is out of bounds'); 
           return
         end
         if ~lrate,
            lrate = DEFAULT_LRATE;
         end
      elseif strcmp(Keyword,'block') | strcmp(Keyword,'blocksize')
         if isstr(Value)
            fprintf('runica(): block size value must be a number')
            return
         end
         block = Value;
         if ~block,
           block = DEFAULT_BLOCK; 
         end
      elseif strcmp(Keyword,'stop') | strcmp(Keyword,'nochange') ...
                    | strcmp(Keyword,'stopping')
         if isstr(Value)
            fprintf('runica(): stop wchange value must be a number')
            return
         end
         nochange = Value;
      elseif strcmp(Keyword,'maxsteps') | strcmp(Keyword,'steps')
         if isstr(Value)
            fprintf('runica(): maxsteps value must be an integer')
            return
         end
         maxsteps = Value;
         if ~maxsteps,