Description of runica

0001 % runica() - Perform Independent Component Analysis (ICA) decomposition
0002 %            of input data using the logistic infomax ICA algorithm of
0003 %            Bell & Sejnowski (1995) with the natural gradient feature
0004 %            of Amari, Cichocki & Yang, or optionally the extended-ICA
0005 %            algorithm of Lee, Girolami & Sejnowski, with optional PCA
0006 %            dimension reduction. Annealing based on weight changes is
0007 %            used to automate the separation process.
0008 % Usage:
0009 %         >> [weights,sphere] = runica(data); % train using defaults
0010 %    else
0011 %         >> [weights,sphere,compvars,bias,signs,lrates,activations] ...
0012 %                             = runica(data,'Key1',Value1',...);
0013 % Input:
0014 %    data     = input data (chans,frames*epochs).
0015 %               Note that if data consists of multiple discontinuous epochs,
0016 %               each epoch should be separately baseline-zero'd using
0017 %                  >> data = rmbase(data,frames,basevector);
0018 %
0019 % Optional keywords [argument]:
0020 % 'extended'  = [N] perform tanh() "extended-ICA" with sign estimation
0021 %               N training blocks. If N > 0, automatically estimate the
0022 %               number of sub-Gaussian sources. If N < 0, fix number of
0023 %               sub-Gaussian comps to -N [faster than N>0] (default|0 -> off)
0024 % 'pca'       = [N] decompose a principal component     (default -> 0=off)
0025 %               subspace of the data. Value is the number of PCs to retain.
0026 % 'ncomps'    = [N] number of ICA components to compute (default -> chans or 'pca' arg)
0027 %               using rectangular ICA decomposition
0028 % 'sphering'  = ['on'/'off'] flag sphering of data      (default -> 'on')
0029 % 'weights'   = [W] initial weight matrix               (default -> eye())
0030 %                            (Note: if 'sphering' 'off', default -> spher())
0031 % 'lrate'     = [rate] initial ICA learning rate (<< 1) (default -> heuristic)
0032 % 'block'     = [N] ICA block size (<< datalength)      (default -> heuristic)
0033 % 'anneal'    = annealing constant (0,1] (defaults -> 0.90, or 0.98, extended)
0034 %                         controls speed of convergence
0035 % 'annealdeg' = [N] degrees weight change for annealing (default -> 70)
0036 % 'stop'      = [f] stop training when weight-change < this (default -> 1e-6
0037 %               if less than 33 channel and 1E-7 otherwise)
0038 % 'maxsteps'  = [N] max number of ICA training steps    (default -> 512)
0039 % 'bias'      = ['on'/'off'] perform bias adjustment    (default -> 'on')
0040 % 'momentum'  = [0<f<1] training momentum               (default -> 0)
0041 % 'specgram'  = [srate loHz hiHz frames winframes] decompose a complex time/frequency
0042 %               transform of the data (Note: winframes must divide frames)
0043 %                            (defaults [srate 0 srate/2 size(data,2) size(data,2)])
0044 % 'posact'    = make all component activations net-positive(default 'on'}
0045 % 'verbose'   = give ascii messages ('on'/'off')        (default -> 'on')
0046 %
0047 % Outputs:    [Note: RO means output in reverse order of projected mean variance
0048 %                    unless starting weight matrix passed ('weights' above)]
0049 % weights     = ICA weight matrix (comps,chans)      [RO]
0050 % sphere      = data sphering matrix (chans,chans) = spher(data)
0051 %               Note that unmixing_matrix = weights*sphere {if sphering off -> eye(chans)}
0052 % compvars    = back-projected component variances   [RO]
0053 % bias        = vector of final (ncomps) online bias [RO]    (default = zeros())
0054 % signs       = extended-ICA signs for components    [RO]    (default = ones())
0055 %                   [ -1 = sub-Gaussian; 1 = super-Gaussian]
0056 % lrates      = vector of learning rates used at each training step [RO]
0057 % activations = activation time courses of the output components (ncomps,frames*epochs)
0058 %
0059 % Authors: Scott Makeig with contributions from Tony Bell, Te-Won Lee,
0060 % Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky, Delorme Arnaud,
0061 % CNL/The Salk Institute, La Jolla, 1996-
0062 
0063 % Uses: posact()
0064 
0065 % Reference (please cite):
0066 %
0067 % Makeig, S., Bell, A.J., Jung, T-P and Sejnowski, T.J.,
0068 % "Independent component analysis of electroencephalographic data,"
0069 % In: D. Touretzky, M. Mozer and M. Hasselmo (Eds). Advances in Neural
0070 % Information Processing Systems 8:145-151, MIT Press, Cambridge, MA (1996).
0071 %
0072 % Toolbox Citation:
0073 %
0074 % Makeig, Scott et al. "EEGLAB: ICA Toolbox for Psychophysiological Research".
0075 % WWW Site, Swartz Center for Computational Neuroscience, Institute of Neural
0076 % Computation, University of San Diego California
0077 % <www.sccn.ucsd.edu/eeglab/>, 2000. [World Wide Web Publication].
0078 %
0079 % For more information:
0080 % http://www.sccn.ucsd.edu/eeglab/icafaq.html - FAQ on ICA/EEG
0081 % http://www.sccn.ucsd.edu/eeglab/icabib.html - mss. on ICA & biosignals
0082 % http://www.cnl.salk.edu/~tony/ica.html - math. mss. on ICA
0083 
0084 % Copyright (C) 1996 Scott Makeig et al, SCCN/INC/UCSD, scott@sccn.ucsd.edu
0085 %
0086 % This program is free software; you can redistribute it and/or modify
0087 % it under the terms of the GNU General Public License as published by
0088 % the Free Software Foundation; either version 2 of the License, or
0089 % (at your option) any later version.
0090 %
0091 % This program is distributed in the hope that it will be useful,
0092 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0093 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0094 % GNU General Public License for more details.
0095 %
0096 % You should have received a copy of the GNU General Public License
0097 % along with this program; if not, write to the Free Software
0098 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
0099 
0100 % $Log: runica.m,v $
0101 % Revision 1.23  2004/10/06 20:41:10  scott
0102 % rm debug -sm
0103 %
0104 % Revision 1.22  2004/10/06 20:40:32  scott
0105 % same -sm
0106 %
0107 % Revision 1.21  2004/10/06 20:39:20  scott
0108 % fixed block size to integer -sm
0109 %
0110 % Revision 1.20  2004/06/29 16:44:09  scott
0111 % randomize data shuffling by clock state
0112 %
0113 % Revision 1.19  2004/05/16 01:13:19  scott
0114 % typo
0115 %
0116 % Revision 1.18  2004/05/16 01:09:30  scott
0117 % typo
0118 %
0119 % Revision 1.17  2004/05/16 01:07:44  scott
0120 % disable new buggy fprintf (eval(ps)) feature - for emergency use -sm
0121 %
0122 % Revision 1.16  2004/05/09 22:49:45  scott
0123 % made wchange printout hold enough decimal places when 'stop' is small
0124 %
0125 % Revision 1.15  2004/05/09 18:18:36  scott
0126 % added printout of frames per weight -sm
0127 %
0128 % Revision 1.14  2003/12/15 23:28:34  arno
0129 % debug nochangeupdate
0130 %
0131 % Revision 1.13  2003/12/11 17:51:11  arno
0132 % stoping rule debug if more than 32 channels
0133 %
0134 % Revision 1.12  2003/10/23 15:48:45  arno
0135 % indents
0136 %
0137 % Revision 1.11  2003/10/19 17:14:15  scott
0138 % cosmetics, help msg
0139 %
0140 % Revision 1.10  2003/10/03 18:21:25  arno
0141 % releasing constraint pca<nchans-1
0142 %
0143 % Revision 1.9  2003/09/19 01:42:56  arno
0144 % documenting stop at 1E-7 for more than 32 channels
0145 %
0146 % Revision 1.8  2003/09/18 23:43:50  arno
0147 % debug 'weight' input (would not sort component by variance and make some function crash
0148 %
0149 % Revision 1.7  2003/08/19 18:56:14  scott
0150 % added third output arg compvar -sm
0151 %
0152 % Revision 1.6  2003/08/07 18:33:15  arno
0153 % same
0154 %
0155 % Revision 1.5  2003/08/07 18:25:27  arno
0156 % default lrate for more than 32 channels
0157 %
0158 % Revision 1.4  2003/05/23 15:31:53  arno
0159 % adding more details for extended option of runica
0160 %
0161 % Revision 1.3  2003/01/15 22:08:21  arno
0162 % typo
0163 %
0164 % Revision 1.2  2002/10/23 18:09:54  arno
0165 % new interupt button
0166 %
0167 % Revision 1.1  2002/04/05 17:36:45  jorn
0168 % Initial revision
0169 %
0170 
0171 %%%%%%%%%%%%%%%%%%%%%%%%%%% Edit history %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0172 %
0173 %  runica()  - by Scott Makeig with contributions from Tony Bell, Te-Won Lee
0174 %              Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky et al.
0175 %                            CNL / Salk Institute 1996-00
0176 %  04-30-96 built from icatest.m and ~jung/.../wtwpwica.m -sm
0177 %  07-28-97 new runica(), adds bias (default on), momentum (default off),
0178 %           extended-ICA (Lee & Sejnowski, 1997), cumulative angledelta
0179 %           (until lrate drops), keywords, signcount for speeding extended-ICA
0180 %  10-07-97 put acos() outside verbose loop; verbose 'off' wasn't stopping -sm
0181 %  11-11-97 adjusted help msg -sm
0182 %  11-30-97 return eye(chans) if sphering 'off' or 'none' (undocumented option) -sm
0183 %  02-27-98 use pinv() instead of inv() to rank order comps if ncomps < chans -sm
0184 %  04-28-98 added 'posact' and 'pca' flags  -sm
0185 %  07-16-98 reduced length of randperm() for kurtosis subset calc. -se & sm
0186 %  07-19-98 fixed typo in weights def. above -tl & sm
0187 %  12-21-99 added 'specgram' option suggested by Michael Zibulevsky, UNM -sm
0188 %  12-22-99 fixed rand() sizing inefficiency on suggestion of Mike Spratling, UK -sm
0189 %  01-11-00 fixed rand() sizing bug on suggestion of Jack Foucher, Strasbourg -sm
0190 %  12-18-00 test for existence of Sig Proc Tlbx function 'specgram'; improve
0191 %           'specgram' option arguments -sm
0192 %  01-25-02 reformated help & license -ad
0193 %  01-25-02 lowered default lrate and block -ad
0194 %
0195 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0196 
0197 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)
0198 
0199 if nargin < 1
0200   help runica  
0201   return
0202 end
0203 
0204 [chans frames] = size(data); % determine the data size
0205 urchans = chans;  % remember original data channels
0206 datalength = frames;
0207 if chans<2 
0208    fprintf('\nrunica() - data size (%d,%d) too small.\n\n', chans,frames);
0209    return
0210 end
0211 %
0212 %%%%%%%%%%%%%%%%%%%%%% Declare defaults used below %%%%%%%%%%%%%%%%%%%%%%%%
0213 %
0214 MAX_WEIGHT           = 1e8;       % guess that weights larger than this have blown up
0215 DEFAULT_STOP         = 0.000001;  % stop training if weight changes below this
0216 DEFAULT_ANNEALDEG    = 60;        % when angle change reaches this value,
0217 DEFAULT_ANNEALSTEP   = 0.90;      %     anneal by multiplying lrate by this
0218 DEFAULT_EXTANNEAL    = 0.98;      %     or this if extended-ICA
0219 DEFAULT_MAXSTEPS     = 512;       % ]top training after this many steps
0220 DEFAULT_MOMENTUM     = 0.0;       % default momentum weight
0221 
0222 DEFAULT_BLOWUP       = 1000000000.0;   % = learning rate has 'blown up'
0223 DEFAULT_BLOWUP_FAC   = 0.8;       % when lrate 'blows up,' anneal by this fac
0224 DEFAULT_RESTART_FAC  = 0.9;       % if weights blowup, restart with lrate
0225                                   % lower by this factor
0226 MIN_LRATE            = 0.000001;  % if weight blowups make lrate < this, quit
0227 MAX_LRATE            = 0.1;       % guard against uselessly high learning rate
0228 DEFAULT_LRATE        = 0.00065/log(chans); 
0229                                   % heuristic default - may need adjustment
0230                                   %   for large or tiny data sets!
0231 % DEFAULT_BLOCK        = floor(sqrt(frames/4));  % heuristic default
0232 DEFAULT_BLOCK          = ceil(min(5*log(frames),0.3*frames)); % heuristic
0233                                   % - may need adjustment!
0234 % Extended-ICA option:
0235 DEFAULT_EXTENDED     = 0;         % default off
0236 DEFAULT_EXTBLOCKS    = 1;         % number of blocks per kurtosis calculation
0237 DEFAULT_NSUB         = 1;         % initial default number of assumed sub-Gaussians
0238                                   % for extended-ICA
0239 DEFAULT_EXTMOMENTUM  = 0.5;       % momentum term for computing extended-ICA kurtosis
0240 MAX_KURTSIZE         = 6000;      % max points to use in kurtosis calculation
0241 MIN_KURTSIZE         = 2000;      % minimum good kurtosis size (flag warning)
0242 SIGNCOUNT_THRESHOLD  = 25;        % raise extblocks when sign vector unchanged
0243                                   % after this many steps
0244 SIGNCOUNT_STEP       = 2;         % extblocks increment factor
0245 
0246 DEFAULT_SPHEREFLAG   = 'on';      % use the sphere matrix as the default
0247                                   %   starting weight matrix
0248 DEFAULT_PCAFLAG      = 'off';     % don't use PCA reduction
0249 DEFAULT_POSACTFLAG   = 'on';      % use posact()
0250 DEFAULT_VERBOSE      = 1;         % write ascii info to calling screen
0251 DEFAULT_BIASFLAG     = 1;         % default to using bias in the ICA update rule
0252 %
0253 %%%%%%%%%%%%%%%%%%%%%%% Set up keyword default values %%%%%%%%%%%%%%%%%%%%%%%%%
0254 %
0255 if nargout < 2, 
0256     fprintf('runica() - needs at least two output arguments.\n');
0257     return
0258 end
0259 epochs = 1;                             % do not care how many epochs in data
0260 
0261 pcaflag    = DEFAULT_PCAFLAG;
0262 sphering   = DEFAULT_SPHEREFLAG;     % default flags
0263 posactflag = DEFAULT_POSACTFLAG;
0264 verbose    = DEFAULT_VERBOSE;
0265 
0266 block      = DEFAULT_BLOCK;          % heuristic default - may need adjustment!
0267 lrate      = DEFAULT_LRATE;
0268 annealdeg  = DEFAULT_ANNEALDEG;
0269 annealstep = 0;                      % defaults declared below
0270 nochange   = NaN;
0271 momentum   = DEFAULT_MOMENTUM;
0272 maxsteps   = DEFAULT_MAXSTEPS;
0273 
0274 weights    = 0;                      % defaults defined below
0275 ncomps     = chans;
0276 biasflag   = DEFAULT_BIASFLAG;
0277 
0278 extended   = DEFAULT_EXTENDED;
0279 extblocks  = DEFAULT_EXTBLOCKS;
0280 kurtsize   = MAX_KURTSIZE;
0281 signsbias  = 0.02;                   % bias towards super-Gaussian components
0282 extmomentum= DEFAULT_EXTMOMENTUM;    % exp. average the kurtosis estimates
0283 nsub       = DEFAULT_NSUB;
0284 wts_blowup = 0;                      % flag =1 when weights too large
0285 wts_passed = 0;                      % flag weights passed as argument
0286 %
0287 %%%%%%%%%% Collect keywords and values from argument list %%%%%%%%%%%%%%%
0288 %
0289    if (nargin> 1 & rem(nargin,2) == 0)
0290       fprintf('runica(): Even number of input arguments???')
0291       return
0292    end
0293    for i = 3:2:nargin % for each Keyword
0294       Keyword = eval(['p',int2str((i-3)/2 +1)]);
0295       Value = eval(['v',int2str((i-3)/2 +1)]);
0296       if ~isstr(Keyword)
0297          fprintf('runica(): keywords must be strings')
0298          return
0299       end
0300       Keyword = lower(Keyword); % convert upper or mixed case to lower
0301 
0302       if strcmp(Keyword,'weights') | strcmp(Keyword,'weight')
0303          if isstr(Value)
0304             fprintf(...
0305       'runica(): weights value must be a weight matrix or sphere')
0306             return
0307          else
0308            weights = Value;
0309            wts_passed =1;
0310          end
0311       elseif strcmp(Keyword,'ncomps')
0312          if isstr(Value)
0313             fprintf('runica(): ncomps value must be an integer')
0314             return
0315          end
0316          if ncomps < urchans & ncomps ~= Value
0317             fprintf('runica(): Use either PCA or ICA dimension reduction');
0318             return
0319          end
0320          ncomps = Value;
0321          if ~ncomps,
0322             ncomps = chans;
0323          end
0324       elseif strcmp(Keyword,'pca') 
0325          if ncomps < urchans & ncomps ~= Value
0326             fprintf('runica(): Use either PCA or ICA dimension reduction');
0327             return
0328          end
0329          if isstr(Value)
0330             fprintf(...
0331 'runica(): pca value should be the number of principal components to retain')
0332             return
0333          end
0334          pcaflag = 'on';
0335          ncomps = Value;
0336          if ncomps > chans | ncomps < 1,
0337             fprintf('runica(): pca value must be in range [1,%d]\n',chans)
0338             return
0339          end
0340          chans = ncomps;
0341       elseif strcmp(Keyword,'posact') 
0342          if ~isstr(Value)
0343            fprintf('runica(): posact value must be on or off')
0344            return
0345          else 
0346            Value = lower(Value);
0347            if ~strcmp(Value,'on') & ~strcmp(Value,'off'),
0348              fprintf('runica(): posact value must be on or off')
0349              return
0350            end
0351            posactflag = Value;
0352          end
0353       elseif strcmp(Keyword,'lrate')
0354          if isstr(Value)
0355             fprintf('runica(): lrate value must be a number')
0356             return
0357          end
0358          lrate = Value;
0359          if lrate>MAX_LRATE | lrate <0,
0360            fprintf('runica(): lrate value is out of bounds'); 
0361            return
0362          end
0363          if ~lrate,
0364             lrate = DEFAULT_LRATE;
0365          end
0366       elseif strcmp(Keyword,'block') | strcmp(Keyword,'blocksize')
0367          if isstr(Value)
0368             fprintf('runica(): block size value must be a number')
0369             return
0370          end
0371          block = floor(Value);
0372          if ~block,
0373            block = DEFAULT_BLOCK; 
0374          end
0375       elseif strcmp(Keyword,'stop') | strcmp(Keyword,'nochange') ...
0376                     | strcmp(Keyword,'stopping')
0377          if isstr(Value)
0378             fprintf('runica(): stop wchange value must be a number')
0379             return
0380          end
0381          nochange = Value;
0382       elseif strcmp(Keyword,'maxsteps') | strcmp(Keyword,'steps')
0383          if isstr(Value)
0384             fprintf('runica(): maxsteps value must be an integer')
0385             return
0386          end
0387          maxsteps = Value;
0388          if ~maxsteps,
0389             maxsteps   = DEFAULT_MAXSTEPS;
0390          end
0391          if maxsteps < 0
0392             fprintf('runica(): maxsteps value (%d) must be a positive integer',maxsteps)
0393             return
0394          end
0395       elseif strcmp(Keyword,'anneal') | strcmp(Keyword,'annealstep')
0396          if isstr(Value)
0397             fprintf('runica(): anneal step value (%2.4f) must be a number (0,1)',Value)
0398             return
0399          end
0400          annealstep = Value;
0401          if annealstep <=0 | annealstep > 1,
0402             fprintf('runica(): anneal step value (%2.4f) must be (0,1]',annealstep)
0403             return
0404          end
0405       elseif strcmp(Keyword,'annealdeg') | strcmp(Keyword,'degrees')
0406          if isstr(Value)
0407             fprintf('runica(): annealdeg value must be a number')
0408             return
0409          end
0410          annealdeg = Value;
0411          if ~annealdeg,
0412              annealdeg = DEFAULT_ANNEALDEG;
0413          elseif annealdeg > 180 | annealdeg < 0
0414           fprintf('runica(): annealdeg (%3.1f) is out of bounds [0,180]',...
0415                 annealdeg);
0416           return
0417                                               
0418          end
0419       elseif strcmp(Keyword,'momentum')
0420          if isstr(Value)
0421             fprintf('runica(): momentum value must be a number')
0422             return
0423          end
0424          momentum = Value;
0425          if momentum > 1.0 | momentum < 0
0426           fprintf('runica(): momentum value is out of bounds [0,1]')
0427           return
0428          end
0429       elseif strcmp(Keyword,'sphering') | strcmp(Keyword,'sphereing') ...
0430                 | strcmp(Keyword,'sphere')
0431          if ~isstr(Value)
0432            fprintf('runica(): sphering value must be on, off, or none')
0433            return
0434          else 
0435            Value = lower(Value);
0436            if ~strcmp(Value,'on') & ~strcmp(Value,'off') & ~strcmp(Value,'none'),
0437              fprintf('runica(): sphering value must be on or off')
0438              return
0439            end
0440            sphering = Value;
0441          end
0442       elseif strcmp(Keyword,'bias')
0443          if ~isstr(Value)
0444            fprintf('runica(): bias value must be on or off')
0445            return
0446          else 
0447            Value = lower(Value);
0448            if strcmp(Value,'on') 
0449               biasflag = 1;
0450            elseif strcmp(Value,'off'),
0451               biasflag = 0;
0452            else
0453               fprintf('runica(): bias value must be on or off')
0454               return
0455            end
0456          end
0457       elseif strcmp(Keyword,'specgram') | strcmp(Keyword,'spec')
0458 
0459          if ~exist('specgram') < 2 % if ~exist or defined workspace variable
0460            fprintf(...
0461    'runica(): MATLAB Sig. Proc. Toolbox function "specgram" not found.\n')
0462            return
0463          end
0464          if isstr(Value)
0465            fprintf('runica(): specgram argument must be a vector')
0466            return
0467          end
0468          srate = Value(1);
0469          if (srate < 0)
0470              fprintf('runica(): specgram srate (%4.1f) must be >=0',srate)
0471              return
0472            end
0473          if length(Value)>1
0474            loHz = Value(2);
0475            if (loHz < 0 | loHz > srate/2)
0476              fprintf('runica(): specgram loHz must be >=0 and <= srate/2 (%4.1f)',srate/2)
0477              return
0478            end
0479          else
0480            loHz = 0; % default
0481          end
0482          if length(Value)>2
0483            hiHz = Value(3);
0484            if (hiHz < loHz | hiHz > srate/2)
0485              fprintf('runica(): specgram hiHz must be >=loHz (%4.1f) and <= srate/2 (%4.1f)',loHz,srate/2)
0486              return
0487            end
0488          else
0489            hiHz = srate/2; % default
0490          end
0491          if length(Value)>3
0492            Hzframes = Value(5);
0493            if (Hzframes<0 | Hzframes > size(data,2))
0494              fprintf('runica(): specgram frames must be >=0 and <= data length (%d)',size(data,2))
0495              return
0496            end
0497          else
0498            Hzframes = size(data,2); % default
0499          end
0500          if length(Value)>4
0501            Hzwinlen = Value(4);
0502            if rem(Hzframes,Hzwinlen) % if winlen doesn't divide frames
0503              fprintf('runica(): specgram Hzinc must divide frames (%d)',Hzframes)
0504              return
0505            end
0506          else
0507            Hzwinlen = Hzframes; % default
0508          end
0509          Specgramflag = 1; % set flag to perform specgram()
0510 
0511       elseif strcmp(Keyword,'extended') | strcmp(Keyword,'extend')
0512          if isstr(Value)
0513            fprintf('runica(): extended value must be an integer (+/-)')
0514            return
0515          else
0516            extended = 1;      % turn on extended-ICA
0517            extblocks = fix(Value); % number of blocks per kurt() compute
0518            if extblocks < 0
0519                 nsub = -1*fix(extblocks);  % fix this many sub-Gauss comps
0520            elseif ~extblocks,
0521                 extended = 0;             % turn extended-ICA off
0522            elseif kurtsize>frames,   % length of kurtosis calculation
0523                 kurtsize = frames;
0524                 if kurtsize < MIN_KURTSIZE
0525                    fprintf(...
0526    'runica() warning: kurtosis values inexact for << %d points.\n',...
0527                                                          MIN_KURTSIZE);
0528                 end
0529            end
0530          end
0531       elseif strcmp(Keyword,'verbose') 
0532          if ~isstr(Value)
0533             fprintf('runica(): verbose flag value must be on or off')
0534             return
0535          elseif strcmp(Value,'on'),
0536              verbose = 1; 
0537          elseif strcmp(Value,'off'),
0538              verbose = 0; 
0539          else
0540              fprintf('runica(): verbose flag value must be on or off')
0541              return
0542          end
0543       else
0544          fprintf('runica(): unknown flag')
0545          return
0546       end
0547    end
0548 
0549 %
0550 %%%%%%%%%%%%%%%%%%%%%%%% Initialize weights, etc. %%%%%%%%%%%%%%%%%%%%%%%%
0551 %
0552 if ~annealstep,
0553   if ~extended,
0554     annealstep = DEFAULT_ANNEALSTEP;     % defaults defined above
0555   else
0556     annealstep = DEFAULT_EXTANNEAL;       % defaults defined above
0557   end
0558 end % else use annealstep from commandline
0559 
0560 if ~annealdeg, 
0561     annealdeg  = DEFAULT_ANNEALDEG - momentum*90; % heuristic
0562     if annealdeg < 0,
0563         annealdeg = 0;
0564     end
0565 end
0566 if ncomps >  chans | ncomps < 1
0567     fprintf('runica(): number of components must be 1 to %d.\n',chans);
0568     return
0569 end
0570 
0571 if weights ~= 0,                    % initialize weights
0572   % starting weights are being passed to runica() from the commandline
0573     if verbose,
0574        fprintf('Using starting weight matrix named in argument list ...\n')
0575     end
0576     if  chans>ncomps & weights ~=0,
0577         [r,c]=size(weights);
0578         if r~=ncomps | c~=chans,
0579      fprintf(...
0580       'runica(): weight matrix must have %d rows, %d columns.\n', ...
0581                               chans,ncomps);
0582             return;
0583         end
0584     end
0585 end;   
0586 %
0587 %%%%%%%%%%%%%%%%%%%%% Check keyword values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0588 %
0589 if frames<chans,
0590     fprintf('runica(): data length (%d) < data channels (%d)!\n',frames,chans)
0591     return
0592 elseif block < 2,
0593     fprintf('runica(): block size %d too small!\n',block)
0594     return
0595 elseif block > frames, 
0596     fprintf('runica(): block size exceeds data length!\n');
0597     return
0598 elseif floor(epochs) ~= epochs,
0599     fprintf('runica(): data length is not a multiple of the epoch length!\n');
0600     return
0601 elseif nsub > ncomps
0602     fprintf('runica(): there can be at most %d sub-Gaussian components!\n',ncomps);
0603     return
0604 end;
0605 
0606 %
0607 % adjust nochange if necessary
0608 %
0609 if isnan(nochange) 
0610     if ncomps > 32
0611         nochange = 1E-7;
0612         nochangeupdated = 1; % for fprinting purposes
0613     else
0614         nochangeupdated = 1; % for fprinting purposes
0615         nochange = DEFAULT_STOP;
0616     end;
0617 else 
0618     nochangeupdated = 0;
0619 end;
0620 
0621 %
0622 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Process the data %%%%%%%%%%%%%%%%%%%%%%%%%%
0623 %
0624 if verbose,
0625    fprintf( ...
0626             '\nInput data size [%d,%d] = %d channels, %d frames/n', ...
0627                         chans,frames,chans,frames);
0628    if strcmp(pcaflag,'on')
0629       fprintf('After PCA dimension reduction,\n  finding ');
0630    else
0631       fprintf('Finding ');
0632    end
0633    if ~extended
0634       fprintf('%d ICA components using logistic ICA.\n',ncomps);
0635    else % if extended
0636       fprintf('%d ICA components using extended ICA.\n',ncomps);
0637       if extblocks > 0
0638          fprintf(...
0639  'Kurtosis will be calculated initially every %d blocks using %d data points.\n',...
0640                                             extblocks,     kurtsize);
0641       else
0642          fprintf(...
0643  'Kurtosis will not be calculated. Exactly %d sub-Gaussian components assumed.\n',...
0644                                         nsub);
0645       end
0646    end
0647    fprintf('Decomposing %d frames per ICA weight ((%d)^2 = %d weights, %d frames)\n',...
0648            floor(frames/ncomps.^2),ncomps.^2,frames);
0649    fprintf('Initial learning rate will be %g, block size %d.\n',...
0650            lrate,block);
0651    if momentum>0,
0652       fprintf('Momentum will be %g.\n',momentum);
0653    end
0654    fprintf( ...
0655 'Learning rate will be multiplied by %g whenever angledelta >= %g deg.\n', ...
0656                annealstep,annealdeg);
0657 
0658    if nochangeupdated 
0659        fprintf('More than 32 channels: default stopping weight change 1E-7\n');
0660    end;
0661    fprintf('Training will end when wchange < %g or after %d steps.\n', ...
0662                     nochange,maxsteps);
0663   if biasflag,
0664     fprintf('Online bias adjustment will be used.\n');
0665   else
0666     fprintf('Online bias adjustment will not be used.\n');
0667   end
0668 end
0669 %
0670 %%%%%%%%%%%%%%%%%%%%%%%%% Remove overall row means %%%%%%%%%%%%%%%%%%%%%%%%
0671 %
0672 if verbose,
0673     fprintf('Removing mean of each channel ...\n');
0674 end
0675 data = data - mean(data')'*ones(1,frames);      % subtract row means
0676 
0677 if verbose,
0678    fprintf('Final training data range: %g to %g\n', ...
0679                           min(min(data)),max(max(data)));
0680 end
0681 %
0682 %%%%%%%%%%%%%%%%%%% Perform PCA reduction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0683 %
0684 if strcmp(pcaflag,'on')
0685    fprintf('Reducing the data to %d principal dimensions...\n',ncomps);
0686    [eigenvectors,eigenvalues,data] = pcsquash(data,ncomps);
0687    % make data its projection onto the ncomps-dim principal subspace
0688 end
0689 %
0690 %%%%%%%%%%%%%%%%%%% Perform specgram transformation %%%%%%%%%%%%%%%%%%%%%%%
0691 %
0692 if exist('Specgramflag') == 1
0693   % [P F T] = SPECGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) % MATLAB Sig Proc Toolbox
0694   % Hzwinlen =  fix(srate/Hzinc); % CHANGED FROM THIS 12/18/00 -sm
0695 
0696   Hzfftlen = 2^(ceil(log(Hzwinlen)/log(2)));   % make FFT length next higher 2^k
0697   Hzoverlap = 0; % use sequential windows
0698   %
0699   % Get freqs and times from 1st channel analysis
0700   %
0701   [tmp,freqs,tms] = specgram(data(1,:),Hzfftlen,srate,Hzwinlen,Hzoverlap);
0702 
0703   fs = find(freqs>=loHz & freqs <= hiHz);
0704   if isempty(fs)
0705     fprintf('runica(): specified frequency range too narrow!\n');
0706     return
0707   end;
0708     
0709   specdata = reshape(tmp(fs,:),1,length(fs)*size(tmp,2));
0710   specdata = [real(specdata) imag(specdata)];
0711      % fprintf('   size(fs) = %d,%d\n',size(fs,1),size(fs,2));
0712      % fprintf('   size(tmp) = %d,%d\n',size(tmp,1),size(tmp,2));
0713   %
0714   % Loop through remaining channels
0715   %
0716   for ch=2:chans
0717     [tmp] = specgram(data(ch,:),Hzwinlen,srate,Hzwinlen,Hzoverlap);
0718       tmp = reshape((tmp(fs,:)),1,length(fs)*size(tmp,2));
0719     specdata = [specdata;[real(tmp) imag(tmp)]]; % channels are rows
0720   end
0721   %
0722   % Print specgram confirmation and details
0723   %
0724   fprintf(...
0725  'Converted data to %d channels by %d=2*%dx%d points spectrogram data.\n',...
0726                     chans,2*length(fs)*length(tms),length(fs),length(tms));
0727   if length(fs) > 1
0728     fprintf(...
0729  '   Low Hz %g, high Hz %g, Hz incr %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),freqs(fs(2))-freqs(fs(1)),Hzwinlen);
0730   else
0731     fprintf(...
0732  '   Low Hz %g, high Hz %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),Hzwinlen);
0733   end
0734   %
0735   % Replace data with specdata
0736   %
0737   data = specdata;
0738   datalength=size(data,2);
0739 end
0740 %
0741 %%%%%%%%%%%%%%%%%%% Perform sphering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0742 %
0743 
0744 if strcmp(sphering,'on'), %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0745   if verbose,
0746       fprintf('Computing the sphering matrix...\n');
0747   end
0748   sphere = 2.0*inv(sqrtm(cov(data'))); % find the "sphering" matrix = spher()
0749   if ~weights,
0750       if verbose,
0751           fprintf('Starting weights are the identity matrix ...\n');
0752       end
0753       weights = eye(ncomps,chans); % begin with the identity matrix
0754   else % weights given on commandline
0755       if verbose,
0756           fprintf('Using starting weights named on commandline ...\n');
0757       end
0758   end
0759   if verbose,
0760       fprintf('Sphering the data ...\n');
0761   end
0762   data = sphere*data;      % actually decorrelate the electrode signals
0763 
0764 elseif strcmp(sphering,'off') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0765   if ~weights
0766       if verbose,
0767        fprintf('Using the sphering matrix as the starting weight matrix ...\n');
0768        fprintf('Returning the identity matrix in variable "sphere" ...\n');
0769       end
0770       sphere = 2.0*inv(sqrtm(cov(data'))); % find the "sphering" matrix = spher()
0771       weights = eye(ncomps,chans)*sphere; % begin with the identity matrix
0772       sphere = eye(chans);                 % return the identity matrix
0773   else % weights ~= 0
0774       if verbose,
0775        fprintf('Using starting weights named on commandline ...\n');
0776        fprintf('Returning the identity matrix in variable "sphere" ...\n');
0777       end
0778       sphere = eye(chans);                 % return the identity matrix
0779   end
0780 elseif strcmp(sphering,'none')
0781   sphere = eye(chans);                     % return the identity matrix
0782   if ~weights
0783       if verbose,
0784        fprintf('Starting weights are the identity matrix ...\n');
0785        fprintf('Returning the identity matrix in variable "sphere" ...\n');
0786       end
0787       weights = eye(ncomps,chans); % begin with the identity matrix
0788   else % weights ~= 0
0789       if verbose,
0790        fprintf('Using starting weights named on commandline ...\n');
0791        fprintf('Returning the identity matrix in variable "sphere" ...\n');
0792       end
0793   end
0794   sphere = eye(chans,chans);
0795   if verbose,
0796       fprintf('Returned variable "sphere" will be the identity matrix.\n');
0797   end
0798 end
0799 %
0800 %%%%%%%%%%%%%%%%%%%%%%%% Initialize ICA training %%%%%%%%%%%%%%%%%%%%%%%%%
0801 %
0802   lastt=fix((datalength/block-1)*block+1);
0803   BI=block*eye(ncomps,ncomps);
0804   delta=zeros(1,chans*ncomps);
0805   changes = [];
0806   degconst = 180./pi;
0807   startweights = weights;
0808   prevweights = startweights;
0809   oldweights = startweights;
0810   prevwtchange = zeros(chans,ncomps);
0811   oldwtchange = zeros(chans,ncomps);
0812   lrates = zeros(1,maxsteps); 
0813   onesrow = ones(1,block); 
0814   bias = zeros(ncomps,1);
0815   signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
0816   for k=1:nsub
0817       signs(k) = -1;
0818   end
0819   if extended & extblocks < 0 & verbose,
0820     fprintf('Fixed extended-ICA sign assignments:  ');
0821     for k=1:ncomps
0822        fprintf('%d ',signs(k));    
0823     end; fprintf('\n');
0824   end
0825   signs = diag(signs); % make a diagonal matrix
0826   oldsigns = zeros(size(signs));;
0827   signcount = 0;              % counter for same-signs
0828   signcounts = [];
0829   urextblocks = extblocks;    % original value, for resets
0830   old_kk = zeros(1,ncomps);   % for kurtosis momemtum
0831 %
0832 %%%%%%%% ICA training loop using the logistic sigmoid %%%%%%%%%%%%%%%%%%%
0833 %
0834   if verbose,
0835     fprintf('Beginning ICA training ...');
0836     if extended,
0837         fprintf(' first training step may be slow ...\n');
0838     else
0839         fprintf('\n');
0840     end
0841   end
0842   step=0;
0843   laststep=0; 
0844   blockno = 1;  % running block counter for kurtosis interrupts
0845 
0846   rand('state',sum(100*clock));  % set the random number generator state to
0847                                  % a position dependent on the system clock
0848   while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0849       permute=randperm(datalength); % shuffle data order at each step
0850 
0851       for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%%
0852           pause(0);
0853           if ~isempty(get(0, 'currentfigure')) & strcmp(get(gcf, 'tag'), 'stop')
0854               close; error('USER ABORT');
0855           end;
0856           if biasflag                                                   
0857               u=weights*data(:,permute(t:t+block-1)) + bias*onesrow;      
0858           else                                                             
0859               u=weights*data(:,permute(t:t+block-1));                      
0860           end                                                              
0861           if ~extended
0862               %%%%%%%%%%%%%%%%%%% Logistic ICA weight update %%%%%%%%%%%%%%%%%%%
0863               y=1./(1+exp(-u));                                                %
0864               weights = weights + lrate*(BI+(1-2*y)*u')*weights;               %
0865               %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0866           else % Tanh extended-ICA weight update
0867                %%%%%%%%%%%%%%%%%%% Extended-ICA weight update %%%%%%%%%%%%%%%%%%%
0868                y=tanh(u);                                                       %
0869                weights = weights + lrate*(BI-signs*y*u'-u*u')*weights;          %
0870                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0871           end
0872           if biasflag 
0873               if ~extended
0874                   %%%%%%%%%%%%%%%%%%%%%%%% Logistic ICA bias %%%%%%%%%%%%%%%%%%%%%%%
0875                   bias = bias + lrate*sum((1-2*y)')'; % for logistic nonlin. %
0876                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0877               else % extended
0878                    %%%%%%%%%%%%%%%%%%% Extended-ICA bias %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0879                    bias = bias + lrate*sum((-2*y)')';  % for tanh() nonlin.   %
0880                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0881               end                                    
0882           end
0883 
0884       if momentum > 0 %%%%%%%%% Add momentum %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0885         weights = weights + momentum*prevwtchange;                
0886         prevwtchange = weights-prevweights;                      
0887         prevweights = weights;                                  
0888       end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0889 
0890       if max(max(abs(weights))) > MAX_WEIGHT
0891           wts_blowup = 1;
0892           change = nochange;
0893       end
0894       if extended & ~wts_blowup
0895        %
0896        %%%%%%%%%%% Extended-ICA kurtosis estimation %%%%%%%%%%%%%%%%%%%%%
0897        %
0898        if extblocks > 0 & rem(blockno,extblocks) == 0, 
0899                                   % recompute signs vector using kurtosis
0900         if kurtsize < frames % 12-22-99 rand() size suggestion by M. Spratling
0901            rp = fix(rand(1,kurtsize)*datalength);  % pick random subset
0902            % Accout for the possibility of a 0 generation by rand
0903            ou = find(rp == 0);
0904            while ~isempty(ou) % 1-11-00 suggestion by J. Foucher
0905               rp(ou) = fix(rand(1,length(ou))*datalength);
0906               ou = find(rp == 0);
0907            end
0908            partact=weights*data(:,rp(1:kurtsize));
0909         else                                        % for small data sets,
0910             partact=weights*data;                   % use whole data
0911         end
0912         m2=mean(partact'.^2).^2; 
0913         m4= mean(partact'.^4);
0914         kk= (m4./m2)-3.0;                           % kurtosis estimates
0915         if extmomentum
0916             kk = extmomentum*old_kk + (1.0-extmomentum)*kk; % use momentum
0917             old_kk = kk;
0918         end
0919         signs=diag(sign(kk+signsbias));             % pick component signs
0920         if signs == oldsigns,
0921            signcount = signcount+1;
0922         else
0923            signcount = 0;
0924         end
0925         oldsigns = signs;
0926         signcounts = [signcounts signcount];
0927         if signcount >= SIGNCOUNT_THRESHOLD,
0928          extblocks = fix(extblocks * SIGNCOUNT_STEP);% make kurt() estimation
0929          signcount = 0;                             % less frequent if sign
0930         end                                         % is not changing
0931        end % extblocks > 0 & . . .
0932       end % if extended %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0933       blockno = blockno + 1;
0934       if wts_blowup
0935          break
0936       end
0937     end % training block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0938 
0939     if ~wts_blowup
0940       oldwtchange = weights-oldweights;
0941       step=step+1; 
0942       %
0943       %%%%%%% Compute and print weight and update angle changes %%%%%%%%%
0944       %
0945       lrates(1,step) = lrate;
0946       angledelta=0.;
0947       delta=reshape(oldwtchange,1,chans*ncomps);
0948       change=delta*delta'; 
0949     end
0950     %
0951     %%%%%%%%%%%%%%%%%%%%%% Restart if weights blow up %%%%%%%%%%%%%%%%%%%%
0952     %
0953     if wts_blowup | isnan(change)|isinf(change),  % if weights blow up,
0954         fprintf('');
0955         step = 0;                          % start again
0956         change = nochange;
0957         wts_blowup = 0;                    % re-initialize variables
0958         blockno = 1;
0959         lrate = lrate*DEFAULT_RESTART_FAC; % with lower learning rate
0960         weights = startweights;            % and original weight matrix
0961         oldweights = startweights;            
0962         change = nochange;
0963         oldwtchange = zeros(chans,ncomps);
0964         delta=zeros(1,chans*ncomps);
0965         olddelta = delta;
0966         extblocks = urextblocks;
0967         prevweights = startweights;
0968         prevwtchange = zeros(chans,ncomps);
0969         lrates = zeros(1,maxsteps);
0970         bias = zeros(ncomps,1);
0971         if extended
0972            signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
0973            for k=1:nsub
0974                signs(k) = -1;
0975            end
0976            signs = diag(signs); % make a diagonal matrix
0977            oldsigns = zeros(size(signs));;
0978         end
0979         if lrate> MIN_LRATE
0980           r = rank(data);
0981           if r<ncomps
0982            fprintf('Data has rank %d. Cannot compute %d components.\n',...
0983                                   r,ncomps);
0984            return
0985           else
0986            fprintf(...
0987             'Lowering learning rate to %g and starting again.\n',lrate);
0988           end
0989         else
0990           fprintf( ...
0991           'runica(): QUITTING - weight matrix may not be invertible!\n');
0992           return;
0993         end
0994     else % if weights in bounds
0995      %
0996      %%%%%%%%%%%%% Print weight update information %%%%%%%%%%%%%%%%%%%%%%
0997      %
0998      if step> 2 
0999          angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
1000      end
1001      if verbose,
1002       places = -floor(log10(nochange));
1003       if step > 2, 
1004           if ~extended,
1005            ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, angledelta %4.1f deg\n', ...
1006                        step,      lrate,        places+1,places,   degconst*angledelta);
1007           else
1008            ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, angledelta %4.1f deg, %d subgauss\n',...
1009                        step,      lrate,        degconst*angledelta,...
1010                                                      places+1,places,           (ncomps-sum(diag(signs)))/2);
1011           end
1012       elseif ~extended
1013          ps = sprintf('step %d - lrate %5f, wchange %%%d.%df\n',...
1014                      step,      lrate,       places+1,places  );
1015       else
1016          ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, %d subgauss\n',...
1017                      step,      lrate,        places+1,places, (ncomps-sum(diag(signs)))/2);
1018       end % step > 2
1019       fprintf('step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ...
1020                        step,      lrate,     change, degconst*angledelta);
1021       % fprintf(ps,change);  % <---- BUG !!!!
1022      end; % if verbose
1023   %
1024   %%%%%%%%%%%%%%%%%%%% Save current values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1025   %
1026       changes = [changes change];
1027       oldweights = weights;
1028   %
1029   %%%%%%%%%%%%%%%%%%%% Anneal learning rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1030   %
1031       if degconst*angledelta > annealdeg,  
1032         lrate = lrate*annealstep;          % anneal learning rate
1033         olddelta   = delta;                % accumulate angledelta until
1034         oldchange  = change;               %  annealdeg is reached
1035       elseif step == 1                     % on first step only
1036         olddelta   = delta;                % initialize
1037         oldchange  = change;               
1038       end
1039   %
1040   %%%%%%%%%%%%%%%%%%%% Apply stopping rule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1041   %
1042       if step >2 & change < nochange,      % apply stopping rule
1043         laststep=step;            
1044            step=maxsteps;                  % stop when weights stabilize
1045       elseif change > DEFAULT_BLOWUP,      % if weights blow up,
1046         lrate=lrate*DEFAULT_BLOWUP_FAC;    % keep trying
1047       end;                                 % with a smaller learning rate
1048     end; % end if weights in bounds
1049 
1050   end; % end training %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1051 
1052   if ~laststep
1053     laststep = step;
1054   end;
1055   lrates = lrates(1,1:laststep);           % truncate lrate history vector
1056   %
1057   %%%%%%%%%%%%%% Orient components towards max positive activation %%%%%%
1058   %
1059   if strcmp(posactflag,'on')
1060       [activations,winvout,weights] = posact(data,weights);
1061        % changes signs of activations and weights to make activations
1062        % net rms-positive
1063   else
1064        activations = weights*data;
1065   end
1066   %
1067   %%%%%%%%%%%%%% If pcaflag, compose PCA and ICA matrices %%%%%%%%%%%%%%%
1068   %
1069   if strcmp(pcaflag,'on')
1070     fprintf('Composing the eigenvector, weights, and sphere matrices\n');
1071     fprintf('  into a single rectangular weights matrix; sphere=eye(%d)\n'...
1072                                                                   ,chans);
1073     weights= weights*sphere*eigenvectors(:,1:ncomps)'; 
1074     sphere = eye(urchans);
1075   end
1076   %
1077   %%%%%% Sort components in descending order of max projected variance %%%%
1078   %
1079   if verbose,
1080     fprintf(...
1081    'Sorting components in descending order of mean projected variance ...\n');
1082   end
1083   %
1084   %%%%%%%%%%%%%%%%%%%% Find mean variances %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1085   %
1086   meanvar  = zeros(ncomps,1);      % size of the projections
1087   if ncomps == urchans % if weights are square . . .
1088       winv = inv(weights*sphere);
1089   else
1090       fprintf('Using pseudo-inverse of weight matrix to rank order component projections.\n');
1091       winv = pinv(weights*sphere);
1092   end
1093   for s=1:ncomps
1094       if verbose,
1095           fprintf('%d ',s);         % construct single-component data matrix
1096       end
1097       % project to scalp, then add row means
1098       compproj = winv(:,s)*activations(s,:);
1099       meanvar(s) = mean(sum(compproj.*compproj)/(size(compproj,1)-1));
1100       % compute mean variance
1101   end                                   % at all scalp channels
1102   if verbose,
1103       fprintf('\n');
1104   end
1105   %
1106   %%%%%%%%%%%%%% Sort components by mean variance %%%%%%%%%%%%%%%%%%%%%%%%
1107   %
1108   [sortvar, windex] = sort(meanvar);
1109   windex = windex(ncomps:-1:1); % order large to small
1110   meanvar = meanvar(windex);
1111   %
1112   %%%%%%%%%%%%%%%%%%%%% Filter data using final weights %%%%%%%%%%%%%%%%%%
1113   %
1114   if nargout>6, % if activations are to be returned
1115       if verbose,
1116           fprintf('Permuting the activation wave forms ...\n');
1117       end
1118       activations = activations(windex,:);
1119   else
1120       clear activations
1121   end
1122   weights = weights(windex,:);% reorder the weight matrix
1123   bias  = bias(windex);       % reorder them
1124   signs = diag(signs);        % vectorize the signs matrix
1125   signs = signs(windex);      % reorder them
1126 
1127 %
1128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1129 %
1130     
1131 return
1132 
1133 %
1134 %%%%%%%%%%%%%%%%%% return nonlinearly-transformed data  %%%%%%%%%%%%%%%%
1135 %
1136 if nargout > 7
1137   u=weights*data + bias*ones(1,frames);      
1138   y = zeros(size(u));
1139   for c=1:chans
1140     for f=1:frames
1141        y(c,f) = 1/(1+exp(-u(c,f)));
1142     end
1143   end
1144 end
runica

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
