% TIMEF - Returns estimates and plots of mean event-related spectral
% perturbation (ERSP) and inter-trial coherence (ITC) changes
% across event-related trials (epochs) of a single input time series.
% * Uses either fixed-window, zero-padded FFTs (fastest), wavelet
% 0-padded DFTs (both Hanning-tapered), OR multitaper spectra ('mtaper').
% * For the wavelet and FFT methods, output frequency spacing
% is the lowest frequency ('srate'/'winsize') divided by 'padratio'.
% NaN input values (such as returned by EVENTLOCK) are ignored.
% * If 'alpha' is given, then bootstrap statistics are computed
% (from a distribution of 'naccu' surrogate data trials) and
% non-significant features of the output plots are zeroed out
% (i.e., plotted in green).
% * Given a 'topovec' scalp map weights vector and an 'elocs' electrode
% location file or structure, the figure also shows a TOPOPLOT
% image of the specified scalp map.
%
% * Note: Left-click on subplots to view and zoom in separate windows.
% Usage:
% >> [ersp,itc,powbase,times,freqs,erspboot,itcboot,itcphase] = ...
% timef(data,frames,tlimits,srate,cycles,...
% 'key1',value1,'key2',value2, ... );
% NOTE:
% * For more detailed information about TIMEF, >> timef details
% * Default values may differ when called from POP_TIMEF
%
% Required inputs:
% data = Single-channel data vector (1,frames*ntrials) (required)
% frames = Frames per trial {def|[]: datalength}
% tlimits = [mintime maxtime] (ms) Epoch time limits
% {def|[]: from frames,srate}
% srate = data sampling rate (Hz) {def:250}
% cycles = If 0 -> Use FFTs (with constant window length) {0 = FFT}
% If >0 -> Number of cycles in each analysis wavelet
% If [wavecycles factor] -> wavelet cycles increase with
% frequency beginning at wavecyles (0<factor<1; factor=1
% -> no increase, standard wavelets; factor=0 -> fixed epoch
% length, as in FFT. Else, 'mtaper' -> multitaper decomp.
%
% Optional Inter-Irial Coherence (ITC) type:
% 'type' = ['coher'|'phasecoher'] Compute either linear coherence
% ('coher') or phase coherence ('phasecoher') also known
% as the phase coupling factor {'phasecoher'}.
% Optional detrending:
% 'detret' = ['on'|'off'], Detrend data in time. {'off'}
% 'detrep' = ['on'|'off'], Detrend data across trials {'off'}
%
% Optional FFT/DFT parameters:
% 'winsize' = If cycles==0: data subwindow length (fastest, 2^n<frames);
% If cycles >0: *longest* window length to use. This
% determines the lowest output frequency {~frames/8}
% 'timesout' = Number of output times (int<frames-winframes) {200}
% 'padratio' = FFT-length/winframes (2^k) {2}
% Multiplies the number of output frequencies by
% dividing their spacing. When cycles==0, frequency
% spacing is (low_freq/padratio).
% 'maxfreq' = Maximum frequency (Hz) to plot (& to output if cycles>0)
% If cycles==0, all FFT frequencies are output. {50}
% 'baseline' = Spectral baseline window center end-time (in ms). {0}
% 'powbase' = Baseline spectrum (power, not dB) to normalize the data.
% {def|NaN->from data}
%
% Optional multitaper parameters:
% 'mtaper' = If [N W], performs multitaper decomposition.
% (N is the time resolution and W the frequency resolution;
% maximum taper number is 2NW-1). Overwrites 'winsize' and
% 'padratio'.
% If [N W K], uses K Slepian tapers (if possible).
% Phase is calculated using standard methods.
% The use of mutitaper with wavelets (cycles>0) is not
% recommended (as multiwavelets are not implemented).
% Uses Matlab functions DPSS, PMTM. {no multitaper}
%
% Optional bootstrap parameters:
% 'alpha' = If non-0, compute two-tailed bootstrap significance prob.
% level. Show non-signif. output values in green {0}
% 'naccu' = Number of bootstrap replications to accumulate {200}
% 'baseboot' = Bootstrap baseline to subtract (1 -> use 'baseline'(above)
% 0 -> use whole trial) {1}
% Optional scalp map:
% 'topovec' = Scalp topography (map) to plot {none}
% 'elocs' = Electrode location file for scalp map
% File should be ascii in format of >> topoplot example
% May also be an EEG.chanlocs struct.
% {default: file named in icadefs.m}
% Optional plotting parameters:
% 'hzdir' = ['up'|'down'] Direction of the frequency axes; reads default
% from icadefs.m {'up'}
% 'plotersp' = ['on'|'off'] Plot power spectral perturbations {'on'}
% 'plotitc' = ['on'|'off'] Plot inter trial coherence {'on'}
% 'plotphase' = ['on'|'off'] Plot sign of the phase in the ITC panel, i.e.
% green->red, pos.-phase ITC, green->blue, neg.-phase ITC {'on'}
% 'erspmax' = [real dB] set the ERSP max. for the scale (min= -max){auto}
% 'itcmax' = [real<=1] set the ITC maximum for the scale {auto}
% 'title' = Optional figure title {none}
% 'marktimes' = Non-0 times to mark with a dotted vertical line (ms) {none}
% 'linewidth' = Line width for 'marktimes' traces (thick=2, thin=1) {2}
% 'pboot' = Bootstrap power limits (e.g., from TIMEF) {from data}
% 'rboot' = Bootstrap ITC limits (e.g., from TIMEF) {from data}
% 'axesfont' = Axes text font size {10}
% 'titlefont' = Title text font size {8}
% 'vert' = [times_vector] -> plot vertical dashed lines at given ms.
% 'verbose' = ['on'|'off'] print text {'on'}
%
% Outputs:
% ersp = Matrix (nfreqs,timesout) of log spectral diffs. from
% baseline (in dB). NB: Not masked for significance.
% Must do this using erspboot
% itc = Matrix of inter-trial coherencies (nfreqs,timesout)
% (range: [0 1]) NB: Not masked for significance.
% Must do this using itcboot
% powbase = Baseline power spectrum (NOT in dB, used to norm. the ERSP)
% times = Vector of output times (sub-window centers) (in ms)
% freqs = Vector of frequency bin centers (in Hz)
% erspboot = Matrix (2,nfreqs) of [lower;upper] ERSP significance diffs
% itcboot = Matrix (2,nfreqs) of [lower;upper] ITC thresholds (not diffs)
% itcphase = Matrix (nfreqs,timesout) of ITC phase (in radians)
%
% Plot description:
% Assuming both 'plotersp' and 'plotitc' options are 'on' (= default).
% The upper panel presents the data ERSP (Event-Related Spectral Perturbation)
% in dB, with mean baseline spectral activity (in dB) subtracted. Use
% "'baseline', NaN" to prevent TIMEF from removing the baseline.
% The lower panel presents the data ITC (Inter-Trial Coherence).
% Click on any plot axes to pop up a new window (using 'AXCOPY')
% -- Upper left marginal panel presents the mean spectrum during the baseline
% period (blue), and when significance is set, the significance threshold
% at each frequency (dotted green-black trace).
% -- The marginal panel under the ERSP image shows the maximum (green) and
% minimum (blue) ERSP values relative to baseline power at each frequency.
% -- The lower left marginal panel shows mean ITC across the imaged time range
% (blue), and when significance is set, the significance threshold (dotted
% green-black).
% -- The marginal panel under the ITC image shows the ERP (which is produced by
% ITC across the data spectral pass band).
%
% Author: Sigurd Enghoff, Arnaud Delorme & Scott Makeig
% CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
% Known problems:
% Significance masking currently fails for linear coherence.
%
% See also: CROSSF
% Copyright (C) 1998 Sigurd Enghoff, Scott Makeig, Arnaud Delorme,
% CNL / Salk Institute 8/1/98-8/28/01
%
% 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.
% 10-19-98 avoided division by zero (using MIN_ABS) -sm
% 10-19-98 improved usage message and commandline info printing -sm
% 10-19-98 made valid [] values for tvec and g.elocs -sm
% 04-01-99 added missing freq in freqs and plots, fixed log scaling bug -se & -tpj
% 06-29-99 fixed frequency indexing for constant-Q -se
% 08-24-99 reworked to handle NaN input values -sm
% 12-07-99 adjusted ERPtimes to plot ERP under ITC -sm
% 12-22-99 debugged ERPtimes, added BASE_BOOT -sm
% 01-10-00 debugged BASE_BOOT=0 -sm
% 02-28-00 added NOTE on formula derivation below -sm
% 03-16-00 added AXCOPY feature -sm & tpj
% 04-16-00 added multiple marktimes loop -sm
% 04-20-00 fixed ITC cbar limits when specified in input -sm
% 07-29-00 changed frequencies displayed msg -sm
% 10-12-00 fixed bug in freqs when cycles>0 -sm
% 02-07-01 fixed inconsistency in BASE_BOOT use -sm
% 08-28-01 matlab 'key' value arguments -ad
% 08-28-01 multitaper decomposition -ad
% 01-25-02 reformated help & license -ad
% 03-08-02 debug & compare to old timef function -ad
% 03-16-02 timeout automatically adjusted if too high -ad
% 04-02-02 added 'coher' option -ad
function [P,R,mbase,times,freqs,Pboot,Rboot,Rphase,PA] = timef(X,frames,tlimits,Fs,varwin,varargin);
% Note: undocumented arg PA is output of 'phsamp','on'
%varwin,winsize,g.timesout,g.padratio,g.maxfreq,g.topovec,g.elocs,g.alpha,g.marktimes,g.powbase,g.pboot,g.rboot)
% ITC: Normally, R = |Sum(Pxy)| / (Sum(|Pxx|)*Sum(|Pyy|)) is linear coherence.
% But here, we consider: Phase(Pyy) = 0 and |Pyy| = 1 -> Pxy = Pxx
% Giving, R = |Sum(Pxx)|/Sum(|Pxx|), the inter-trial coherence (ITC)
% Also called 'phase-locking factor' by Tallon-Baudry et al. (1996),
% the ITC is the phase coherence between the data time series and the
% time-locking event time series.
% Read system-wide / dir-wide constants:
icadefs
% Constants set here:
ERSP_CAXIS_LIMIT = 0; % 0 -> use data limits; else positive value
% giving symmetric +/- caxis limits.
ITC_CAXIS_LIMIT = 0; % 0 -> use data limits; else positive value
% giving symmetric +/- caxis limits.
% Commandline arg defaults:
DEFAULT_EPOCH = NaN; % Frames per trial
DEFAULT_TIMLIM = NaN; % Time range of g.frames (ms)
DEFAULT_FS = 250; % Sampling frequency (Hz)
DEFAULT_NWIN = 200; % Number of windows = horizontal resolution
DEFAULT_VARWIN = 0; % Fixed window length or fixed number of cycles.
% =0: fix window length to that determined by nwin
% >0: set window length equal to varwin cycles
% Bounded above by winsize, which determines
% the min. freq. to be computed.
DEFAULT_OVERSMP = 2; % Number of times to oversample frequencies
DEFAULT_MAXFREQ = 50; % Maximum frequency to display (Hz)
DEFAULT_TITLE = ''; % Figure title
DEFAULT_ELOC = 'chan.locs'; % Channel location file
DEFAULT_ALPHA = NaN; % Percentile of bins to keep
DEFAULT_MARKTIME= NaN;
% Font sizes:
AXES_FONT = 10; % axes text FontSize
TITLE_FONT = 8;
if nargout>7
Rphase = []; % initialize in case Rphase asked for, but ITC not computed
end
if (nargin < 1)
help timef
return
end
if ischar(X) && strcmp(X,'details')
more on
help timefdetails
more off
return
end
if (min(size(X))~=1 || length(X)<2)
error('Data must be a row or column vector.');
end
if nargin < 2 || isempty(frames) || isnan(frames)
frames = DEFAULT_EPOCH;
elseif (~isnumeric(frames) || length(frames)~=1 || frames~=round(frames))
error('Value of frames must be an integer.');
elseif (frames <= 0)
error('Value of frames must be positive.');
elseif (rem(length(X),frames) ~= 0)
error('Length of data vector must be divisible by frames.');
end
if isnan(frames) || isempty(frames)
frames = length(X);
end
if nargin < 3 || isempty(tlimits) || isnan(tlimits(1))
tlimits = DEFAULT_TIMLIM;
elseif (~isnumeric(tlimits) || sum(size(tlimits))~=3)
error('Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
error('tlimits interval must be ascending.');
end
if (nargin < 4)
Fs = DEFAULT_FS;
elseif (~isnumeric(Fs) || length(Fs)~=1)
error('Value of srate must be a number.');
elseif (Fs <= 0)
error('Value of srate must be positive.');
end
if isempty(tlimits) || isnan(tlimits(1))
hlim = 1000*frames/Fs; % fit default tlimits to srate and frames
tlimits = [0 hlim];
end
framesdiff = frames - Fs*(tlimits(2)-tlimits(1))/1000;
if abs(framesdiff) > 1
error('Given time limits, frames and sampling rate are incompatible');
elseif framesdiff ~= 0
tlimits(1) = tlimits(1) - 0.5*framesdiff*1000/Fs;
tlimits(2) = tlimits(2) + 0.5*framesdiff*1000/Fs;
fprintf('Adjusted time limits slightly, to [%.1f,%.1f] ms, to match frames and srate.\n',tlimits(1),tlimits(2));
end
if (nargin < 5)
varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) || length(varwin)>2)
error('Value of cycles must be a number.');
elseif (varwin < 0)
error('Value of cycles must be zero or positive.');
end
% consider structure for these arguments
% --------------------------------------
if ~isempty(varargin)
try, g = struct(varargin{:});
catch, error('Argument error in the {''param'', value} sequence'); end;
end
g.tlimits = tlimits;
g.frames = frames;
g.srate = Fs;
g.cycles = varwin(1);
if length(varwin)>1
g.cyclesfact = varwin(2);
else
g.cyclesfact = 1;
end
try, g.title; catch, g.title = DEFAULT_TITLE; end
try, g.winsize; catch, g.winsize = max(pow2(nextpow2(g.frames)-3),4); end
try, g.pad; catch, g.pad = max(pow2(nextpow2(g.winsize)),4); end
try, g.timesout; catch, g.timesout = DEFAULT_NWIN; end
try, g.padratio; catch, g.padratio = DEFAULT_OVERSMP; end
try, g.maxfreq; catch, g.maxfreq = DEFAULT_MAXFREQ; end
try, g.topovec; catch, g.topovec = []; end
try, g.elocs; catch, g.elocs = DEFAULT_ELOC; end
try, g.alpha; catch, g.alpha = DEFAULT_ALPHA; end;
try, g.marktimes; catch, g.marktimes = DEFAULT_MARKTIME; end
try, g.powbase; catch, g.powbase = NaN; end
try, g.pboot; catch, g.pboot = NaN; end
try, g.rboot; catch, g.rboot = NaN; end
try, g.plotersp; catch, g.plotersp = 'on'; end
try, g.plotitc; catch, g.plotitc = 'on'; end
try, g.detrep; catch, g.detrep = 'off'; end
try, g.detret; catch, g.detret = 'off'; end
try, g.baseline; catch, g.baseline = 0; end
try, g.baseboot; catch, g.baseboot = 1; end
try, g.linewidth; catch, g.linewidth = 2; end
try, g.naccu; catch, g.naccu = 200; end
try, g.mtaper; catch, g.mtaper = []; end
try, g.vert; catch, g.vert = []; end
try, g.type; catch, g.type = 'phasecoher'; end
try, g.phsamp; catch, g.phsamp = 'off'; end
try, g.plotphase; catch, g.plotphase = 'on'; end
try, g.itcmax; catch, g.itcmax = []; end
try, g.erspmax; catch, g.erspmax = []; end
try, g.verbose; catch, g.verbose = 'on'; end
try, g.chaninfo; catch, g.chaninfo = []; end
try, g.hzdir; catch, g.hzdir = HZDIR; end; % default from icadefs
lasterr('');
% testing arguments consistency
% -----------------------------
if strcmp(g.hzdir,'up')
g.hzdir = 'normal';
elseif strcmp(g.hzdir,'down')
g.hzdir = 'reverse';
else
error('unknown ''hzdir'' value - not ''up'' or ''down''');
end
switch lower(g.verbose)
case { 'on', 'off' }, ;
otherwise error('verbose must be either on or off');
end
if (~ischar(g.title))
error('Title must be a string.');
end
if (~isnumeric(g.winsize) || length(g.winsize)~=1 || g.winsize~=round(g.winsize))
error('Value of winsize must be an integer number.');
elseif (g.winsize <= 0)
error('Value of winsize must be positive.');
elseif (g.cycles == 0 && pow2(nextpow2(g.winsize)) ~= g.winsize)
error('Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (g.winsize > g.frames)
error('Value of winsize must be less than frames per epoch.');
end
if (~isnumeric(g.timesout) || length(g.timesout)~=1 || g.timesout~=round(g.timesout))
error('Value of timesout must be an integer number.');
elseif (g.timesout <= 0)
error('Value of timesout must be positive.');
end
if (g.timesout > g.frames-g.winsize)
g.timesout = g.frames-g.winsize;
disp(['Value of timesout must be <= frames-winsize, timeout adjusted to ' int2str(g.timesout) ]);
end
if (~isnumeric(g.padratio) || length(g.padratio)~=1 || g.padratio~=round(g.padratio))
error('Value of padratio must be an integer.');
elseif (g.padratio <= 0)
error('Value of padratio must be positive.');
elseif (pow2(nextpow2(g.padratio)) ~= g.padratio)
error('Value of padratio must be an integer power of two [1,2,4,8,16,...]');
end
if (~isnumeric(g.maxfreq) || length(g.maxfreq)~=1)
error('Value of maxfreq must be a number.');
elseif (g.maxfreq <= 0)
error('Value of maxfreq must be positive.');
elseif (g.maxfreq > Fs/2)
myprintf(g.verbose,['Warning: value of maxfreq reduced to Nyquist rate' ...
' (%3.2f)\n\n'], Fs/2);
g.maxfreq = Fs/2;
end
if isempty(g.topovec)
g.topovec = [];
if isempty(g.elocs)
error('Channel location file must be specified.');
end
end
if isempty(g.elocs)
g.elocs = DEFAULT_ELOC;
elseif (~ischar(g.elocs)) && ~isstruct(g.elocs)
error('Channel location file must be a valid text file.');
end
if (~isnumeric(g.alpha) || length(g.alpha)~=1)
error('timef(): Value of g.alpha must be a number.\n');
elseif (round(g.naccu*g.alpha) < 2)
myprintf(g.verbose,'Value of g.alpha is out of the normal range [%g,0.5]\n',2/g.naccu);
g.naccu = round(2/g.alpha);
myprintf(g.verbose,' Increasing the number of bootstrap iterations to %d\n',g.naccu);
end
if g.alpha>0.5 || g.alpha<=0
error('Value of g.alpha is out of the allowed range (0.00,0.5).');
end
if ~isnan(g.alpha)
if g.baseboot > 0
myprintf(g.verbose,'Bootstrap analysis will use data in baseline (pre-0 centered) subwindows only.\n')
else
myprintf(g.verbose,'Bootstrap analysis will use data in all subwindows.\n')
end
end
if ~isnumeric(g.vert)
error('vertical line(s) option must be a vector');
else
if ~isempty(g.vert)
if min(g.vert(:)) < g.tlimits(1) || max(g.vert(:)) > g.tlimits(2)
error('vertical line(s) time out-of-bound');
end
end
end
if ~isnan (g.rboot)
if size(g.rboot) == [1,1]
if g.cycles == 0
g.rboot = g.rboot*ones(g.winsize*g.padratio/2);
end
end
end
if ~isempty(g.mtaper) % mutitaper, inspired from Bijan Pesaran matlab function
if length(g.mtaper) < 3
%error('mtaper argument must be [N W] or [N W K]');
if g.mtaper(1) * g.mtaper(2) < 1
error('mtaper 2 first arguments'' product must be higher than 1');
end
if length(g.mtaper) == 2
g.mtaper(3) = floor( 2*g.mtaper(2)*g.mtaper(1) - 1);
end
if length(g.mtaper) == 3
if g.mtaper(3) > 2 * g.mtaper(1) * g.mtaper(2) -1
error('mtaper number too high (maximum (2*N*W-1))');
end
end
disp(['Using ' num2str(g.mtaper(3)) ' tapers.']);
NW = g.mtaper(1)*g.mtaper(2); % product NW
N = g.mtaper(1)*g.srate;
[e,v] = dpss(N, NW, 'calc');
e=e(:,1:g.mtaper(3));
g.alltapers = e;
else
g.alltapers = g.mtaper;
disp('mtaper argument not [N W] or [N W K]; considering raw taper matrix');
end
g.winsize = size(g.alltapers, 1);
g.pad = max(pow2(nextpow2(g.winsize)),256); % pad*nextpow
%nfk = floor([0 g.maxfreq]./g.srate.*g.pad); % not used any more
%g.padratio = 2*nfk(2)/g.winsize;
g.padratio = g.pad/g.winsize;
%compute number of frequencies
%nf = max(256, g.pad*2^nextpow2(g.winsize+1));
%nfk = floor([0 g.maxfreq]./g.srate.*nf);
%freqs = linspace( 0, g.maxfreq, diff(nfk)); % this also work in the case of a FFT
end;
switch lower(g.plotphase)
case { 'on', 'off' }, ;
otherwise error('plotphase must be either on or off');
end
switch lower(g.plotersp)
case { 'on', 'off' }, ;
otherwise error('plotersp must be either on or off');
end
switch lower(g.plotitc)
case { 'on', 'off' }, ;
otherwise error('plotitc must be either on or off');
end
switch lower(g.detrep)
case { 'on', 'off' }, ;
otherwise error('detrep must be either on or off');
end
switch lower(g.detret)
case { 'on', 'off' }, ;
otherwise error('detret must be either on or off');
end
switch lower(g.phsamp)
case { 'on', 'off' }, ;
otherwise error('phsamp must be either on or off');
end
if ~isnumeric(g.linewidth)
error('linewidth must be numeric');
end
if ~isnumeric(g.naccu)
error('naccu must be numeric');
end
if ~isnumeric(g.baseline)
error('baseline must be numeric');
end
switch g.baseboot
case {0,1}, ;
otherwise, error('baseboot must be 0 or 1');
end
switch g.type
case { 'coher', 'phasecoher', 'phasecoher2' },;
otherwise error('Type must be either ''coher'' or ''phasecoher''');
end;
if isnan(g.baseline)
g.unitpower = 'uV/Hz';
else
g.unitpower = 'dB';
end
if (g.cycles == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
freqs = linspace(0, g.srate/2, g.padratio*g.winsize/2+1);
freqs = freqs(2:end);
win = hanning(g.winsize);
P = zeros(g.padratio*g.winsize/2,g.timesout); % summed power
PP = zeros(g.padratio*g.winsize/2,g.timesout); % power
R = zeros(g.padratio*g.winsize/2,g.timesout); % mean coherence
RR = zeros(g.padratio*g.winsize/2,g.timesout); % (coherence)
Pboot = zeros(g.padratio*g.winsize/2,g.naccu); % summed bootstrap power
Rboot = zeros(g.padratio*g.winsize/2,g.naccu); % summed bootstrap coher
Rn = zeros(1,g.timesout);
Rbn = 0;
switch g.type
case { 'coher' 'phasecoher2' },
cumulX = zeros(g.padratio*g.winsize/2,g.timesout);
cumulXboot = zeros(g.padratio*g.winsize/2,g.naccu);
case 'phasecoher'
switch g.phsamp
case 'on'
cumulX = zeros(g.padratio*g.winsize/2,g.timesout);
end
end;
else % %%%%%%%%%%%%%%%%%% cycles>0, Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%
freqs = g.srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
dispf = find(freqs <= g.maxfreq);
freqs = freqs(dispf);
win = dftfilt(g.winsize,g.maxfreq/g.srate,g.cycles,g.padratio,g.cyclesfact);
P = zeros(size(win,2),g.timesout); % summed power
R = zeros(size(win,2),g.timesout); % mean coherence
PP = repmat(NaN,size(win,2),g.timesout); % initialize with NaN
RR = repmat(NaN,size(win,2),g.timesout); % initialize with NaN
Pboot = zeros(size(win,2),g.naccu); % summed bootstrap power
Rboot = zeros(size(win,2),g.naccu); % summed bootstrap coher
Rn = zeros(1,g.timesout);
Rbn = 0;
switch g.type
case { 'coher' 'phasecoher2' },
cumulX = zeros(size(win,2),g.timesout);
cumulXboot = zeros(size(win,2),g.naccu);
case 'phasecoher'
switch g.phsamp
case 'on'
cumulX = zeros(size(win,2),g.timesout);
end
end;
end
switch g.phsamp
case 'on'
PA = zeros(size(P,1),size(P,1),g.timesout); % NB: (freqs,freqs,times)
end % phs amp
wintime = 1000/g.srate*(g.winsize/2); % (1000/g.srate)*(g.winsize/2);
times = [g.tlimits(1)+wintime:(g.tlimits(2)-g.tlimits(1)-2*wintime)/(g.timesout-1):g.tlimits(2)-wintime];
ERPtimes = [g.tlimits(1):(g.tlimits(2)-g.tlimits(1))/(g.frames-1):g.tlimits(2)+0.000001];
ERPindices = [];
for ti=times
[tmp indx] = min(abs(ERPtimes-ti));
ERPindices = [ERPindices indx];
end
ERPtimes = ERPtimes(ERPindices); % subset of ERP frames on t/f window centers
if ~isempty(find(times < g.baseline))
baseln = find(times < g.baseline); % subtract means of pre-0 (centered) windows
else
baseln = 1:length(times); % use all times as baseline
end
if ~isnan(g.alpha) && length(baseln)==0
myprintf(g.verbose,'timef(): no window centers in baseline (times<%g) - shorten (max) window length.\n', g.baseline)
return
elseif ~isnan(g.alpha) && g.baseboot
myprintf(g.verbose,' %d bootstrap windows in baseline (center times < %g).\n',...
length(baseln), g.baseline)
end
dispf = find(freqs <= g.maxfreq);
stp = (g.frames-g.winsize)/(g.timesout-1);
myprintf(g.verbose,'Computing Event-Related Spectral Perturbation (ERSP) and\n');
switch g.type
case 'phasecoher', myprintf(g.verbose,' Inter-Trial Phase Coherence (ITC) images based on %d trials\n',length(X)/g.frames);
case 'phasecoher2', myprintf(g.verbose,' Inter-Trial Phase Coherence 2 (ITC) images based on %d trials\n',length(X)/g.frames);
case 'coher', myprintf(g.verbose,' Linear Inter-Trial Coherence (ITC) images based on %d trials\n',length(X)/g.frames);
end
myprintf(g.verbose,' of %d frames sampled at %g Hz.\n',g.frames,g.srate);
myprintf(g.verbose,'Each trial contains samples from %d ms before to\n',g.tlimits(1));
myprintf(g.verbose,' %.0f ms after the timelocking event.\n',g.tlimits(2));
myprintf(g.verbose,'The window size used is %d samples (%g ms) wide.\n',g.winsize,2*wintime);
myprintf(g.verbose,'The window is applied %d times at an average step\n',g.timesout);
myprintf(g.verbose,' size of %g samples (%g ms).\n',stp,1000*stp/g.srate);
myprintf(g.verbose,'Results are oversampled %d times; the %d frequencies\n',g.padratio,length(dispf));
myprintf(g.verbose,' displayed are from %2.1f Hz to %3.1f Hz.\n',freqs(dispf(1)),freqs(dispf(end)));
if ~isnan(g.alpha)
myprintf(g.verbose,'Only significant values (bootstrap p<%g) will be colored;\n',g.alpha)
myprintf(g.verbose,' non-significant values will be plotted in green\n');
end
trials = length(X)/g.frames;
baselength = length(baseln);
myprintf(g.verbose,'\nOf %d trials total, processing trial:',trials);
% detrend over epochs (trials) if requested
% -----------------------------------------
switch g.detrep
case 'on'
X = reshape(X, g.frames, length(X)/g.frames);
X = X - mean(X,2)*ones(1, length(X(:))/g.frames);
X = X(:)';
end;
for i=1:trials
if (rem(i,100)==0)
myprintf(g.verbose,'\n');
end
if (rem(i,10) == 0)
myprintf(g.verbose,'%d',i);
elseif (rem(i,2) == 0)
myprintf(g.verbose,'.');
end
ERP = blockave(X,g.frames); % compute the ERP trial average
Wn = zeros(1,g.timesout);
for j=1:g.timesout,
tmpX = X([1:g.winsize]+floor((j-1)*stp)+(i-1)*g.frames);
% pull out data g.frames
tmpX = tmpX - mean(tmpX); % remove the mean for that window
switch g.detret, case 'on', tmpX = detrend(tmpX); end
if ~any(isnan(tmpX))
if (g.cycles == 0) % FFT
if ~isempty(g.mtaper) % apply multitaper (no hanning window)
tmpXMT = fft(g.alltapers .* ...
(tmpX(:) * ones(1,size(g.alltapers,2))), g.pad);
%tmpXMT = tmpXMT(nfk(1)+1:nfk(2),:);
tmpXMT = tmpXMT(2:g.padratio*g.winsize/2+1,:);
PP(:,j) = mean(abs(tmpXMT).^2, 2);
% power; can also ponderate multitaper by their eigenvalues v
tmpX = win .* tmpX(:);
tmpX = fft(tmpX, g.pad);
tmpX = tmpX(2:g.padratio*g.winsize/2+1);
else
% TF and MC (12/2006): Calculation changes made so that
% power can be correctly calculated from ERSP.
tmpX = win .* tmpX(:);
tmpX = fft(tmpX,g.padratio*g.winsize);
tmpX = tmpX / g.winsize; % TF and MC (12/11/2006): normalization, divide by g.winsize
tmpX = tmpX(2:g.padratio*g.winsize/2+1);
PP(:,j) = 2/0.375*abs(tmpX).^2; % power
% TF and MC (12/14/2006): multiply by 2 account for negative frequencies,
% Counteract the reduction by a factor 0.375
% that occurs as a result of cosine (Hann) tapering. Refer to Bug 446
end;
else % wavelet
if ~isempty(g.mtaper) % apply multitaper
tmpXMT = g.alltapers .* (tmpX(:) * ones(1,size(g.alltapers,2)));
tmpXMT = transpose(win) * tmpXMT;
PP(:,j) = mean(abs(tmpXMT).^2, 2); % power
tmpX = transpose(win) * tmpX(:);
else
tmpX = transpose(win) * tmpX(:);
PP(:,j) = abs(tmpX).^2; % power
end
end
if abs(tmpX) < eps % If less than smallest possible machine value
% (i.e. if it's zero) then call it 0.
RR(:,j) = zeros(size(RR(:,j)));
else
switch g.type
case { 'coher' },
RR(:,j) = tmpX;
cumulX(:,j) = cumulX(:,j)+abs(tmpX).^2;
case { 'phasecoher2' },
RR(:,j) = tmpX;
cumulX(:,j) = cumulX(:,j)+abs(tmpX);
case 'phasecoher',
RR(:,j) = tmpX ./ abs(tmpX); % normalized cross-spectral vector
switch g.phsamp
case 'on'
cumulX(:,j) = cumulX(:,j)+abs(tmpX); % accumulate for PA
end
end
end
Wn(j) = 1;
end
switch g.phsamp
case 'on' % PA (freq x freq x time)
PA(:,:,j) = PA(:,:,j) + (tmpX ./ abs(tmpX)) * ((PP(:,j)))';
% cross-product: unit phase (column)
% times amplitude (row)
end
end % window
if ~isnan(g.alpha) % save surrogate data for bootstrap analysis
j = 1;
goodbasewins = find(Wn==1);
if g.baseboot % use baseline windows only
goodbasewins = find(goodbasewins<=baselength);
end
ngdbasewins = length(goodbasewins);
if ngdbasewins>1
while j <= g.naccu
i=ceil(rand*ngdbasewins);
i=goodbasewins(i);
Pboot(:,j) = Pboot(:,j) + PP(:,i);
Rboot(:,j) = Rboot(:,j) + RR(:,i);
switch g.type
case 'coher', cumulXboot(:,j) = cumulXboot(:,j)+abs(tmpX).^2;
case 'phasecoher2', cumulXboot(:,j) = cumulXboot(:,j)+abs(tmpX);
end
j = j+1;
end
Rbn = Rbn + 1;
end
end % bootstrap
Wn = find(Wn>0);
if length(Wn)>0
P(:,Wn) = P(:,Wn) + PP(:,Wn); % add non-NaN windows
R(:,Wn) = R(:,Wn) + RR(:,Wn);
Rn(Wn) = Rn(Wn) + ones(1,length(Wn)); % count number of addends
end
end % trial
% if coherence, perform the division
% ----------------------------------
switch g.type
case 'coher',
R = R ./ ( sqrt( trials*cumulX ) );
if ~isnan(g.alpha)
Rboot = Rboot ./ ( sqrt( trials*cumulXboot ) );
end
case 'phasecoher2',
R = R ./ ( cumulX );
if ~isnan(g.alpha)
Rboot = Rboot ./ cumulXboot;
end;
case 'phasecoher',
R = R ./ (ones(size(R,1),1)*Rn);
end;
switch g.phsamp
case 'on'
tmpcx(1,:,:) = cumulX; % allow ./ below
for j=1:g.timesout
PA(:,:,j) = PA(:,:,j) ./ repmat(PP(:,j)', [size(PP,1) 1]);
end
end
if min(Rn) < 1
myprintf(g.verbose,'timef(): No valid timef estimates for windows %s of %d.\n',...
int2str(find(Rn==0)),length(Rn));
Rn(find(Rn<1))==1;
return
end
P = P ./ (ones(size(P,1),1) * Rn);
if isnan(g.powbase)
myprintf(g.verbose,'\nComputing the mean baseline spectrum\n');
mbase = mean(P(:,baseln),2)';
else
myprintf(g.verbose,'Using the input baseline spectrum\n');
mbase = g.powbase;
end
if ~isnan( g.baseline(1) ) && ~isnan( mbase(1) )
P = 10 * (log10(P) - repmat(log10(mbase(1:size(P,1)))',[1 g.timesout])); % convert to (10log10) dB
else
P = 10 * log10(P);
end
Rsign = sign(imag(R));
if nargout > 7
for lp = 1:size(R,1)
Rphase(lp,:) = rem(angle(R(lp,:)),2*pi); % replaced obsolete phase() -sm 2/1/6
end
Rphase(find(Rphase>pi)) = 2*pi-Rphase(find(Rphase>pi));
Rphase(find(Rphase<-pi)) = -2*pi-Rphase(find(Rphase<-pi));
end
R = abs(R); % convert coherence vector to magnitude
if ~isnan(g.alpha) % if bootstrap analysis included . . .
if Rbn>0
i = round(g.naccu*g.alpha);
if isnan(g.pboot)
Pboot = Pboot / Rbn; % normalize
if ~isnan( g.baseline )
Pboot = 10 * (log10(Pboot) - repmat(log10(mbase)',[1 g.naccu]));
else
Pboot = 10 * log10(Pboot);
end;
Pboot = sort(Pboot');
Pboot = [mean(Pboot(1:i,:)) ; mean(Pboot(g.naccu-i+1:g.naccu,:))];
else
Pboot = g.pboot;
end
if isnan(g.rboot)
Rboot = abs(Rboot) / Rbn;
Rboot = sort(Rboot');
Rboot = mean(Rboot(g.naccu-i+1:g.naccu,:));
else
Rboot = g.rboot;
end
else
myprintf(g.verbose,'No valid bootstrap trials...!\n');
end
end
switch lower(g.plotitc)
case 'on',
switch lower(g.plotersp),
case 'on', ordinate1 = 0.67; ordinate2 = 0.1; height = 0.33; g.plot = 1;
case 'off', ordinate2 = 0.1; height = 0.9; g.plot = 1;
end;
case 'off', ordinate1 = 0.1; height = 0.9;
switch lower(g.plotersp),
case 'on', ordinate1 = 0.1; height = 0.9; g.plot = 1;
case 'off', g.plot = 0;
end;
end;
if g.plot
myprintf(g.verbose,'\nNow plotting...\n');
set(gcf,'DefaultAxesFontSize',AXES_FONT)
colormap(jet(256));
pos = get(gca,'position');
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)];
end
switch lower(g.plotersp)
case 'on'
%
%%%%%%% image the ERSP %%%%%%%%%%%%%%%%%%%%%%%%%%
%
h(1) = subplot('Position',[.1 ordinate1 .9 height].*s+q);
PP = P; % PP will be ERSP power after
if ~isnan(g.alpha) % zero out nonsignif. power differences
PP(find((PP > repmat(Pboot(1,:)',[1 g.timesout])) ...
& (PP < repmat(Pboot(2,:)',[1 g.timesout])))) = 0;
end
if ERSP_CAXIS_LIMIT == 0
ersp_caxis = [-1 1]*1.1*max(max(abs(P(dispf,:))));
else
ersp_caxis = ERSP_CAXIS_LIMIT*[-1 1];
end
if ~isnan( g.baseline )
imagesc(times,freqs(dispf),PP(dispf,:),ersp_caxis);
else
imagesc(times,freqs(dispf),PP(dispf,:));
end
set(gca,'ydir',g.hzdir); % make frequency ascend or descend
if ~isempty(g.erspmax)
caxis([-g.erspmax g.erspmax]);
end
hold on
plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth); % plot time 0
if ~isnan(g.marktimes) % plot marked time
for mt = g.marktimes(:)'
plot([mt mt],[0 freqs(max(dispf))],'--k','LineWidth',g.linewidth);
end
end
hold off
set(h(1),'YTickLabel',[],'YTick',[])
set(h(1),'XTickLabel',[],'XTick',[])
if ~isempty(g.vert)
for index = 1:length(g.vert)
line([g.vert(index), g.vert(index)], [min(freqs(dispf)) max(freqs(dispf))], 'linewidth', 1, 'color', 'm');
end
end
h(2) = gca;
h(3) = cbar('vert'); % ERSP colorbar axes
set(h(2),'Position',[.1 ordinate1 .8 height].*s+q)
set(h(3),'Position',[.95 ordinate1 .05 height].*s+q)
title([ 'ERSP(' g.unitpower ')' ])
E = [min(P(dispf,:));max(P(dispf,:))];
h(4) = subplot('Position',[.1 ordinate1-0.1 .8 .1].*s+q); % plot marginal ERSP means
% below the ERSP image
plot(times,E,[0 0],...
[min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3], ...
'--m','LineWidth',g.linewidth)
axis([min(times) max(times) ...
min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3])
tick = get(h(4),'YTick');
set(h(4),'YTick',[tick(1) ; tick(end)])
set(h(4),'YAxisLocation','right')
set(h(4),'TickLength',[0.020 0.025]);
xlabel('Time (ms)')
ylabel( g.unitpower )
E = 10 * log10(mbase(dispf));
h(5) = subplot('Position',[0 ordinate1 .1 height].*s+q); % plot mean spectrum
% to left of ERSP image
plot(freqs(dispf),E,'LineWidth',g.linewidth)
if ~isnan(g.alpha)
hold on;
plot(freqs(dispf),Pboot(:,dispf)+[E;E],'g', 'LineWidth',g.linewidth);
plot(freqs(dispf),Pboot(:,dispf)+[E;E],'k:','LineWidth',g.linewidth)
end
axis([freqs(1) freqs(max(dispf)) min(E)-max(abs(E))/3 max(E)+max(abs(E))/3])
tick = get(h(5),'YTick');
if (length(tick)>1)
set(h(5),'YTick',[tick(1) ; tick(end-1)])
end
set(h(5),'TickLength',[0.020 0.025]);
set(h(5),'View',[90 90])
xlabel('Frequency (Hz)')
ylabel( g.unitpower )
if strcmp(g.hzdir,'normal')
freqdir = 'reverse';
else
freqdir = 'normal';
end
set(h(5),'xdir',freqdir); % make frequency ascend or descend
end
switch lower(g.plotitc)
case 'on'
%
%%%%%%%%%%%% Image the ITC %%%%%%%%%%%%%%%%%%
%
h(6) = subplot('Position',[.1 ordinate2 .9 height].*s+q); % ITC image
RR = R; % RR is the masked ITC (R)
if ~isnan(g.alpha)
RR(find(RR < repmat(Rboot(1,:)',[1 g.timesout]))) = 0;
end
if ITC_CAXIS_LIMIT == 0
coh_caxis = min(max(max(R(dispf,:))),1)*[-1 1]; % 1 WAS 0.4 !
else
coh_caxis = ITC_CAXIS_LIMIT*[-1 1];
end
if exist('Rsign') && strcmp(g.plotphase, 'on')
imagesc(times,freqs(dispf),Rsign(dispf,:).*RR(dispf,:),coh_caxis); % <---
else
imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % <---
end
if ~isempty(g.itcmax)
caxis([-g.itcmax g.itcmax]);
end
tmpcaxis = caxis;
set(gca,'ydir',g.hzdir); % make frequency ascend or descend
hold on
plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
if ~isnan(g.marktimes)
for mt = g.marktimes(:)'
plot([mt mt],[0 freqs(max(dispf))],'--k','LineWidth',g.linewidth);
end
end
hold off
set(h(6),'YTickLabel',[],'YTick',[])
set(h(6),'XTickLabel',[],'XTick',[])
if ~isempty(g.vert)
for index = 1:length(g.vert)
line([g.vert(index), g.vert(index)], ...
[min(freqs(dispf)) max(freqs(dispf))], ...
'linewidth', 1, 'color', 'm');
end
end
h(7) = gca;
h(8) = cbar('vert');
%h(9) = get(h(8),'Children');
set(h(7),'Position',[.1 ordinate2 .8 height].*s+q)
set(h(8),'Position',[.95 ordinate2 .05 height].*s+q)
set(h(8),'YLim',[0 tmpcaxis(2)]);
title('ITC')
%
%%%%% plot the ERP below the ITC image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% E = mean(R(dispf,:));
ERPmax = max(ERP);
ERPmin = min(ERP);
ERPmax = ERPmax + 0.1*(ERPmax-ERPmin);
ERPmin = ERPmin - 0.1*(ERPmax-ERPmin);
h(10) = subplot('Position',[.1 ordinate2-0.1 .8 .1].*s+q); % ERP
plot(ERPtimes,ERP(ERPindices),...
[0 0],[ERPmin ERPmax],'--m','LineWidth',g.linewidth);
hold on; plot([times(1) times(length(times))],[0 0], 'k');
axis([min(ERPtimes) max(ERPtimes) ERPmin ERPmax]);
tick = get(h(10),'YTick');
set(h(10),'YTick',[tick(1) ; tick(end)])
set(h(10),'TickLength',[0.02 0.025]);
set(h(10),'YAxisLocation','right')
xlabel('Time (ms)')
ylabel('\muV')
if (~isempty(g.topovec))
if length(g.topovec) ~= 1, ylabel(''); end; % ICA component
end
E = mean(R(dispf,:)');
h(11) = subplot('Position',[0 ordinate2 .1 height].*s+q); % plot the marginal mean
% ITC left of the ITC image
if ~isnan(g.alpha)
plot(freqs(dispf),E,'LineWidth',g.linewidth); hold on;
plot(freqs(dispf),Rboot(dispf),'g', 'LineWidth',g.linewidth);
plot(freqs(dispf),Rboot(dispf),'k:','LineWidth',g.linewidth);
axis([freqs(1) freqs(max(dispf)) 0 max([E Rboot(dispf)])+max(E)/3])
else
plot(freqs(dispf),E,'LineWidth',g.linewidth)
axis([freqs(1) freqs(max(dispf)) min(E)-max(E)/3 max(E)+max(E)/3])
end
tick = get(h(11),'YTick');
set(h(11),'YTick',[tick(1) ; tick(length(tick))])
set(h(11),'View',[90 90])
set(h(11),'TickLength',[0.020 0.025]);
xlabel('Frequency (Hz)')
ylabel('ERP')
if strcmp(g.hzdir,'normal')
freqdir = 'reverse';
else
freqdir = 'normal';
end
set(gca,'xdir',freqdir); % make frequency ascend or descend
%
%%%%%%%%%%%%%%% plot a topoplot() %%%%%%%%%%%%%%%%%%%%%%%
%
if (~isempty(g.topovec))
h(12) = subplot('Position',[-.1 .43 .2 .14].*s+q);
if length(g.topovec) == 1
topoplot(g.topovec,g.elocs,'electrodes','off', ...
'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
else
topoplot(g.topovec,g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
end
axis('square')
end
end; % switch
if g.plot
try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end
if (length(g.title) > 0)
axes('Position',pos,'Visible','Off');
h(13) = text(-.05,1.01,g.title);
set(h(13),'VerticalAlignment','bottom')
set(h(13),'HorizontalAlignment','left')
set(h(13),'FontSize',TITLE_FONT);
end
axcopy(gcf);
end
% symmetric Hanning tapering function
% -----------------------------------
function w = hanning(n)
if ~rem(n,2)
w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
w = [w; w(end:-1:1)];
else
w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
w = [w; w(end-1:-1:1)];
end
function myprintf(verbose, varargin)
if strcmpi(verbose, 'on')
fprintf(varargin{:});
end
Datasets
Models
