% SPECTOPO - Plot the power spectral density (PSD) of winsize length segments of data
% epochs at all channels as a bundle of traces. At specified frequencies,
% plot the relative topographic distribution of PSD. If available, uses
% PWELCH from the Matlab signal processing toolbox, else the EEGLAB SPEC
% function. Plots the mean spectrum for all of the supplied data, not just
% the pre-stimulus baseline.
% Usage:
% >> spectopo(data, frames, srate);
% >> [spectra,freqs,speccomp,contrib,specstd] = ...
% spectopo(data, frames, srate, 'key1','val1', 'key2','val2' ...);
% Inputs:
% data = If 2-D (nchans,time_points); % may be a continuous single epoch,
% else a set of concatenated data epochs, else a 3-D set of data
% epochs (nchans,frames,epochs)
% frames = frames per epoch {default|0 -> data length}
% srate = sampling rate per channel (Hz)
%
% Optional 'keyword',[argument] input pairs:
% 'freq' = [float vector (Hz)] vector of frequencies at which to plot power
% scalp maps, or else a single frequency at which to plot component
% contributions at a single channel (see also 'plotchan')
% 'chanlocs' = [electrode locations filename or EEG.chanlocs structure]
% For format, see >> topoplot example
% 'limits' = [xmin xmax ymin ymax cmin cmax] axis limits. Sets x, y, and color
% axis limits. May omit final values or use NaNs.
% Ex: [0 60 NaN NaN -10 10], [0 60], ...
% Default color limits are symmetric around 0 and are different
% for each scalp map {default|all NaN's: from the data limits}
% 'title' = [quoted string] plot title {default: none}
% 'freqfac' = [integer] ntimes to oversample (to adjust frequency resolution) {default: 1}
% 'nfft' = [integer] Data points to zero-pad data windows to (overwrites 'freqfac')
% 'winsize' = [integer] window size in data points {default: Sampling Rate}
% 'overlap' = [integer] window overlap in data points {default: 0}
% 'percent' = [float 0 to 100] percent of the data to sample for computing the
% spectra. Values < 100 speed up the computation. {default: 100}
% 'freqrange' = [min max] frequency range to plot. Changes x-axis limits {default:
% 1 Hz for the min and Nyquist (srate/2) for the max. If specified
% power distribution maps are plotted, the highest mapped frequency
% determines the max freq}.
% 'wintype' = ['hamming','blackmanharris'] Window type used on the power spectral
% density estimation. The Blackman-Harris windows offers better attenuation
% than Hamming windows, but lower spectral resolution. {default: 'hamming'}
% 'blckhn' = [integer] Parameter to scale the windows length when Blackman-Harris
% windows is used in 'wintype' {default: 2}
% 'reref' = ['averef'|'off'] convert data to average reference {default: 'off'}
% 'mapnorm' = [float vector] If 'data' contain the activity of an independent
% component, this parameter should contain its scalp map. In this case
% the spectrum amplitude will be scaled to component RMS scalp power.
% Useful for comparing component strengths {default: none}
% 'boundaries' = data point indices of discontinuities in the signal {default: none}
% 'plot' = ['on'|'off'] 'off' -> disable plotting {default: 'on'}
% 'rmdc' = ['on'|'off'] 'on' -> remove DC {default: 'off'}
% 'plotmean' = ['on'|'off'] 'on' -> plot the mean channel spectrum {default: 'off'}
% 'plotchans' = [integer array] plot only specific channels {default: all}
% 'verbose' = ['on'|'off'] 'on' shows information on command line {default: 'on'}
%
% Optionally plot component contributions:
% 'weights' = ICA unmixing matrix. Here, 'freq' (above) must be a single frequency.
% ICA maps of the N ('nicamaps') components that account for the most
% power at the selected frequency ('freq') are plotted along with
% the spectra of the selected channel ('plotchan') and components
% ('icacomps').
% 'plotchan' = [integer] channel at which to compute independent conmponent
% contributions at the selected frequency ('freq'). If 0, plot RMS
% power at all channels. {defatul|[] -> channel with highest power
% at specified 'freq' (above)). Do not confuse with
% 'plotchans' which select channels for plotting.
% 'mapchans' = [int vector] channels to plot in topoplots {default: all}
% 'mapframes'= [int vector] frames to plot {default: all}
% 'nicamaps' = [integer] number of ICA component maps to plot {default: 4}.
% 'icacomps' = [integer array] indices of ICA component spectra to plot ([] -> all).
% 'icamode' = ['normal'|'sub'] in 'sub' mode, instead of computing the spectra of
% individual ICA components, the function computes the spectrum of
% the data minus their contributions {default: 'normal'}
% 'icamaps' = [integer array] force plotting of selected ICA component maps
% {default: [] = the 'nicamaps' largest contributing components}.
% 'icawinv' = [float array] inverse component weight or mixing matrix. Normally,
% this is computed by inverting the ICA unmixing matrix 'weights' (above).
% However, if any components were removed from the supplied 'weights'mapchans
% then the component maps will not be correctly drawn and the 'icawinv'
% matrix should be supplied here {default: from component 'weights'}
% 'memory' = ['low'|'high'] a 'low' setting will use less memory for computing
% component activities, will take longer {default: 'high'}
%
% Replotting options:
% 'specdata' = [freq x chan array ] spectral data
% 'freqdata' = [freq] array of frequencies
%
% Topoplot options:
% other 'key','val' options are propagated to TOPOPLOT for map display
% (See >> help topoplot)
%
% Outputs:
% spectra = (nchans,nfreqs) power spectra (mean power over epochs), in dB
% freqs = frequencies of spectra (Hz)
% speccomp = component spectra (ncomps,nfreqs). Warning, only the function
% only computes the spectrum of the components given as input using
% the 'icacomps' parameter. Other component spectrum are filled
% with 0.
% contrib = contribution of each component (if 'icamode' is 'normal', relative
% variance, if 'icamode' is 'sub', percent variance accounted for)
% specstd = spectrum standard deviation in dB
%
% Notes: The original input format is still functional for backward compatibility.
% PSD has been replaced by PWELCH (see Matlab note 24750 on their web site)
%
% Non-backward compatible change (Nov 15 2015):
% Default winsize was set to the sampling rate (giving a default window
% length of 1 sec). Also, the y-axis label in the plot was corrected
% to read, "Log Power Spectral Density 10*log_{10}(\muV^{2}/Hz)'
% Finally, when winsize is not a power of 2, it is no longer promoted to
% the next higher power of 2. Thanks to Andreas Widmann for his comments.
%
%
% Authors: Scott Makeig, Arnaud Delorme & Marissa Westerfield,
% SCCN/INC/UCSD, La Jolla, 3/01
%
% See also: TIMTOPO, ENVTOPO, TFTOPO, TOPOPLOT
% Copyright (C) 3/01 Scott Makeig & Arnaud Delorme & Marissa Westerfield, SCCN/INC/UCSD,
% scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
% 3-20-01 added limits arg -sm
% 01-25-02 reformated help & license -ad
% 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)
% formerly: ... headfreqs,chanlocs,limits,titl,freqfac, percent, varargin)
icadefs;
LOPLOTHZ = 1; % low Hz to plot
FREQFAC = 1; % approximate frequencies/Hz (default)
allcolors = { [0 0.7500 0.7500]
[1 0 0]
[0 0.5000 0]
[0 0 1]
[0.2500 0.2500 0.2500]
[0.7500 0.7500 0]
[0.7500 0 0.7500] }; % colors from real plots };
if nargin<3
help spectopo
return
end
if nargin <= 3 || ischar(varargin{1})
% 'key' 'val' sequence
fieldlist = { 'freq' 'real' [] [] ;
'specdata' 'real' [] [] ;
'freqdata' 'real' [] [] ;
'chanlocs' '' [] [] ;
'freqrange' 'real' [0 srate/2] [] ;
'memory' 'string' {'low','high'} 'high' ;
'plot' 'string' {'on','off'} 'on' ;
'plotmean' 'string' {'on','off'} 'off' ;
'title' 'string' [] '';
'limits' 'real' [] [nan nan nan nan nan nan];
'freqfac' 'integer' [] FREQFAC;
'percent' 'real' [0 100] 100 ;
'reref' 'string' { 'averef','off','no' } 'off' ;
'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)] [] ;
'plotchans' 'integer' [1:size(data,1)] [] ;
'nicamaps' 'integer' [] 4 ;
'blckhn' 'integer' [] 2 ;
'icawinv' 'real' [] [] ;
'icacomps' 'integer' [] [] ;
'icachansind' 'integer' [] [1:size(data,1)] ;
'icamaps' 'integer' [] [] ;
'rmdc' 'string' {'on','off'} 'off';
'verbose' 'string' {'on','off'} 'on';
'wintype' 'string' {} 'hamming';
'mapchans' 'integer' [1:size(data,1)] []
'mapframes' 'integer' [1:size(data,2)] []};
[g, varargin] = finputcheck( varargin, fieldlist, 'spectopo', 'ignore');
if ischar(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 = 'off';
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(): needs channel location information');
end
if isempty(g.weights) && ~isempty(g.plotchans)
data = data(g.plotchans,:);
if ~isempty(g.chanlocs)
g.chanlocs = g.chanlocs(g.plotchans);
end
end
if strcmpi(g.rmdc, 'on')
data = data - repmat(mean(data,2), [ 1 size(data,2) 1]);
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
g.chanlocs2 = g.chanlocs;
if ~isempty(g.specdata)
eegspecdB = g.specdata;
freqs = g.freqdata;
else
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 = ones(1,epochs);
if g.percent ~= 1 && epochs == 1
fprintf('Selecting the first %2.1f%% of data for analysis...\n', g.percent*100);
frames = round(size(data,2)*g.percent);
data = data(:, 1:frames);
g.boundaries(find(g.boundaries > frames)) = [];
if ~isempty(g.boundaries)
g.boundaries(end+1) = frames;
end;
end
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);
else
epoch_subset = find(epoch_subset == 1);
end
if isempty(g.weights)
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute data spectra
%%%%%%%%%%%%%%%%%%%%%%%%%%%
myfprintf(g.verbose, '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
myfprintf(g.verbose, 'Scaling spectrum by component RMS of scalp map power\n');
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
tmpc = find(eegspecdB(:,1)); % > 0 power chans
tmpindices = find(eegspecdB(:,1) == 0);
if ~isempty(tmpindices)
zchans = int2str(tmpindices); % 0-power chans
else zchans = [];
end
if length(tmpc) ~= size(eegspecdB,1)
myfprintf(g.verbose, '\nWarning: channels [%s] have 0 values, so will be omitted from the display', ...
zchans);
eegspecdB = eegspecdB(tmpc,:);
if ~isempty(specstd), specstd = specstd(tmpc,:); end
if ~isempty(g.chanlocs)
g.chanlocs2 = g.chanlocs(tmpc);
end
end
eegspecdB = 10*log10(eegspecdB);
specstd = 10*log10(specstd);
myfprintf(g.verbose, '\n');
else
% compute data spectrum
if isempty(g.plotchan) || g.plotchan == 0
myfprintf(g.verbose, 'Computing spectra')
[eegspecdB freqs specstd] = spectcomp( data, frames, srate, epoch_subset, g);
myfprintf(g.verbose, '\n'); % log below
else
myfprintf(g.verbose, 'Computing spectra at specified channel')
g.reref = 'off';
[eegspecdB freqs specstd] = spectcomp( data(g.plotchan,:), frames, srate, epoch_subset, g);
myfprintf(g.verbose, '\n'); % log below
end
g.reref = 'off';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select channels and spectra
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(g.plotchan) % find channel of minimum power
[tmp indexfreq] = min(abs(g.freq-freqs));
[tmp g.plotchan] = min(eegspecdB(:,indexfreq));
myfprintf(g.verbose, 'Channel %d has maximum power at %g\n', g.plotchan, g.freq);
eegspecdBtoplot = eegspecdB(g.plotchan, :);
elseif g.plotchan == 0
myfprintf(g.verbose, 'Computing RMS power at all channels\n');
eegspecdBtoplot = sqrt(mean(eegspecdB.^2,1)); % RMS before log as for components
else
eegspecdBtoplot = eegspecdB;
end
tmpc = find(eegspecdB(:,1)); % > 0 power chans
zchans = int2str(find(eegspecdB(:,1) == 0)); % 0-power chans
if length(tmpc) ~= size(eegspecdB,1)
myfprintf(g.verbose, '\nWarning: channels [%s] have 0 values, so will be omitted from the display', ...
zchans);
eegspecdB = eegspecdB(tmpc,:);
if ~isempty(specstd), specstd = specstd(tmpc,:); end
if ~isempty(g.chanlocs)
g.chanlocs2 = g.chanlocs(tmpc);
end
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')
myfprintf(g.verbose, 'Computing component spectra: ')
[compeegspecdB freqs] = spectcomp( newweights*data(:,:), frames, srate, epoch_subset, g);
else % in case out of memory error, multiply conmponent sequentially
if strcmp(g.icamode, 'sub') && ~isempty(g.plotchan) && g.plotchan == 0
% scan all electrodes
myfprintf(g.verbose, '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
myfprintf(g.verbose, 'Computing component spectra: ')
[compeegspecdB freqs] = spectcomp( data, frames, srate, epoch_subset, g, newweights);
end
end
myfprintf(g.verbose, '\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rescale spectra with respect to projection (rms = root mean square)
% all channels: 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
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 % no freq specified
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)'); % old command
if strcmpi(g.plotmean, 'on'), specdata = mean(eegspecdB,1); % average channels
else specdata = eegspecdB;
end
for index = 1:size(specdata,1) % scan channels
tmpcol = allcolors{mod(index, length(allcolors))+1};
command = [ 'disp(''Channel ' int2str(index) ''')' ];
pl(index)=plot(freqs(1:maxfreqidx),specdata(index,1:maxfreqidx)', ...
'color', tmpcol, 'ButtonDownFcn', command); hold on;
end
else
for index = 1:size(eegspecdBtoplot,1)
tmpcol = allcolors{mod(index, length(allcolors))+1};
command = [ 'disp(''Channel ' int2str(g.plotchan(index)) ''')' ];
pl(index)=plot(freqs(1:maxfreqidx),eegspecdBtoplot(index,1:maxfreqidx)', ...
'color', tmpcol, 'ButtonDownFcn', command); hold on;
end
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',AXES_FONTSIZE_L);
% yl=ylabel('Rel. Power (dB)');
yl=ylabel('Log Power Spectral Density 10*log_{10}(\muV^{2}/Hz)');%yl=ylabel('Power 10*log_{10}(\muV^{2}/Hz)');
set(yl,'fontsize',AXES_FONTSIZE_L);
set(gca,'fontsize',AXES_FONTSIZE_L)
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};
command = [ 'disp(''Component ' int2str(g.icamaps(f)) ''')' ];
pl2(index)=plot(freqs(1:maxfreqidx),compeegspecdB(g.icamaps(f),1:maxfreqidx)', ...
'color', colr, 'ButtonDownFcn', command); hold on;
end
othercomps = setdiff_bc(1:size(compeegspecdB,1), g.icamaps);
if ~isempty(othercomps)
for index = 1:length(othercomps)
tmpcol = allcolors{mod(index, length(allcolors))+1};
command = [ 'disp(''Component ' int2str(othercomps(index)) ''')' ];
pl(index)=plot(freqs(1:maxfreqidx),compeegspecdB(othercomps(index),1:maxfreqidx)', ...
'color', tmpcol, 'ButtonDownFcn', command); hold on;
end
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));
myfprintf(g.verbose, '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);
myfprintf(g.verbose, 'Component %d percent variance accounted for: %6.2f ± %3.2f\n', ...
g.icacomps(index), resvar(index), stdvar(index));
else
resvar(index) = 100 - 100*exp(-(maxdatadb-compeegspecdB(index, indexfreq))/10*log(10));
myfprintf(g.verbose, 'Component %d percent variance accounted for: %6.2f\n', g.icacomps(index), resvar(index));
end
end
end
% for icamode=sub and plotchan == 0 -> take the RMS power of all channels
% -----------------------------------------------------------------------
if ndims(compeegspecdB) == 3
compeegspecdB = exp( compeegspecdB/10*log(10) );
compeegspecdB = sqrt(mean(compeegspecdB.^2,3)); % RMS before log (dim1=comps, dim2=freqs, dim3=chans)
compeegspecdB = 10*log10( compeegspecdB );
end
end
if ~isempty(g.freq) && strcmpi(g.plot, 'on')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot vertical lines through channel trace bundle at each headfreq
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(g.weights)
for f=1:length(g.freq)
hold on;
plot([freqs(freqidx(f)) freqs(freqidx(f))], ...
[min(eegspecdB(:,freqidx(f))) max(eegspecdB(:,freqidx(f)))],...
'k','LineWidth',2.5);
end
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot vertical line at comp analysis freq
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mincompdB = min([min(eegspecdB(:,freqidx(1))) min(compeegspecdB(:,freqidx(1)))]);
maxcompdB = max([max(eegspecdB(:,freqidx(1))) max(compeegspecdB(:,freqidx(1)))]);
hold on;
plot([freqs(freqidx(1)) freqs(freqidx(1))],[mincompdB maxcompdB],'k', 'LineWidth',2.5);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% create axis for topoplot
%%%%%%%%%%%%%%%%%%%%%%%%%%
tmpmainpos = get(gca, 'position');
headax = zeros(1,length(g.freq));
for f=1:length(g.freq)+length(g.icamaps)
headax(f) = sbplot(3,length(g.freq)+length(g.icamaps),f, 'ax', mainfig);
axis([-1 1 -1 1]);
%axis x coords and use
tmppos = get(headax(f), 'position');
allaxcoords(f) = tmppos(1);
allaxuse(f) = 0;
end
large = sbplot(1,1,1, 'ax', mainfig);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute relative positions on plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(g.weights)
freqnormpos = tmpmainpos(1) + tmpmainpos(3)*(freqs(freqidx(1))-g.limits(1))/(g.limits(2)-g.limits(1));
for index = 1:length(g.icamaps)+1
[realpos(index) allaxuse] = closestplot( freqnormpos, allaxcoords, allaxuse );
end
% put the channel plot a little bit higher
tmppos = get(headax(realpos(1)), 'position');
tmppos(2) = tmppos(2)+0.04;
set(headax(realpos(1)), 'position', tmppos);
else
realpos = 1:length(g.freq); % indices giving order of plotting positions
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot connecting lines using changeunits()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for f=1:length(g.freq)+length(g.icamaps)
if ~isempty(g.weights)
if f>length(g.freq) % special case of ICA components
from = changeunits([freqs(freqidx(1)),g.icafreqsval(f-1)],specaxes,large);
%g.icafreqsval contains the sorted frequency values at the specified frequency
else
from = changeunits([freqs(freqidx(f)),maxcompdB],specaxes,large);
end
else
from = changeunits([freqs(freqidx(f)),max(eegspecdB(:,freqidx(f)))],specaxes,large);
end
pos = get(headax(realpos(f)),'position');
to = changeunits([0,0],headax(realpos(f)),large)+[0 -min(pos(3:4))/2.5];
hold on;
if f<=length(g.freq)
colr = 'k';
else
colr = colrs{mod((f-2),5)+1};
end
li(realpos(f)) = plot([from(1) to(1)],[from(2) to(2)],colr,'LineWidth',PLOT_LINEWIDTH_S);
axis([0 1 0 1]);
axis off;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot selected channel head using topoplot()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
myfprintf(g.verbose, 'Plotting scalp distributions: ')
for f=1:length(g.freq)
axes(headax(realpos(f)));
topodata = eegspecdB(:,freqidx(f))-nan_mean(eegspecdB(:,freqidx(f)));
if isnan(g.limits(5)),
maplimits = 'absmax';
else
maplimits = [g.limits(5) g.limits(6)];
end
%
% If 1 channel in g.plotchan
%
if ~isempty(g.plotchan) && g.plotchan ~= 0
% if ~isempty(varargin) % if there are extra topoplot() flags
% topoplot(g.plotchan,g.chanlocs,'electrodes','off', ...
% 'style', 'blank', 'emarkersize1chan', 10, varargin{:});
% else
topoplot(g.plotchan,g.chanlocs,'electrodes','off', ...
'style', 'blank', 'emarkersize1chan', 10);
% end
if isstruct(g.chanlocs)
tl=title(g.chanlocs(g.plotchan).labels);
else
tl=title([ 'c' int2str(g.plotchan)]);
end
else % plot all channels in g.plotchans
if isempty(g.mapframes) || g.mapframes(1) == 0
g.mapframes = 1:size(eegspecdB,1); % default to plotting all chans
end
if ~isempty(varargin)
topoplot(topodata(g.mapframes),g.chanlocs2,'maplimits',maplimits, varargin{:});
else
topoplot(topodata(g.mapframes),g.chanlocs2,'maplimits',maplimits);
end
if f<length(g.freq)
tl=title([num2str(freqs(freqidx(f)), '%3.1f')]);
else
tl=title([num2str(freqs(freqidx(f)), '%3.1f') ' Hz']);
end
end
set(tl,'fontsize',AXES_FONTSIZE_L);
axis square;
drawnow
myfprintf(g.verbose, '.');
end
myfprintf(g.verbose, '\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot independent components
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(g.weights)
% use headaxe from 2 to end (reserved earlier)
if realpos(1) == max(realpos), plotcolbar(g); end
% make the line with the scalp topoplot thicker than others
set(li(realpos(1)), 'linewidth', 2.5);
if isempty(g.mapchans) || g.mapchans == 0
g.mapchans = 1:size(g.icawinv,1); % default to plotting all chans
end
for index = 1:length(g.icamaps)
axes(headax(realpos(index+1)));
compnum = g.icamaps(index);
topoplot(g.icawinv(g.mapchans,compnum).^2,g.chanlocs,varargin{:});
tl=title(int2str(g.icacomps(compnum)));
set(tl,'fontsize',16);
axis square;
drawnow
try,
if strcmpi(g.icamode, 'normal')
set(gca, 'userdata', ['text(-0.6, -0.6, ''Rel. Var.: ' sprintf('%6.2f', resvar(g.icacomps(compnum))) ''');'] );
else
set(gca, 'userdata', ['text(-0.6, -0.6, ''PVAF: ' sprintf('%6.2f', resvar(g.icacomps(compnum))) ''');'] );
end
catch, end
if realpos(index+1) == max(realpos), plotcolbar(g); end
end
else
plotcolbar(g);
end
end
%%%%%%%%%%%%%%%%
% Draw title
%%%%%%%%%%%%%%%%
if ~isempty(g.title) && strcmpi(g.plot, 'on')
tl = textsc(g.title);
set(tl,'fontsize',15)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% return component spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(g.weights) && nargout >= 3
tmp = compeegspecdB;
compeegspecdB = zeros(ncompsori, size(tmp,2));
compeegspecdB(g.icacomps,:) = tmp;
end
%%%%%%%%%%%%%%%%
% Turn on axcopy (disabled to allow to click on curves)
%%%%%%%%%%%%%%%%
if strcmpi(g.plot, 'on')
disp('Click on each trace for channel/component index');
axcopy(gcf, 'if ~isempty(get(gca, ''''userdata'''')), eval(get(gca, ''''userdata'''')); end;');
% will not erase the commands for the curves
end
%%%%%%%%%%%%%%%%
% Plot color bar
%%%%%%%%%%%%%%%%
function plotcolbar(g)
cb=cbar;
pos = get(cb,'position');
set(cb,'position',[pos(1) pos(2) 0.03 pos(4)]);
set(cb,'fontsize',12);
try
if isnan(g.limits(5))
ticks = get(cb,'ytick');
[tmp zi] = find(ticks == 0);
ticks = [ticks(1) ticks(zi) ticks(end)];
set(cb,'ytick',ticks);
set(cb,'yticklabel',strvcat('-',' ','+'));
end
catch, end; % in a single channel is plotted
%%%%%%%%%%%%%%%%%%%%%%%
% function closest plot
%%%%%%%%%%%%%%%%%%%%%%%
function [index, usedplots] = closestplot(xpos, xcentercoords, usedplots);
notused = find(usedplots == 0);
xcentercoords = xcentercoords(notused);
[tmp index] = min(abs(xcentercoords-xpos));
index = notused(index);
usedplots(index) = 1;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function computing spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [eegspecdB, freqs, specstd] = spectcomp( data, frames, srate, epoch_subset, g, newweights);
usepwelch = license('checkout','Signal_Toolbox') && exist('pwelch');
if exist('newweights') == 1
nchans = size(newweights,1);
else
nchans = size(data,1);
end
%fftlength = 2^round(log(srate)/log(2))*g.freqfac;
if isempty(g.winsize)
% winlength = max(pow2(nextpow2(frames)-3),4); %*2 since diveded by 2 later
% winlength = min(winlength, 512);
% winlength = max(winlength, 256);
winlength = min(round(srate), frames);
else
winlength = g.winsize;
end
if strcmpi(g.wintype,'blackmanharris')
if usepwelch
winlength = blackmanharris(round(winlength/g.blckhn));
else
g.wintype = 'hamming';
fprintf('\nSignal processing toolbox (SPT) absent: unable to use Blackman-Harris window\n');
fprintf(' using pwelch function from Octave\n');
end
end
if isempty(g.nfft) && strcmp(g.wintype,'hamming')
%fftlength = 2^(nextpow2(winlength))*g.freqfac;
fftlength = winlength*g.freqfac;
elseif ~isempty(g.nfft) && strcmp(g.wintype,'hamming')
fftlength = g.nfft;
elseif strcmp(g.wintype,'blackmanharris')
fftlength = 2^(nextpow2(length(winlength)))*g.freqfac;
end
if ~usepwelch
myfprintf(g.verbose, '\nSignal processing toolbox (SPT) absent: spectrum computed using the pwelch()\n');
myfprintf(g.verbose, 'function from Octave which is supposedly 100%% compatible with the Matlab pwelch function\n');
end
myfprintf(g.verbose,' (window length %d; fft length: %d; overlap %d):\n', winlength, fftlength, g.overlap);
for c=1:nchans % scan channels or components
if exist('newweights') == 1
if strcmp(g.icamode, 'normal')
tmpdata = newweights(c,:)*data; % component activity
else % data - component contribution
tmpdata = data(g.plotchan,:) - (g.icawinv(g.plotchan,c)*newweights(c,:))*data;
end
else
tmpdata = data(c,:); % channel activity
end
if strcmp(g.reref, 'averef')
tmpdata = averef(tmpdata);
end
for e=epoch_subset
if isempty(g.boundaries)
if usepwelch
[tmpspec,freqs] = pwelch(matsel(tmpdata,frames,0,1,e),...
winlength,g.overlap,fftlength,srate);
else
[tmpspec,freqs] = spec(matsel(tmpdata,frames,0,1,e),fftlength,srate,...
winlength,g.overlap);
end
%[tmpspec,freqs] = psd(matsel(tmpdata,frames,0,1,e),fftlength,srate,...
% winlength,g.overlap);
if c==1 && e==epoch_subset(1)
eegspec = zeros(nchans,length(freqs));
specstd = zeros(nchans,length(freqs));
end
eegspec(c,:) = eegspec(c,:) + tmpspec';
specstd(c,:) = specstd(c,:) + tmpspec'.^2;
else
g.boundaries = round(g.boundaries);
for n=1:length(g.boundaries)-1
if g.boundaries(n+1) - g.boundaries(n) >= winlength % ignore segments of less than winlength
if usepwelch
[tmpspec,freqs] = pwelch(tmpdata(e,g.boundaries(n)+1:g.boundaries(n+1)),...
winlength,g.overlap,fftlength,srate);
else
[tmpspec,freqs] = spec(tmpdata(e,g.boundaries(n)+1:g.boundaries(n+1)),...
fftlength,srate,winlength,g.overlap);
end
if exist('eegspec') ~= 1
eegspec = zeros(nchans,length(freqs));
specstd = zeros(nchans,length(freqs));
end
eegspec(c,:) = eegspec(c,:) + tmpspec'* ...
((g.boundaries(n+1)-g.boundaries(n)+1)/g.boundaries(end));
specstd(c,:) = eegspec(c,:) + tmpspec'.^2 * ...
((g.boundaries(n+1)-g.boundaries(n)+1)/g.boundaries(end));
end
end
end
end
myfprintf(g.verbose,'.');
end
n = length(epoch_subset);
eegspecdB = eegspec/n; % normalize by the number of sections
if n>1 % normalize standard deviation by the number of sections
specstd = sqrt( (specstd + eegspec.^2/n)/(n-1) );
else specstd = [];
end
return;
function myfprintf(verbose, varargin)
if strcmpi(verbose, 'on')
fprintf(varargin{:});
end
Datasets
Models
