spectopo.m

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

% 02-15-02 scaling by epoch number line 108 - ad, sm & lf
% 03-15-02 add all topoplot options -ad
% 03-18-02 downsampling factor to speed up computation -ad
% 03-27-02 downsampling factor exact calculation -ad
% 04-03-02 added axcopy -sm

% Uses: MATLAB pwelch(), changeunits(), topoplot(), textsc()

function [eegspecdB,freqs,compeegspecdB,resvar,specstd] = spectopo(data,frames,srate,varargin) 
	%headfreqs,chanlocs,limits,titl,freqfac, percent, varargin)

LOPLOTHZ = 1;  % low  Hz to plot
FREQFAC  = 2;  % approximate frequencies/Hz (default)

if nargin<3
   help spectopo
   return
end
if nargin <= 3 | isstr(varargin{1})
	% 'key' 'val' sequence
	fieldlist = { 'freq'          'real'     []                        [] ;
				  'chanlocs'      ''         []                        [] ;
				  'freqrange'     'real'     [0 srate/2]               [] ;
				  'memory'        'string'   {'low' 'high'}           'high' ;
				  'plot'          'string'   {'on' 'off'}             'on' ;
				  'title'         'string'   []                       '';
				  'limits'        'real'     []                       [nan nan nan nan nan nan];
				  'freqfac'       'integer'  []                        FREQFAC;
				  'percent'       'real'     [0 100]                  100 ;
				  'reref'         'string'   { 'averef' 'no' }         'no' ;
				  'boundaries'    'integer'  []                       [] ;
				  'nfft'          'integer'  [1 Inf]                  [] ;
				  'winsize'       'integer'  [1 Inf]                  [] ;
				  'overlap'       'integer'  [1 Inf]                  0 ;
				  'icamode'       'string'   { 'normal' 'sub' }        'normal' ;
				  'weights'       'real'     []                       [] ;
				  'mapnorm'       'real'     []                       [] ;
				  'plotchan'      'integer'  [1:size(data,1)]         [] ;
				  'nicamaps'      'integer'  []                       4 ;
				  'icawinv'       'real'     []                       [] ;
				  'icacomps'      'integer'  []                       [] ;
				  'icamaps'       'integer'  []                       [] };
	
	[g varargin] = finputcheck( varargin, fieldlist, 'spectopo', 'ignore');
	if isstr(g), error(g); end;
	if ~isempty(g.freqrange), g.limits(1:2) = g.freqrange; end;
	if ~isempty(g.weights)
		if isempty(g.freq) | length(g.freq) > 2
            if ~isempty(get(0,'currentfigure')) & strcmp(get(gcf, 'tag'), 'spectopo'), close(gcf); end;
			error('spectopo(): for computing component contribution, one must specify a (single) frequency');
		end;
	end;
else
        if ~isnumeric(data)
           error('spectopo(): Incorrect call format (see >> help spectopo).')
        end
        if ~isnumeric(frames) | round(frames) ~= frames
           error('spectopo(): Incorrect call format (see >> help spectopo).')
        end
        if ~isnumeric(srate)  % 3rd arg must be the sampling rate in Hz
           error('spectopo(): Incorrect call format (see >> help spectopo).')
        end

	if nargin > 3,    g.freq = varargin{1};
	else              g.freq = [];
	end;
	if nargin > 4,	  g.chanlocs = varargin{2};
	else              g.chanlocs = [];
	end;
	if nargin > 5,    g.limits = varargin{3};
	else              g.limits = [nan nan nan nan nan nan];
	end;
	if nargin > 6,    g.title = varargin{4};
	else              g.title = '';
	end;
	if nargin > 7,    g.freqfac = varargin{5};
	else              g.freqfac = FREQFAC;
	end;
	if nargin > 8,    g.percent = varargin{6};
	else              g.percent = 100;
	end;
	if nargin > 10,    g.reref = 'averef';
	else               g.reref = 'no';
	end;
	g.weights = [];
	g.icamaps = [];
end;
if g.percent > 1
	g.percent = g.percent/100; % make it from 0 to 1
end;
if ~isempty(g.freq) & isempty(g.chanlocs)
	error('spectopo: need channel location file');
end;

data = reshape(data, size(data,1), size(data,2)*size(data,3));
if frames == 0
  frames = size(data,2); % assume one epoch
end
%if ~isempty(g.plotchan) & g.plotchan == 0 & strcmpi(g.icamode, 'sub')
%    if ~isempty(get(0,'currentfigure')) & strcmp(get(gcf, 'tag'), 'spectopo'), close(gcf); end;
%    error('Cannot plot data component at all channels (option not implemented)');
%end;
if ~isempty(g.freq) & min(g.freq)<0
    if ~isempty(get(0,'currentfigure')) & strcmp(get(gcf, 'tag'), 'spectopo'), close(gcf); end;
   fprintf('spectopo(): freqs must be >=0 Hz\n');
   return
end
epochs = round(size(data,2)/frames);
if frames*epochs ~= size(data,2)
   error('Spectopo: non-integer number of epochs');
end
if ~isempty(g.weights)
    ncompsori = size(g.weights,1);
    if isempty(g.icawinv)
        g.icawinv = pinv(g.weights); % maps
    end;
	if ~isempty(g.icacomps)
		g.weights = g.weights(g.icacomps, :);
		g.icawinv = g.icawinv(:,g.icacomps);
	else 
		g.icacomps = [1:size(g.weights,1)];
	end;
end;
compeegspecdB = [];
resvar = NaN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute channel spectra using pwelch()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epoch_subset = 1:epochs;
if g.percent ~= 1 & epochs > 1
    epoch_subset = zeros(1,epochs);
    nb = ceil( g.percent*epochs);
    while nb>0
        index = ceil(rand*epochs);
        if ~epoch_subset(index)
            epoch_subset(index) = 1;
            nb = nb-1;
        end;
    end;        
    epoch_subset = find(epoch_subset == 1);
    fprintf('Randomly selecting %d of %d data epochs for analysis...\n', length(epoch_subset),epochs);
end;
if isempty(g.weights)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% compute data spectrum
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	fprintf('Computing spectra')
	[eegspecdB freqs specstd] = spectcomp( data, frames, srate, epoch_subset, g);
    if ~isempty(g.mapnorm) % normalize by component map RMS power (if data contain 1 component
        disp('Scaling spectrum by component RMS of scalp map power');
        eegspecdB       = sqrt(mean(g.mapnorm.^4)) * eegspecdB;
        % the idea is to take the RMS of the component activity (compact) projected at each channel
        % spec = sqrt( power(g.mapnorm(1)*compact).^2 + power(g.mapnorm(2)*compact).^2 + ...)
        % spec = sqrt( g.mapnorm(1)^4*power(compact).^2 + g.mapnorm(1)^4*power(compact).^2 + ...)
        % spec = sqrt( g.mapnorm(1)^4 + g.mapnorm(1)^4 + ... )*power(compact)
    end;
	eegspecdB = 10*log10(eegspecdB);
    specstd   = 10*log10(specstd);
    fprintf('\n');
else
	% compute data spectrum
	if isempty(g.plotchan) | g.plotchan == 0
		fprintf('Computing spectra')
		[eegspecdB freqs specstd] = spectcomp( data, frames, srate, epoch_subset, g);	
        fprintf('\n'); % log below
	else
		fprintf('Computing spectra at specified channel')
		g.reref = 'no';
		[eegspecdB freqs specstd] = spectcomp( data(g.plotchan,:), frames, srate, epoch_subset, g);
        fprintf('\n'); % log below
	end;
	g.reref = 'no';
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% select channel and spectrum
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if isempty(g.plotchan) % find channel of minimum power
		[tmp indexfreq] = min(abs(g.freq-freqs));
		[tmp g.plotchan] = min(eegspecdB(:,indexfreq));
		fprintf('Channel %d has maximum power at %g\n', g.plotchan, g.freq);
		eegspecdBtoplot = eegspecdB(g.plotchan, :);		
	elseif g.plotchan == 0
		fprintf('Computing RMS power at all channels\n');
		eegspecdBtoplot = sqrt(mean(eegspecdB.^2,1)); % RMS before log as for components
	else 
		eegspecdBtoplot = eegspecdB;
	end;
    specstd   = 10*log10(specstd);
    eegspecdB = 10*log10(eegspecdB);
	eegspecdBtoplot = 10*log10(eegspecdBtoplot);
    
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% compute component spectra
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
    newweights = g.weights;
	if strcmp(g.memory, 'high') & strcmp(g.icamode, 'normal')
        fprintf('Computing component spectra: ')
		[compeegspecdB freqs] = spectcomp( newweights*data, frames, srate, epoch_subset, g);
	else % in case out of memory error, multiply conmponent sequencially
        if strcmp(g.icamode, 'sub') & ~isempty(g.plotchan) & g.plotchan == 0
            % scan all electrodes
            fprintf('Computing component spectra at each channel: ')
            for index = 1:size(data,1)
                g.plotchan = index;
                [compeegspecdB(:,:,index) freqs] = spectcomp( data, frames, srate, epoch_subset, g, newweights);
            end;
            g.plotchan = 0;
        else
            fprintf('Computing component spectra: ')
            [compeegspecdB freqs] = spectcomp( data, frames, srate, epoch_subset, g, newweights);
        end;
	end;
	fprintf('\n');
    
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% rescale spectra with respect to projection (rms = root mean square)
    % all channel: component_i power = rms(inverseweigths(component_i)^2)*power(activation_component_i);
    % one channel: component_i power = inverseweigths(channel_j,component_i)^2)*power(activation_component_i);
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if strcmpi(g.icamode, 'normal')
        for index = 1:size(compeegspecdB,1) 
            if g.plotchan == 0 % normalize by component scalp map power
                compeegspecdB(index,:) = 10*log10( sqrt(mean(g.icawinv(:,index).^4)) * compeegspecdB(index,:) );
            else 
                compeegspecdB(index,:) = 10*log10( g.icawinv(g.plotchan,index)^2 * compeegspecdB(index,:) );
            end;
        end;
    else % already spectrum of data-components
        compeegspecdB = 10*log10( compeegspecdB );
    end;
        
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% select components to plot
	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if isempty(g.icamaps)
		[tmp indexfreq] = min(abs(g.freq-freqs));
		g.icafreqsval   = compeegspecdB(:, indexfreq);
		[g.icafreqsval g.icamaps] = sort(g.icafreqsval);
		if strcmp(g.icamode, 'normal')
			g.icamaps = g.icamaps(end:-1:1);
			g.icafreqsval = g.icafreqsval(end:-1:1);
		end;
		if g.nicamaps < length(g.icamaps), g.icamaps = g.icamaps(1:g.nicamaps); end;
	else 
		[tmp indexfreq] = min(abs(g.freq-freqs));
		g.icafreqsval   = compeegspecdB(g.icamaps, indexfreq);
	end;
end;
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute axis and caxis g.limits
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(g.limits)<1 | isnan(g.limits(1))
   g.limits(1) = LOPLOTHZ;
end

if ~isempty(g.freq)
	if length(g.limits)<2 | isnan(g.limits(2))
		maxheadfreq = max(g.freq);
		if rem(maxheadfreq,5) ~= 0
			g.limits(2) = 5*ceil(maxheadfreq/5);
		else
			g.limits(2) = maxheadfreq*1.1;
		end
	end
	
	g.freq = sort(g.freq);          % Determine topoplot frequencies
	freqidx = zeros(1,length(g.freq)); % Do not interpolate between freqs
	for f=1:length(g.freq)
		[tmp fi] = min(abs(freqs-g.freq(f)));
		freqidx(f)=fi;
	end
else
    if isnan(g.limits(2))
        g.limits(2) = srate/2;
    end;
end;

[tmp maxfreqidx] = min(abs(g.limits(2)-freqs)); % adjust max frequency
[tmp minfreqidx] = min(abs(g.limits(1)-freqs)); % adjust min frequency

if length(g.limits)<3|isnan(g.limits(3))
	reallimits(1) = min(min(eegspecdB(:,minfreqidx:maxfreqidx)));
else 
	reallimits(1) = g.limits(3);
end
if length(g.limits)<4|isnan(g.limits(4))
	reallimits(2) = max(max(eegspecdB(:,minfreqidx:maxfreqidx)));
	dBrange = reallimits(2)-reallimits(1);   % expand range a bit beyond data g.limits
	reallimits(1) = reallimits(1)-dBrange/7;
	reallimits(2) = reallimits(2)+dBrange/7;
else 
	reallimits(2) = g.limits(4);
end

if length(g.limits)<5 % default caxis plotting g.limits
  g.limits(5) = nan;
end
if length(g.limits)<6 
  g.limits(6) = nan;
end

if isnan(g.limits(5))+isnan(g.limits(6)) == 1
   fprintf('spectopo(): limits 5 and 6 must either be given or nan\n');
   return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot spectrum of each channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmpi(g.plot, 'on')
    mainfig = gca; axis off;
    if ~isempty(g.freq)
        specaxes = sbplot(3,4,[5 12], 'ax', mainfig); 
    end;
    
    if isempty(g.weights)
        pl=plot(freqs(1:maxfreqidx),eegspecdB(:,1:maxfreqidx)');
    else 
        pl=plot(freqs(1:maxfreqidx),eegspecdBtoplot(:,1:maxfreqidx)');
    end;
    set(pl,'LineWidth',2);
    set(gca,'TickLength',[0.02 0.02]);
    try, 
        axis([freqs(minfreqidx) freqs(maxfreqidx) reallimits(1) reallimits(2)]);
    catch, disp('Could not adjust axis'); end;
    xl=xlabel('Frequency (Hz)');
    set(xl,'fontsize',16);
    % yl=ylabel('Rel. Power (dB)');
    yl=ylabel('Power 10*log_{10}(\muV^{2}/Hz)');
    set(yl,'fontsize',16);
    set(gca,'fontsize',16)
    box off;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   plot component contribution   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
colrs = {'r','b','g','m','c'}; % component spectra trace colors
if ~isempty(g.weights)
    if strcmpi(g.plot, 'on')
        if strcmpi(g.icamode, 'sub')
            set(pl,'LineWidth',5, 'color', 'k');
        else 
            set(pl, 'linewidth', 2, 'color', 'k');
        end;
        
        hold on;
        for f=1:length(g.icamaps)
            colr = colrs{mod((f-1),5)+1};
            pl2=plot(freqs(1:maxfreqidx),compeegspecdB(g.icamaps(f),1:maxfreqidx)',colr);
        end
        othercomps = setdiff(1:size(compeegspecdB,1), g.icamaps);
        if ~isempty(othercomps)
            pl2=plot(freqs(1:maxfreqidx),compeegspecdB(othercomps,1:maxfreqidx)');
        end;
        if length(g.limits)<3|isnan(g.limits(3))
            newaxis = axis;
            newaxis(3) = min(newaxis(3), min(min(compeegspecdB(:,1:maxfreqidx))));
            newaxis(4) = max(newaxis(4), max(max(compeegspecdB(:,1:maxfreqidx))));
            axis(newaxis);
        end;
	end;
    
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% indicate component contribution %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    maxdatadb = max(eegspecdBtoplot(:,freqidx(1)));
    [tmp indexfreq] = min(abs(g.freq-freqs));
    for index = 1:length(g.icacomps)
        if strcmp(g.icamode, 'normal')
            % note: maxdatadb = eegspecdBtoplot (RMS power of data)
            resvar(index)  = 100*exp(-(maxdatadb-compeegspecdB(index, indexfreq))/10*log(10));
            fprintf('Component %d percent relative variance:%6.2f\n', g.icacomps(index), resvar(index));
        else
            if g.plotchan == 0
                resvartmp = [];
                for chan = 1:size(eegspecdB,1) % scan channels
                    resvartmp(chan)  = 100 - 100*exp(-(eegspecdB(chan,freqidx(1))-compeegspecdB(index, indexfreq, chan))/10*log(10));
                end;
                resvar(index) = mean(resvartmp); % mean contribution for all channels
                stdvar(index) = std(resvartmp);
                fprintf('Component %d percent variance accounted for:%6.2f