ext_infomax_runica.m

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

            maxsteps   = DEFAULT_MAXSTEPS;
         end
         if maxsteps < 0
            fprintf('runica(): maxsteps value (%d) must be a positive integer',maxsteps)
            return
         end
      elseif strcmp(Keyword,'anneal') | strcmp(Keyword,'annealstep')
         if isstr(Value)
            fprintf('runica(): anneal step value (%2.4f) must be a number (0,1)',Value)
            return
         end
         annealstep = Value;
         if annealstep <=0 | annealstep > 1,
            fprintf('runica(): anneal step value (%2.4f) must be (0,1]',annealstep)
            return
         end
      elseif strcmp(Keyword,'annealdeg') | strcmp(Keyword,'degrees')
         if isstr(Value)
            fprintf('runica(): annealdeg value must be a number')
            return
         end
         annealdeg = Value;
         if ~annealdeg,
             annealdeg = DEFAULT_ANNEALDEG;
         elseif annealdeg > 180 | annealdeg < 0
          fprintf('runica(): annealdeg (%3.1f) is out of bounds [0,180]',...
                annealdeg);
          return
                                              
         end
      elseif strcmp(Keyword,'momentum')
         if isstr(Value)
            fprintf('runica(): momentum value must be a number')
            return
         end
         momentum = Value;
         if momentum > 1.0 | momentum < 0
          fprintf('runica(): momentum value is out of bounds [0,1]')
          return
         end
      elseif strcmp(Keyword,'sphering') | strcmp(Keyword,'sphereing') ...
                | strcmp(Keyword,'sphere')
         if ~isstr(Value)
           fprintf('runica(): sphering value must be on, off, or none')
           return
         else 
           Value = lower(Value);
           if ~strcmp(Value,'on') & ~strcmp(Value,'off') & ~strcmp(Value,'none'),
             fprintf('runica(): sphering value must be on or off')
             return
           end
           sphering = Value;
         end
      elseif strcmp(Keyword,'bias')
         if ~isstr(Value)
           fprintf('runica(): bias value must be on or off')
           return
         else 
           Value = lower(Value);
           if strcmp(Value,'on') 
              biasflag = 1;
           elseif strcmp(Value,'off'),
              biasflag = 0;
           else
              fprintf('runica(): bias value must be on or off')
              return
           end
         end
      elseif strcmp(Keyword,'specgram') | strcmp(Keyword,'spec')

         if ~exist('specgram') < 2 % if ~exist or defined workspace variable
           fprintf(...
   'runica(): MATLAB Sig. Proc. Toolbox function "specgram" not found.\n')
           return
         end
         if isstr(Value)
           fprintf('runica(): specgram argument must be a vector')
           return
         end
         srate = Value(1);
         if (srate < 0)
             fprintf('runica(): specgram srate (%4.1f) must be >=0',srate)
             return
           end
         if length(Value)>1
           loHz = Value(2);
           if (loHz < 0 | loHz > srate/2)
             fprintf('runica(): specgram loHz must be >=0 and <= srate/2 (%4.1f)',srate/2)
             return
           end
         else
           loHz = 0; % default
         end
         if length(Value)>2
           hiHz = Value(3);
           if (hiHz < loHz | hiHz > srate/2)
             fprintf('runica(): specgram hiHz must be >=loHz (%4.1f) and <= srate/2 (%4.1f)',loHz,srate/2)
             return
           end
         else
           hiHz = srate/2; % default
         end
         if length(Value)>3
           Hzframes = Value(5);
           if (Hzframes<0 | Hzframes > size(data,2))
             fprintf('runica(): specgram frames must be >=0 and <= data length (%d)',size(data,2))
             return
           end
         else
           Hzframes = size(data,2); % default
         end
         if length(Value)>4
           Hzwinlen = Value(4);
           if rem(Hzframes,Hzwinlen) % if winlen doesn't divide frames
             fprintf('runica(): specgram Hzinc must divide frames (%d)',Hzframes)
             return
           end
         else
           Hzwinlen = Hzframes; % default
         end
         Specgramflag = 1; % set flag to perform specgram()

      elseif strcmp(Keyword,'extended') | strcmp(Keyword,'extend')
         if isstr(Value)
           fprintf('runica(): extended value must be an integer (+/-)')
           return
         else
           extended = 1;      % turn on extended-ICA
           extblocks = fix(Value); % number of blocks per kurt() compute
           if extblocks < 0
                nsub = -1*fix(extblocks);  % fix this many sub-Gauss comps
           elseif ~extblocks,
                extended = 0;             % turn extended-ICA off
           elseif kurtsize>frames,   % length of kurtosis calculation
                kurtsize = frames;
                if kurtsize < MIN_KURTSIZE
                   fprintf(...
   'runica() warning: kurtosis values inexact for << %d points.\n',...
                                                         MIN_KURTSIZE);
                end
           end
         end
      elseif strcmp(Keyword,'verbose') 
         if ~isstr(Value)
            fprintf('runica(): verbose flag value must be on or off')
            return
         elseif strcmp(Value,'on'),
             verbose = 1; 
         elseif strcmp(Value,'off'),
             verbose = 0; 
         else
             fprintf('runica(): verbose flag value must be on or off')
             return
         end
      else
         fprintf('runica(): unknown flag')
         return
      end
   end

%
%%%%%%%%%%%%%%%%%%%%%%%% Initialize weights, etc. %%%%%%%%%%%%%%%%%%%%%%%%
%
if ~annealstep,
  if ~extended,
    annealstep = DEFAULT_ANNEALSTEP;     % defaults defined above
  else
    annealstep = DEFAULT_EXTANNEAL;       % defaults defined above
  end
end % else use annealstep from commandline

if ~annealdeg, 
    annealdeg  = DEFAULT_ANNEALDEG - momentum*90; % heuristic
    if annealdeg < 0,
        annealdeg = 0;
    end
end
if ncomps >  chans | ncomps < 1
    fprintf('runica(): number of components must be 1 to %d.\n',chans);
    return
end

if weights ~= 0,                    % initialize weights
  % starting weights are being passed to runica() from the commandline
    if verbose,
       fprintf('Using starting weight matrix named in argument list ...\n')
    end
    if  chans>ncomps & weights ~=0,
        [r,c]=size(weights);
        if r~=ncomps | c~=chans,
     fprintf(...
      'runica(): weight matrix must have %d rows, %d columns.\n', ...
                              chans,ncomps);
            return;
        end
    end
end;   
%
%%%%%%%%%%%%%%%%%%%%% Check keyword values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if frames<chans,
    fprintf('runica(): data length (%d) < data channels (%d)!\n',frames,chans)
    return
elseif block < 2,
    fprintf('runica(): block size %d too small!\n',block)
    return
elseif block > frames, 
    fprintf('runica(): block size exceeds data length!\n');
    return
elseif floor(epochs) ~= epochs,
    fprintf('runica(): data length is not a multiple of the epoch length!\n');
    return
elseif nsub > ncomps
    fprintf('runica(): there can be at most %d sub-Gaussian components!\n',ncomps);
    return
end;

% 
% adjust nochange if necessary
%
if isnan(nochange) 
    if ncomps > 32
        nochange = 1E-7;
        nochangeupdated = 1; % for fprinting purposes
    else
        nochangeupdated = 1; % for fprinting purposes
        nochange = DEFAULT_STOP;
    end;
else 
    nochangeupdated = 0;
end;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Process the data %%%%%%%%%%%%%%%%%%%%%%%%%%
%
if verbose,
   fprintf( ...
            '\nInput data size [%d,%d] = %d channels, %d frames/n', ...
                        chans,frames,chans,frames);
   if strcmp(pcaflag,'on')
      fprintf('After PCA dimension reduction,\n  finding ');
   else
      fprintf('Finding ');
   end
   if ~extended
      fprintf('%d ICA components using logistic ICA.\n',ncomps);
   else % if extended
      fprintf('%d ICA components using extended ICA.\n',ncomps);
      if extblocks > 0
         fprintf(...
 'Kurtosis will be calculated initially every %d blocks using %d data points.\n',...
                                            extblocks,     kurtsize);
      else
         fprintf(...
 'Kurtosis will not be calculated. Exactly %d sub-Gaussian components assumed.\n',...
                                        nsub);
      end
   end
   fprintf('Decomposing %d frames per ICA weight ((%d)^2 = %d weights, %d frames)\n',...
           floor(frames/ncomps.^2),ncomps.^2,frames);
   fprintf('Initial learning rate will be %g, block size %d.\n',...
           lrate,block);
   if momentum>0,
      fprintf('Momentum will be %g.\n',momentum);
   end
   fprintf( ...
'Learning rate will be multiplied by %g whenever angledelta >= %g deg.\n', ...
               annealstep,annealdeg);

   if nochangeupdated 
       fprintf('More than 32 channels: default stopping weight change 1E-7\n');
   end;
   fprintf('Training will end when wchange < %g or after %d steps.\n', ...
                    nochange,maxsteps);
  if biasflag,
    fprintf('Online bias adjustment will be used.\n');
  else
    fprintf('Online bias adjustment will not be used.\n');
  end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%% Remove overall row means %%%%%%%%%%%%%%%%%%%%%%%%
%
if verbose,
    fprintf('Removing mean of each channel ...\n');
end
data = data - mean(data')'*ones(1,frames);      % subtract row means

if verbose,
   fprintf('Final training data range: %g to %g\n', ...
                          min(min(data)),max(max(data)));
end
%
%%%%%%%%%%%%%%%%%%% Perform PCA reduction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmp(pcaflag,'on')
   fprintf('Reducing the data to %d principal dimensions...\n',ncomps);
   [eigenvectors,eigenvalues,data] = pcsquash(data,ncomps);
   % make data its projection onto the ncomps-dim principal subspace
end
%
%%%%%%%%%%%%%%%%%%% Perform specgram transformation %%%%%%%%%%%%%%%%%%%%%%%
%
if exist('Specgramflag') == 1
  % [P F T] = SPECGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) % MATLAB Sig Proc Toolbox
  % Hzwinlen =  fix(srate/Hzinc); % CHANGED FROM THIS 12/18/00 -sm

  Hzfftlen = 2^(ceil(log(Hzwinlen)/log(2)));   % make FFT length next higher 2^k
  Hzoverlap = 0; % use sequential windows
  %
  % Get freqs and times from 1st channel analysis
  %
  [tmp,freqs,tms] = specgram(data(1,:),Hzfftlen,srate,Hzwinlen,Hzoverlap);

  fs = find(freqs>=loHz & freqs <= hiHz);
  if isempty(fs)
    fprintf('runica(): specified frequency range too narrow!\n');
    return
  end;
    
  specdata = reshape(tmp(fs,:),1,length(fs)*size(tmp,2));
  specdata = [real(specdata) imag(specdata)];
     % fprintf('   size(fs) = %d,%d\n',size(fs,1),size(fs,2));
     % fprintf('   size(tmp) = %d,%d\n',size(tmp,1),size(tmp,2));
  %
  % Loop through remaining channels
  %
  for ch=2:chans
    [tmp] = specgram(data(ch,:),Hzwinlen,srate,Hzwinlen,Hzoverlap);
      tmp = reshape((tmp(fs,:)),1,length(fs)*size(tmp,2));
    specdata = [specdata;[real(tmp) imag(tmp)]]; % channels are rows
  end
  %
  % Print specgram confirmation and details
  %
  fprintf(...
 'Converted data to %d channels by %d=2*%dx%d points spectrogram data.\n',...
                    chans,2*length(fs)*length(tms),length(fs),length(tms));
  if length(fs) > 1
    fprintf(...
 '   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);
  else
    fprintf(...
 '   Low Hz %g, high Hz %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),Hzwinlen);
  end
  %
  % Replace data with specdata
  %
  data = specdata;
  datalength=size(data,2);
end
%
%%%%%%%%%%%%%%%%%%% Perform sphering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if strcmp(sphering,'on'), %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if verbose,
      fprintf('Computing the sphering matrix...\n');
  end
  sphere = 2.0*inv(sqrtm(cov(data'))); % find the "sphering" matrix = spher()
  if ~weights,
      if verbose,
          fprintf('Starting weights are the identity matrix ...\n');
      end
      weights = eye(ncomps,chans); % begin with the identity matrix
  else % weights given on commandline
      if verbose,
          fprintf('Using starting weights named on commandline ...\n');
      end
  end
  if verbose,
      fprintf('Sphering the data ...\n');
  end
  data = sphere*data;      % actually decorrelate the electrode signals

elseif strcmp(sphering,'off') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ~weights
      if verbose,
       fprintf('Using the sphering matrix as the starting weight matrix ...\n');