envtopo.m

来自「含有多种ICA算法的eeglab工具箱」· M 代码 · 共 1,112 行 · 第 1/3 页

	end

	%icadefs;    % read toolbox defaults
	ENVCOLORS = strvcat('w..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..');

	if isempty(g.colorfile)
	    g.colorfile = ENVCOLORS; % filename read from icadefs
	end
	if isempty(g.voffsets) | ( size(g.voffsets) == [1,1] & g.voffsets(1) == 0 )
	  g.voffsets = zeros(1,MAXTOPOS); 
	end
	if isempty(g.plotchans) | g.plotchans(1)==0 
	   g.plotchans = 1:chans;
	end
	if max(g.plotchans) > chans | min(g.plotchans) < 1
	   error('invalid ''plotchan'' index');
	end
	if isempty(g.compnums) | g.compnums(1) == 0
	    g.compnums = 1:wtcomps; % by default, all components
	end
	if min(g.compnums) < 0
	  if length(g.compnums) > 1
	     fprintf('envtopo(): negative compnums must be a single integer.\n');
	     return
	  end
	  if -g.compnums > MAXTOPOS
	    fprintf('Can only plot a maximum of %d components.\n',MAXTOPOS);
	    return
	  else
	    MAXTOPOS = -g.compnums;
	    g.compnums = 1:wtcomps;
	  end
	end
	ncomps = length(g.compnums);
	for i=1:ncomps-1
          if g.compnums(i) == 0
	       fprintf('Removing component number 0 in compnums.\n');
	       g.compnums(i)=[];
          elseif g.compnums(i)>wtcomps
	       fprintf('compnums(%d) > number of comps (%d)?\n',i,wtcomps);
	       return
          end
	  for j=i+1:ncomps
	    if g.compnums(i)==g.compnums(j)
	       fprintf('Removing repeated component number (%d) in compnums.\n',g.compnums(i));
	       g.compnums(j)=[];
            end
	  end
	end

	limitset = 0;
	if isempty(g.limits)
	  g.limits = 0;
	end
	if length(g.limits)>1
	  limitset = 1;
	end

	%
	%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	if length(g.limits) < 2
            if g.limits == 0,         
	      xmin = 0;
	      xmax = frames-1;  % use dummy times
            else
	      fprintf('envtopo: limits should be 0, [minms maxms], or [minms maxms minuV maxuV].\n');
              if isempty(g.limits) disp empty
              else g.limits
              end
              return
            end
	     ymin = min(min(data));
	     ymax = max(max(data));
        elseif length(g.limits) == 2 
	     xmin = g.limits(1);
	     xmax = g.limits(2);
	     ymin = min(min(data));
	     ymax = max(max(data));
	     datarange = ymax-ymin;
	     ymin = ymin-0.05*datarange;
	     ymax = ymax+0.05*datarange;
	elseif length(g.limits) == 4
	  xmin = g.limits(1);
	  xmax = g.limits(2);
          if g.limits(3)==0 & g.limits(4)==0 % compute y limits from data
	     ymin = min(min(data));
	     ymax = max(max(data));
	     datarange = ymax-ymin;
	     ymin = ymin-0.05*datarange;
	     ymax = ymax+0.05*datarange;
          else
	     ymin = g.limits(3);
	     ymax = g.limits(4);
          end
        else
	    fprintf('envtopo: limits should be 0, [minms maxms], or [minms maxms minuV maxuV].\n');
            if isempty(g.limits) disp empty
            else g.limits
            end
            return
	end;

	if xmax == 0 & xmin == 0,
	  times = (0:1:frames-1);
	  xmin = 0;
	  xmax = frames-1;
	else
	  dt = (xmax-xmin)/(frames-1);
	  times=xmin*ones(1,frames)+dt*(0:frames-1); % compute times-values
	end;
	if xmax<=xmin,
	  fprintf('envtopo() - maxms must be > minms.\n')
	  return
	end

	dataenvelope = envelope(data, g.envmode);
	if ymax == 0 & ymin == 0,
	  ymax=max(max(dataenvelope));
	  ymin=min(min(dataenvelope));
	  datarange = ymax-ymin;
	  ymin = ymin-0.05*datarange;
	  ymax = ymax+0.05*datarange;
	end
	if ymax<=ymin,
	  fprintf('envtopo() - ymax must be > ymin.\n')
	  return
	end
	%
	%%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	if ~isstr(g.colorfile)
	  fprintf('envproj(): color file name must be a string.\n');
	  return
	end
	if strcmpi(g.colorfile,'bold')
	   all_bold = 1;
	   g.colorfile = ENVCOLORS; % filename read from icadefs
	end
	  if exist(g.colorfile) == 2  % if an existing file
		cid = fopen(g.colorfile,'r');
		if cid <3,
			fprintf('envproj(): cannot open file %s.\n',g.colorfile);
			return
		else
			colors = fscanf(cid,'%s',[3 MAXENVPLOTCHANS]);
			colors = colors';
		end;
	  else
		colors = g.colorfile;
	  end
	  [r c] = size(colors);
	  for i=1:r
		for j=1:c
		  if colors(i,j)=='.',
		    if j==1
			fprintf('envtopo(): colors file should have color letter in 1st column.\n');
			return
		    elseif j==2
			colors(i,j)='-';
		    elseif j>2
			colors(i,j)=' ';
		    end;
		  end;
		end;
	  end;
	colors(1,1) = 'k'; % make sure 1st color (for data envelope) is black
	% [rr cc] = size(colors);
	% colors

	%
	%%%%%%%%%%%%%%% Compute plotframes and envdata %%%%%%%%%%%%%%%%%%%%%
	%

	ntopos = length(g.compnums);
	if ntopos > MAXTOPOS
	  ntopos = MAXTOPOS; % limit the number of topoplots to display
	end

	if max(g.compnums) > wtcomps | min(g.compnums)< 1
	  fprintf(...
	'envtopo(): one or more compnums out of range (1,%d).\n',wtcomps);
	  return
	end

	plotframes = ones(ncomps);
	maxproj = zeros(chans,ncomps);
        %
        % first, plot the data envelope
        %
	envdata = zeros(2,frames*(ncomps+1));
	envdata(:,1:frames) = envelope(data(g.plotchans,:), g.envmode); 
	fprintf('\nComparing projection sizes for components:  ');
        if ncomps>32
           fprintf('\n');
        end
	compvars = zeros(1,ncomps);
        %
	%%%%%%%% find max variances and their frame indices %%%%%%%%%%%
        %
	for c = 1:ncomps 

	  if ~rem(c,5) 
	      fprintf('%d ... ',g.compnums(c)); % c is index into compnums
	  end
	  if ~rem(c,100)
	    fprintf('\n');
	  end
	  %proj = icaproj(data,weights,g.compnums(c)); % updated arg list 12/00 -sm
	  if isempty(g.icaact)
	      proj = g.icawinv(:,g.compnums(c))*weights(g.compnums(c),:)*data; % updated -ad 10/2002
	  else 
	      proj = g.icawinv(:,g.compnums(c))*g.icaact(g.compnums(c),:);     % updated -sm 4/2004
	  end;
	  envdata(:,c*frames+1:(c+1)*frames) = envelope(proj(g.plotchans,:), g.envmode);

  [val,i] = max(sum(proj(:,frame1:frame2).*proj(:,frame1:frame2))); % find max variance
  compvars(c)   = val;

  % find variance in interval after removing component
  if strcmpi(g.pvaf, 'on')
      pvaf(c) = mean(mean((data(:,frame1:frame2)-proj(:,frame1:frame2)).^2)); 
  else
      pvaf(c) = mean(mean(proj(:,frame1:frame2).^2));      
  end;

  i = i+frame1-1;
  if envdata(1,c*frames+i) > ymax % if envelop max at max variance clipped in plot
      ix = find(envdata(1,c*frames+1:(c+1)*frames) > ymax);
      [val,ix] = max(envdata(1,c*frames+ix));
      plotframes(c) = ix; % draw line from max non-clipped env maximum
      maxproj(:,c)  = proj(:,ix);
  else  % draw line from max envelope value at max projection time point
      plotframes(c) = i;
      maxproj(:,c)  = proj(:,i);
  end
end %c
fprintf('\n');

% print percent variance accounted for
% ---------------------------------------
% compute pvaf
fprintf('  in the interval %3.0f to %3.0f ms.\n',1000*times(frame1),1000*times(frame2));
vardat = mean(mean((data(:,frame1:frame2).^2))); % find data variance in interval
if strcmpi(g.pvaf, 'on')
    pvaf = 100-100*pvaf/vardat;
    ot   = 'pvaf';
else 
    pvaf = 100*pvaf/vardat;
    ot   = 'rv';
end;
[sortpvaf spx] = sort(pvaf);
sortpvaf = sortpvaf(end:-1:1);
spx      = spx(end:-1:1);
npercol = ceil(ncomps/3);
% for index =1:npercol
%     try, fprintf('   IC%d\t%s: %6.2f%%\t', spx(index)          , ot, sortpvaf(index)); catch, end;
%     try, fprintf('   IC%d\t%s: %6.2f%%\t', spx(index+npercol)  , ot, sortpvaf(index+npercol)); catch, end;
%     try, fprintf('   IC%d\t%s: %6.2f%%\t', spx(index+2*npercol), ot, sortpvaf(index+2*npercol)); catch, end;
%     fprintf('\n');
% end;

%
%%%%%%%%%%%%%%%%%%%%%%%%% Sort by max variance in data %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
sampint = (xmax-xmin)/(frames-1);     % sampling interval = 1000/srate;
times = xmin:sampint:xmax;            % make vector of times-values

%[compvars,compx] = sort(compvars');  % sort compnums on max variance
[tmp,compx]  = sort(pvaf');   % sort compnums on max variance

compx        = compx(ncomps:-1:1);    % reverse order of sort
compvarorder = g.compnums(compx);     % actual component numbers (output var)
compvars     = compvars(ncomps:-1:1)';% reverse order of sort (output var)
plotframes   = plotframes(compx);     % plotted comps have these max frames 
compframes   = plotframes';           % frame of max variance in each comp (output var)
comptimes    = times(plotframes(compx));  % time of max variance in each comp (output var)
compsplotted = compvarorder(1:ntopos);% (output var)
%
%%%%%%%%%%%%%%%%%%%%%%%% Reduce to ntopos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[plotframes,ifx] = sort(plotframes(1:ntopos));% sort plotframes on their temporal order
plottimes  = times(plotframes);        % convert to times in ms
compx      = compx(ifx);              % indices into compnums, in plotting order
maporder   = g.compnums(compx);       % comp. numbers, in plotting order (l->r)
maxproj    = maxproj(:,compx);        % maps in plotting order 

vlen = length(g.voffsets); % extend voffsets if necessary
while vlen< ntopos
  g.voffsets = [g.voffsets g.voffsets(vlen)]; % repeat last offset given
  vlen=vlen+1;
end

head_sep = 1.2;
topowidth = pos(3)/(ntopos+(ntopos-1)/5); % width of each topoplot
if topowidth > 0.20    % adjust for maximum height
    topowidth = 0.2;
end

if rem(ntopos,2) == 1  % odd number of topos
   topoleft = pos(3)/2 - (floor(ntopos/2)*head_sep + 0.5)*topowidth;
else % even number of topos
   topoleft = pos(3)/2 - (floor(ntopos/2)*head_sep)*topowidth;
end

%
%%%%%%%%%%%%%%%%%%%% Print times and frames of comp maxes %%%%%%%%%%%%%%
%

fprintf('\n');
fprintf('Plotting envelopes of %d component projections.\n',ntopos);
if length(g.plotchans) ~= chans
  fprintf('Envelopes computed from %d specified data channels.\n',...
      length(g.plotchans));
end
fprintf('Topo maps will show components: ');
for t=1:ntopos
  fprintf('%4d  ',maporder(t));
end
fprintf('\n');
fprintf('    with max variance at times: ');
for t=1:ntopos
  fprintf('%4.0f  ',1000*plottimes(t));
end
fprintf('\n');

fprintf('               or epoch frames: ');
for t=1:ntopos
   fprintf('%4d  ',frame1-1+plotframes(t));
end
fprintf('\n');

if strcmp(g.pvaf,'on')
   fprintf('    component pvaf in interval:  ');
   for t=1:ntopos
     fprintf('%4.2f ',pvaf(t));
   end
   fprintf('\n');
end

sumproj = zeros(size(data));
for n = 1:ntopos
  if isempty(g.icaact)
      sumproj = sumproj + g.icawinv(:,maporder(n))*weights(maporder(n),:)*data; % updated -ad 10/2002
  else 
      sumproj = sumproj + g.icawinv(:,maporder(n))*g.icaact(maporder(n),:);     % updated -sm 4/2004
  end;
end
varproj = mean(mean((data(g.plotchans,frame1:frame2).^2))); % find data variance in interval
if strcmpi(g.pvaf, 'on')
      sumpvaf = mean(mean((data(g.plotchans,frame1:frame2)-sumproj(g.plotchans,frame1:frame2)).^2)); 
  else
      sumpvaf = mean(mean(sumproj(g.plotchans,frame1:frame2).^2));      
  end;
if strcmpi(g.pvaf, 'on')
    sumpvaf = 100-100*sumpvaf/varproj;
    ot   = 'pvaf';
else 
    sumpvaf = 100*sumpvaf/varproj;
    ot   = 'rv';
end;
fprintf('    summed component %s in interval: %4.2f %%\n',ot,sumpvaf);
%
%%%%%%%%%%%%%%%%%%%%% Plot the data envelopes %%%%%%%%%%%%%%%%%%%%%%%%%
%
BACKCOLOR = [0.7 0.7 0.7];
newaxes=axes('position',pos);
axis off
%set(newaxes,'Units','Normalized','Position',...
%           [0 0 1 1],'FontSize',16,'FontWeight','Bold','Visible','off');
set(newaxes,'FontSize',16,'FontWeight','Bold','Visible','off');