% STD_PAC - Compute or read PAC data (Phase Amplitude Coupling).
%
% Usage:
% >> [X times logfreqs ] = std_pac(EEG, 'key', 'val', ...);
% Inputs:
% EEG - an EEG dataset structure.
%
% Optional inputs:
% 'components1'- [numeric vector] components in the EEG structure used
% for spectral amplitude in PAC {default|[]: all }
% 'components2'- [numeric vector] components in the EEG structure used
% for phase in PAC {default|[]: all }
% 'channels1' - [numeric vector or cell array of channel labels] channels
% in the EEG structure for spectral amplitude in PAC
% {default|[]: no channels}
% 'channels2' - [numeric vector or cell array of channel labels] channels
% in the EEG structure for phase in PAC
% {default|[]: no channels}
% 'freqs' - [minHz maxHz] the PAC frequency range to compute power.
% {default: 12 to EEG sampling rate divided by 2}
% 'cycles' - [wavecycles (factor)]. If 0 -> DFT (constant window length
% across frequencies).
% If >0 -> the number of cycles in each analysis wavelet.
% If [wavecycles factor], wavelet cycles increase with
% frequency, beginning at wavecyles. (0 < factor < 1)
% factor = 0 -> fixed epoch length (DFT, as in FFT).
% factor = 1 -> no increase (standard wavelets)
% {default: [0]}
% 'freqphase' - [valHz] single number for computing the phase at a given
% frequency.
% 'cyclephase' - [valcycle] single cycle number.
% 'timewindow' - [minms maxms] time window (in ms) to plot.
% {default: all output latencies}
% 'padratio' - (power of 2). Multiply the number of output frequencies
% by dividing their frequency spacing through 0-padding.
% Output frequency spacing is (low_freq/padratio).
% 'recompute' - ['on'|'off'] 'on' forces recomputation of PAC.
% {default: 'off'}
%
% Other optional inputs:
% This function will take any of the NEWTIMEF optional inputs (for instance
% to compute log-space frequencies)...
%
% Outputs:
% X - the PAC of the requested ICA components/channels
% in the selected frequency and time range.
% times - vector of time points for which the PAC were computed.
% freqs - vector of frequencies (in Hz) at which the
% PAC was evaluated.
%
% Files written or modified:
% [dataset_filename].icapac <-- saved component PAC
% OR for channels
% [dataset_filename].datpac <-- saved channel PAC
%
% See also: TIMEF, STD_ITC, STD_ERP, STD_SPEC, STD_TOPO, STD_PRECLUST
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, July, 2009-
% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, arno@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.
function [X, times, freqs, parameters] = std_pac(EEG, varargin)
if nargin < 1
help std_pac;
return;
end
options = {};
[g timefargs] = finputcheck(varargin, { ...
'components1' 'integer' [] [];
'channels1' { 'cell','integer' } { [],[] } {};
'components2' 'integer' [] [];
'channels2' { 'cell','integer' } { [],[] } {};
'outputfile' 'string' [] '';
'powbase' 'real' [] [];
'plot' 'string' { 'on','off' } 'off';
'recompute' 'string' { 'on','off' } 'off';
'getparams' 'string' { 'on','off' } 'off';
'timerange' 'real' [] [];
'freqrange' 'real' [] [];
'padratio' 'real' [] 1;
'freqs' 'real' [] [12 EEG.srate/2];
'cycles' 'real' [] [8];
'freqphase' 'real' [] [5];
'cyclephase' 'real' [] [3];
'interp' 'struct' { } struct([]);
'rmcomps' 'integer' [] [];
'freqscale' 'string' [] 'log' }, 'std_pac', 'ignore');
if ischar(g), error(g); end
% checking input parameters
% -------------------------
if isempty(g.components1) && isempty(g.channels1)
if isempty(EEG(1).icaweights)
error('EEG.icaweights not found');
end
g.components1 = 1:size(EEG(1).icaweights,1);
g.components2 = 1:size(EEG(1).icaweights,1);
disp('Computing PAC with default values for all components of the dataset');
end
% select ICA components or data channels
% --------------------------------------
if ~isempty(g.outputfile)
filenamepac = fullfile('', [ g.outputfile '.datpac' ]);
g.indices1 = std_chaninds(EEG, g.channels1);
g.indices2 = std_chaninds(EEG, g.channels2);
prefix = 'chan';
elseif ~isempty(g.components1)
g.indices1 = g.components1;
g.indices2 = g.components2;
prefix = 'comp';
filenamepac = fullfile(EEG.filepath, [ EEG.filename(1:end-3) 'icapac' ]);
if ~isempty(g.channels1)
error('Cannot compute PAC for components and channels at the same time');
end
elseif ~isempty(g.channels1)
g.indices1 = std_chaninds(EEG, g.channels1);
g.indices2 = std_chaninds(EEG, g.channels2);
prefix = 'chan';
filenamepac = fullfile(EEG.filepath, [ EEG.filename(1:end-3) 'datpac' ]);
end
% Compute PAC parameters
% -----------------------
parameters = { 'wavelet', g.cycles, 'padratio', g.padratio, ...
'freqs2', g.freqphase, 'wavelet2', g.cyclephase, 'freqscale', g.freqscale, timefargs{:} };
if length(g.freqs)>0, parameters = { parameters{:} 'freqs' g.freqs }; end
% Check if PAC information found in datasets and if fits requested parameters
% ----------------------------------------------------------------------------
if exist( filenamepac ) && strcmpi(g.recompute, 'off')
fprintf('Use existing file for PAC: %s\n', filenamepac);
if ~isempty(g.components1)
[X, times, freqs, parameters] = std_readpac(EEG, 1, g.indices1, g.indices2, g.timerange, g.freqrange);
else [X, times, freqs, parameters] = std_readpac(EEG, 1, -g.indices1, -g.indices2, g.timerange, g.freqrange);
end
return;
end
% return parameters
% -----------------
if strcmpi(g.getparams, 'on')
X = []; times = []; freqs = [];
return;
end
options = {};
if ~isempty(g.components1)
tmpdata = eeg_getdatact(EEG, 'component', [1:size(EEG(1).icaweights,1)]);
else
EEG.data = eeg_getdatact(EEG, 'channel', [1:EEG.nbchan], 'rmcomps', g.rmcomps);
if ~isempty(g.rmcomps), options = { options{:} 'rmcomps' g.rmcomps }; end
if ~isempty(g.interp),
EEG = eeg_interp(EEG, g.interp, 'spherical');
options = { options{:} 'interp' g.interp };
end
tmpdata = EEG.data;
end;
% Compute PAC
% -----------
all_pac = [];
for k = 1:length(g.indices1) % for each (specified) component
for l = 1:length(g.indices2) % for each (specified) component
tmpparams = parameters;
% Run pac() to get PAC
% --------------------
timefdata1 = tmpdata(g.indices1(k),:,:);
timefdata2 = tmpdata(g.indices2(l),:,:);
if strcmpi(g.plot, 'on'), figure; end
%[logersp,logitc,logbase,times,logfreqs,logeboot,logiboot,alltfX] ...
[pacvals, times, freqs1, freqs2] = pac( timefdata1, timefdata2, EEG(1).srate, 'tlimits', [EEG.xmin EEG.xmax]*1000, tmpparams{1:end});
all_pac = setfield( all_pac, [ prefix int2str(g.indices1(k)) '_' int2str(g.indices2(l)) '_pac' ], squeeze(single(pacvals )));
end
end
% Save PAC into file
% ------------------
all_pac.freqs = freqs1;
all_pac.times = times;
all_pac.datatype = 'PAC';
all_pac.parameters = tmpparams;
if ~isempty(g.channels1)
if ~isempty(EEG(1).chanlocs)
tmpchanlocs = EEG(1).chanlocs;
all_pac.chanlabels1 = { tmpchanlocs(g.indices1).labels };
all_pac.chanlabels2 = { tmpchanlocs(g.indices2).labels };
end
end
std_savedat( filenamepac , all_pac );
if ~isempty(g.components1)
[X, times, freqs, parameters] = std_readpac(EEG, 1, g.indices1, g.indices2, g.timerange, g.freqrange);
else [X, times, freqs, parameters] = std_readpac(EEG, 1, -g.indices1, -g.indices2, g.timerange, g.freqrange);
end
% --------------------------------------------------------
% -------------------- READ PAC DATA ---------------------
% --------------------------------------------------------
function [pacvals, freqs, timevals, params] = std_readpac(ALLEEG, abset, comp1, comp2, timewindow, freqrange);
if nargin < 5
timewindow = [];
end
if nargin < 6
freqrange = [];
end
% multiple entry
% --------------
if length(comp1) > 1 || length(comp2) > 1
for index1 = 1:length(comp1)
for index2 = 1:length(comp2)
[tmppac, freqs, timevals, params] = std_readpac(ALLEEG, abset, comp1(index1), comp2(index2), timewindow, freqrange);
pacvals(index1,index2,:,:,:) = tmppac;
end
end
return;
end
for k = 1: length(abset)
if comp1 < 0
filename = fullfile( ALLEEG(abset(k)).filepath,[ ALLEEG(abset(k)).filename(1:end-3) 'datpac']);
comp1 = -comp1;
comp2 = -comp2;
prefix = 'chan';
else
filename = fullfile( ALLEEG(abset(k)).filepath,[ ALLEEG(abset(k)).filename(1:end-3) 'icapac']);
prefix = 'comp';
end
try
tmppac = load( '-mat', filename, 'parameters', 'times', 'freqs');
catch
error( [ 'Cannot read file ''' filename '''' ]);
end
tmppac.parameters = removedup(tmppac.parameters);
params = struct(tmppac.parameters{:});
params.times = tmppac.times;
params.freqs = tmppac.freqs;
if isempty(comp1)
pacvals = [];
freqs = [];
timevals = [];
return;
end
tmppac = load( '-mat', filename, 'parameters', 'times', 'freqs', ...
[ prefix int2str(comp1) '_' int2str(comp2) '_pac']);
pacall{k} = double(getfield(tmppac, [ prefix int2str(comp1) '_' int2str(comp2) '_pac']));
tlen = length(tmppac.times);
flen = length(tmppac.freqs);
end
% select plotting or clustering time/freq range
% ---------------------------------------------
if ~isempty(timewindow)
if timewindow(1) > tmppac.times(1) || timewindow(end) < tmppac.times(end)
maxind = max(find(tmppac.times <= timewindow(end)));
minind = min(find(tmppac.times >= timewindow(1)));
else
minind = 1;
maxind = tlen;
end
else
minind = 1;
maxind = tlen;
end
if ~isempty(freqrange)
if freqrange(1) > exp(1)^tmppac.freqs(1) || freqrange(end) < exp(1)^tmppac.freqs(end)
fmaxind = max(find(tmppac.freqs <= freqrange(end)));
fminind = min(find(tmppac.freqs >= freqrange(1)));
else
fminind = 1;
fmaxind = flen;
end
else
fminind = 1;
fmaxind = flen;
end
% return parameters
% ----------------
for cond = 1:length(abset)
try
pac = pacall{cond}(fminind:fmaxind,minind:maxind);
catch
pac = pacall{cond}; % for 'method', 'latphase'
end
pacvals(:,:,cond) = pac;
end
freqs = tmppac.freqs(fminind:fmaxind);
timevals = tmppac.times(minind:maxind);
% remove duplicates in the list of parameters
% -------------------------------------------
function cella = removedup(cella)
[tmp indices] = unique_bc(cella(1:2:end));
if length(tmp) ~= length(cella)/2
%fprintf('Warning: duplicate ''key'', ''val'' parameter(s), keeping the last one(s)\n');
end
cella = cella(sort(union(indices*2-1, indices*2)));
Datasets
Models
