erpimage.m

含有多种ICA算法的eeglab工具箱

%
%%%%%%%%%%%%%%%%%%%%%% Remove the ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(Rmerp, 'yes')
    data = data - repmat(nan_mean(data')', [1 size(data,2)]);
end;

%
%%%%%%%%%%%%%%% Sort the data trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if exist('phargs') == 1 % if phase-sort
	if length(phargs) >= 4 % find max frequency in specified band
        if exist('psd') == 2
            [pxx,freqs] = psd(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
        else
            [pxx,freqs] = spec(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
        end;
	%gf = gcf; % figure;plot(freqs,pxx); %xx=axis; %axis([phargs(3) phargs(4) xx(3) xx(4)]); %figure(gf);
		
		pxx = 10*log10(pxx);
		n = find(freqs >= phargs(3) & freqs <= phargs(4));
		if ~length(n)
			freq = phargs(3);
		end
		[dummy maxx] = max(pxx(n));
		freq = freqs(n(maxx));
	else
		freq = phargs(3); % else use specified frequency
	end
	
	[dummy minx] = min(abs(times-phargs(1)));
	winlen = floor(DEFAULT_CYCLES*srate/freq);
	winloc = minx-linspace(floor(winlen/2), floor(-winlen/2), winlen+1); 
        tmprange = find(winloc>0 & winloc<=frames);
        winloc = winloc(tmprange); % sorting window times
    
	[phaseangles phsamp] = phasedet(data,frames,srate,winloc,freq);
	
    if length(tmprange) ~=  winlen+1
        filtersize = DEFAULT_CYCLES * length(tmprange) / (winlen+1);
        timecenter = median(winloc)/srate*1000+times(1); % center of window in ms
        phaseangles = phaseangles + 2*pi*(timecenter-phargs(1))*freq;
        fprintf('Sorting data epochs by phase at frequency %2.1f Hz: \n', freq);
        fprintf('    Data time limits reached -> now uses a %1.1f cycles (%1.0f ms) window centered at %1.0f ms\n', ...
                filtersize, 1000/freq*filtersize, timecenter);
        fprintf(...
  '    Filter length is %d; Phase has been linearly interpolated to latency at %1.0f ms.\n', ...
                        length(winloc), phargs(1));
    else
        fprintf(...
  'Sorting data epochs by phase at %2.1f Hz in a %1.1f-cycle (%1.0f ms) window centered at %1.0f ms.\n',...  
			freq,DEFAULT_CYCLES,1000/freq*DEFAULT_CYCLES,times(minx));
        fprintf('Phase is computed using a filter of length %d frames.\n',length(winloc));
    end;
	%
	% Reject small (or large) phsamp trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	phargs(2) = phargs(2)/100; % convert rejection rate from % to fraction
	[tmp n] = sort(phsamp); % sort amplitudes
	if phargs(2)>=0
	  n = n(ceil(phargs(2)*length(n))+1:end); % if rej 0, select all trials
	  fprintf('Retaining %d epochs (%g percent) with largest power at the analysis frequency,\n',...
	     length(n),100*(1-phargs(2)));
	else % phargs(2) < 0
	   phargs(2) = 1+phargs(2); % subtract from end
	   n = n(1:floor(phargs(2)*length(n)));
	   fprintf('Retaining %d epochs (%g percent) with smallest power at the analysis frequency,\n',...
                      length(n),phargs(2)*100);
	end
	%
	% Remove low|high-amplitude trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	data = data(:,n); % amp-sort the data, removing rejected-amp trials
	phsamp = phsamp(n);           % amp-sort the amps
	phaseangles = phaseangles(n); % amp-sort the phaseangles
	sortvar = sortvar(n);         % amp-sort the trial indices
	ntrials = length(n);          % number of trials retained
	if ~isempty(auxvar)
	   auxvar = auxvar(:,n);
        end
	%
	% Sort remaining data by phase angle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	topphase = topphase/360*2*pi; % convert from degrees to radians
	phaseangles = -phaseangles;
	ip = find(phaseangles>topphase);
	phaseangles(ip) = phaseangles(ip)-2*pi; % rotate so topphase at top of plot
	
	[phaseangles sortidx] = sort(phaseangles); % sort trials on (rotated) phase
	data    =  data(:,sortidx);                % sort data by phase
	phsamp  =  phsamp(sortidx);                % sort amps by phase
	sortvar = sortvar(sortidx);                % sort input sortvar by phase
	phaseangles = -phaseangles; % Note: phsangles now descend from pi 
	if ~isempty(auxvar)
	   auxvar = auxvar(:,sortidx);
	end
	
	fprintf('Size of data = [%d,%d]\n',size(data,1),size(data,2));
	sortidx = n(sortidx); % return original trial indices in final sorted order
%
% %%%%%%%%%%%%%%% Sort data by amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
elseif exist('ampargs') == 1 % if amplitude-sort
	if length(ampargs) == 4 % find max frequency in specified band
        if exist('psd') == 2
            [pxx,freqs] = psd(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
        else
            [pxx,freqs] = spec(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
        end;
		
		pxx = 10*log10(pxx);
		n = find(freqs >= ampargs(3) & freqs <= ampargs(4));
		if ~length(n)
			freq = ampargs(3);
		end
		[dummy maxx] = max(pxx(n));
		freq = freqs(n(maxx));
	else
		freq = ampargs(3); % else use specified frequency
	end
	
	[dummy minx] = min(abs(times-ampargs(1)));
	winlen = floor(DEFAULT_CYCLES*srate/freq);
	%winloc = minx-[winlen:-1:0]; % ending time version
	winloc = minx-linspace(floor(winlen/2), floor(-winlen/2), winlen+1);
        tmprange = find(winloc>0 & winloc<=frames);
        winloc = winloc(tmprange); % sorting window frames
    
	[phaseangles phsamp] = phasedet(data,frames,srate,winloc,freq);
	
    if length(tmprange) ~=  winlen+1
        filtersize = DEFAULT_CYCLES * length(tmprange) / (winlen+1);
        timecenter = median(winloc)/srate*1000+times(1); % center of window in ms
        phaseangles = phaseangles + 2*pi*(timecenter-ampargs(1))*freq;
        fprintf('Sorting data epochs by phase at frequency %2.1f Hz: \n', freq);
        fprintf(...
'    Data time limits reached -> now uses a %1.1f cycles (%1.0f ms) window centered at %1.0f ms\n', ...
                filtersize, 1000/freq*filtersize, timecenter);
        fprintf(...
'    Filter length is %d; Phase has been linearly interpolated to latency et %1.0f ms.\n', ...
           length(winloc), ampargs(1));
    else
        fprintf(...
'Sorting data epochs by phase at %2.1f Hz in a %1.1f-cycle (%1.0f ms) window centered at %1.0f ms.\n',...  
			freq,DEFAULT_CYCLES,1000/freq*DEFAULT_CYCLES,times(minx));
        fprintf('Phase is computed using a filter of length %d frames.\n',length(winloc));
    end;

	%
	% Reject small (or large) phsamp trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	ampargs(2) = ampargs(2)/100; % convert rejection rate from % to fraction
	[tmp n] = sort(phsamp); % sort amplitudes
	if ampargs(2)>=0
	   n = n(ceil(ampargs(2)*length(n))+1:end); % if rej 0, select all trials
	   fprintf('Retaining %d epochs (%g percent) with largest power at the analysis frequency,\n',...
	      length(n),100*(1-ampargs(2)));
	else % ampargs(2) < 0
		ampargs(2) = 1+ampargs(2); % subtract from end
		n = n(1:floor(ampargs(2)*length(n)));
		fprintf(...
            'Retaining %d epochs (%g percent) with smallest power at the analysis frequency,\n',...
                   length(n),ampargs(2)*100);
	end
	%
	% Remove low|high-amplitude trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	data = data(:,n); % amp-sort the data, removing rejected-amp trials
	phsamp = phsamp(n);           % amp-sort the amps
	phaseangles = phaseangles(n); % amp-sort the phaseangles
	sortvar = sortvar(n);         % amp-sort the trial indices
	ntrials = length(n);          % number of trials retained
	if ~isempty(auxvar)
           auxvar = auxvar(:,n);
        end
	%
	% Sort remaining data by amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%
	[phsamp sortidx] = sort(phsamp); % sort trials on amplitude
	data    =  data(:,sortidx);                % sort data by amp
	phaseangles  =  phaseangles(sortidx);      % sort angles by amp
	sortvar = sortvar(sortidx);                % sort input sortvar by amp
	if ~isempty(auxvar)
	   auxvar = auxvar(:,sortidx);
	end
	fprintf('Size of data = [%d,%d]\n',size(data,1),size(data,2));

	sortidx = n(sortidx); % return original trial indices in final sorted order
%
%%%%%%%%%%%%%%%%%%%%%% Don't Sort trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
elseif Nosort == YES
  fprintf('Not sorting data on input sortvar.\n');
  sortidx = 1:ntrials;	
%
%%%%%%%%%%%%%%%%%%%%%% Sort trials on (mean) value %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
elseif exist('valargs')
    [sttime stframe] = min(abs(times-valargs(1)));
    sttime = times(stframe);
    if length(valargs)>1
        [endtime endframe] = min(abs(times-valargs(2)));
        endtime = times(endframe);
    else
        endframe = stframe;
        endtime = times(endframe);
    end
    if length(valargs)==1 | sttime == endtime
        fprintf('Sorting data on value at time %4.0f ms.\n',sttime);
    elseif length(valargs)>1
        fprintf('Sorting data on mean value between %4.0f and %4.0f ms.\n',...
                sttime,endtime);
    end
    if endframe>stframe
        sortval = mean(data(stframe:endframe,:));
    else
        sortval = data(stframe,:);
    end
    [sortval,sortidx] = sort(sortval);
    if length(valargs)>2
        if valargs(3) <0
            sortidx = sortidx(end:-1:1); % plot largest values on top
                                         % if direction < 0   
        end
    end
    data = data(:,sortidx);
    if ~isempty(auxvar)
        auxvar = auxvar(:,sortidx);
    end
    winloc = [stframe,endframe];
%
%%%%%%%%%%%%%%%%%%%%%% Sort trials on sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
else 
    fprintf('Sorting data on input sortvar.\n');
    [sortvar,sortidx] = sort(sortvar);
    data = data(:,sortidx);
    if ~isempty(auxvar)
        auxvar = auxvar(:,sortidx);
    end
end

%if max(sortvar)<0
%   fprintf('Changing the sign of sortvar: making it positive.\n');
%   sortvar = -sortvar;
%end
%
%%%%%%%%%%%%%%%%%%% Adjust decfactor if phargs or ampargs %%%%%%%%%%%%%%%%%%%%%
%
if decfactor < 0
    decfactor = -decfactor;
    invdec = 1;
else
    invdec = 0;
end;
if decfactor > sqrt(ntrials) % if large, output this many trials
    n = 1:ntrials;
    if exist('phargs') & length(phargs)>1
        if phargs(2)>0
            n = n(ceil(phargs(2)*ntrials)+1:end); % trials after rejection
        elseif phargs(2)<0
            n = n(1:floor(phargs(2)*length(n)));  % trials after rejection
        end
    elseif exist('ampargs') & length(ampargs)>1
        if ampargs(2)>0
            n = n(ceil(ampargs(2)*ntrials)+1:end); % trials after rejection
        elseif ampargs(2)<0
            n = n(1:floor(ampargs(2)*length(n)));  % trials after rejection
        end
    end
    if invdec
        decfactor = (length(n)-avewidth)/decfactor;
    else
        decfactor = length(n)/decfactor;
    end;
end

% 
%%%%%%%%%%%%%%%%%% Smooth data using moving average %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
urdata = data; % save data to compute amp, coher on unsmoothed data
if ~Allampsflag & ~exist('data2') % if imaging potential,
    if avewidth > 1 | decfactor > 1
        if Nosort == YES
            fprintf('Smoothing the data using a window width of %g epochs ',avewidth);
        else
            fprintf('Smoothing the sorted epochs with a %g-epoch moving window.',...
                    avewidth);
        end
        fprintf('\n');
        fprintf('  and a decimation factor of %g\n',decfactor);
        if ~exist('phargs') % if not phase-sorted trials
           [data,outtrials]    = movav(data,1:ntrials,avewidth,decfactor); 
           % Note: movav here sorts using square window
           [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor); 
        else % if phase-sorted trials, use circular / wrap-around smoothing
           backhalf  = floor(avewidth/2);
           fronthalf = floor((avewidth-1)/2);
           if avewidth > 2
            [data,outtrials] = movav([data(:,[(end-backhalf+1):end]),...
                                      data,...
                                      data(:,[1:fronthalf])],...
                                      [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor); 
            [outsort,outtrials] = movav([sortvar((end-backhalf+1):end),...
                                        sortvar,...
                                        sortvar(1:fronthalf)],...
                                        1:(ntrials+backhalf+fronthalf),avewidth,decfactor); 
            % outtrials = 1:ntrials;
           else % avewidth==2
            [data,outtrials] = movav([data(:,end),data],...
                                       [1:(ntrials+1)],avewidth,decfactor); 
            % Note: movav here sorts using square window
            [outsort,outtrials] = movav([sortvar(end) sortvar],...
                                        1:(ntrials+1),avewidth,decfactor); 
            % outtrials = 1:ntrials;
           end
        end
        if ~isempty(auxvar)
          if ~exist('phargs') % if not phase-sorted trials
            [auxvar,tmp] = movav(auxvar,1:ntrials,avewidth,decfactor); 
          else % if phase-sorted trials
           if avewidth>2 
            [auxvar,tmp] = movav([auxvar(:,[(end-backhalf+1):end]),...
                                  auxvar,...
                                  auxvar(:,[1:fronthalf])],...
                                 [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor); 
           else % avewidth==2
            [auxvar,tmp] = movav([auxvar(:,end),auxvar],[1:(ntrials+1)],avewidth,decfactor); 
           end
          end
        end
        fprintf('Output data will be %d frames by %d smoothed trials.\n',...
                frames,length(outtrials));
        fprintf('Outtrials: %3.2f to %4.2f\n',min(outtrials),max(outtrials));
    else % don't smooth
        outtrials = 1:ntrials;
        outsort = sortvar;
    end

    %
    %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~isempty(Caxis) 
        mindat = Caxis(1);
        maxdat = Caxis(2);
        fprintf('Using the specified caxis range of [%g,%g].\n', mindat, maxdat);
    else
        mindat = min(min(data));
        maxdat = max(max(data));
        maxdat =  max(abs([mindat maxdat])); % make symmetrical about 0
        mindat = -maxdat;
        if ~isempty(caxfraction)
            adjmax = (1-caxfraction)/2*(maxdat-mindat);
            mindat = mindat+adjmax;
            maxdat = maxdat-adjmax;
            fprintf(...
                'The caxis range will be %g times the sym. abs. data range -> [%g,%g].\n',...
                caxfraction,mindat,maxdat);
        else
            fprintf(...
                'The caxis range will be the sym. abs. data range -> [%g,%g].\n',...
                mindat,maxdat);
        end
    end
end % if ~Allampsflag & ~exist('data2')

%
%%%%%%%%%%%%%%%%%%%%%%%%%% Set time limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if isnan(timelimits(1))
    timelimits = [min(times) max(times)];
end
fprintf('Data will be plotted between %g and %g ms.\n',timelimits(1),timelimits(2));

%
%%%%%%%%%%%%% Image the aligned/sorted/smoothed data %%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(noshow, 'no')  
    if ~any(isnan(coherfreq))       % if plot three time axes
        image_loy = 3*PLOT_HEIGHT;
    elseif Erpflag == YES   % elseif if plot only one time axes
        image_loy = 1*PLOT_HEIGHT;
    else                    % else plot erp-image only
        image_loy = 0*PLOT_HEIGHT;
    end
    gcapos=get(gca,'Position');
    delete(gca)
    if isempty(topomap)
        image_top = 1;
    else
        image_top = 0.9;
    end
    ax1=axes('Position',...
             [gcapos(1) gcapos(2)+image_loy*gcapos(4) ...