% PAC_CONT - compute phase-amplitude coupling (power of first input
% correlation with phase of second). There is no graphical output
% to this function.
%
% Usage:
% >> pac_cont(x,y,srate);
% >> [pac timesout pvals] = pac_cont(x,y,srate,'key1', 'val1', 'key2', val2' ...);
%
% Inputs:
% x = [float array] 1-D data vector of size (1xtimes)
% y = [float array] 1-D data vector of size (1xtimes)
% srate = data sampling rate (Hz)
%
% Optional time information inputs:
% '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. Note that this
% parameter is overwritten if the minimum frequency has been set
% manually and requires a longer time window {~frames/8}
% 'ntimesout' = Number of output times (int<frames-winframes). Enter a
% negative value [-S] to subsample original time by S.
% 'timesout' = Enter an array to obtain spectral decomposition at
% specific time values (note: algorithm find closest time
% point in data and this might result in an unevenly spaced
% time array). Overwrite 'ntimesout'. {def: automatic}
% 'tlimits' = [min max] time limits in ms.
%
% Optional PAC inputs:
% 'method' = ['modulation'|'plv'|'corr'|'glm'] (see reference).
% 'freqphase' = [min max] frequency limits. Default [minfreq 50],
% minfreq being determined by the number of data points,
% cycles and sampling frequency. Use 0 for minimum frequency
% to compute default minfreq. You may also enter an
% array of frequencies for the spectral decomposition
% (for FFT, closest computed frequency will be returned; use
% 'padratio' to change FFT freq. resolution).
% 'freqamp' = [float array] array of frequencies for the second
% argument. 'freqs' is used for the first argument.
% By default it is the same as 'freqs'.
% 'filterfunc' = ['eegfilt'|'iirfilt'|'eegfftfilt'] filtering function.
% Default is iirfilt. Warning, filtering may dramatically
% affect the result. With the 'corr' method, make sure you
% have a large window size because each window is filtered
% independently.
% 'filterphase' = @f_handle. Function handle to filter the data for the
% phase information. For example, @(x)iirfilt(x, 1000, 2,
% 20). Note that 'freqphase' is ignore in this case.
% 'filteramp' = @f_handle. Function handle to filter the data for the
% amplitude information. Note that 'freqamp' is ignore in
% this case.
%
% Inputs for statistics:
% 'alpha' = [float] p-value threshold. Default is none (no statistics).
% 'mcorrect' = ['none'|'fdr'] method to correct for multiple comparison.
% Default is 'none'.
% 'baseline' = [min max] baseline period for the Null distribution. Default
% is the whole data range. Note that this option is ignored
% for instanstaneous statistics.
% 'instantstat' = ['on'|'off'] performs statistics for each time window
% independently. Default is 'off'.
% 'naccu' = [integer] number of accumulations for surrogate
% statistics.
% 'statlim' = ['parametric'|'surrogate'] use a parametric methods to
% asseess the limit of the surrogate distribution or use
% the tail of the distribution ('surrogate' method)
%
% Other inputs:
% 'title' = [string] figure title. Default is none.
% 'vert' = [float array] array of time value for which to plot
% vertical lines. Default is none.
%
% Outputs:
% pac = Phase-amplitude coupling values.
% timesout = vector of time indices
% pvals = Associated p-values
%
% Author: Arnaud Delorme and Makoto Miyakoshi, SCCN/INC, UCSD 2012-
%
% References:
% Methods used here are introduced and compared in:
% Penny, Duzel, Miller, Ojemann. (20089). Testing for Nested Oscilations.
% J Neuro Methods. 174:50-61
%
% Modulation Index is defined in:
% Canolty, Edwards, Dalal, Soltani, Nagarajan, Kirsch, et al. (2006). Modulation index is defined in High Gamma Power Is Phase-Locked to Theta
% Oscillations in Human Neocortex. Science. 313:1626-8.
%
% PLV (Phase locking value) is defined in:
% Lachaux, Rodriguez, Martiniere, Varela. (1999). Measuring phase synchrony
% in brain signal. Hum Brain Mapp. 8:194-208.
%
% corr (correlation) method is defined in:
% Brunce, Eckhorn. (2004). Task-related coupling from high- to
% low-frequency signals among visual cortical areas in human subdural
% recordings. Int J Psychophysiol. 51:97-116.
%
% glm (general linear model) is defined in
% Penny, Duzel, Miller, Ojemann. (20089). Testing for Nested Oscilations.
% J Neuro Methods. 174:50-61
% Copyright (C) 2012 Arnaud Delorme, UCSD
%
% 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 [m_raw pvals indexout] = pac_cont(X, Y, srate, varargin);
if nargin < 1
help pac_cont;
return;
end
% deal with 3-D inputs
% --------------------
if ndims(X) == 3 || ndims(Y) == 3, error('Cannot process 3-D input'); end
if size(X,1) > 1, X = X'; end
if size(Y,1) > 1, Y = Y'; end
if size(X,1) ~= 1 || size(Y,1) ~= 1, error('Cannot only process vector input'); end
frame = size(X,2);
pvals = [];
g = finputcheck(varargin, ...
{ ...
'alpha' 'real' [0 0.2] [];
'baseline' 'float' [] [];
'freqphase' 'real' [0 Inf] [0 srate/2];
'freqamp' 'real' [0 Inf] [];
'mcorrect' 'string' { 'none' 'fdr' } 'none';
'method' 'string' { 'plv' 'modulation' 'glm' 'corr' } 'modulation';
'naccu' 'integer' [1 Inf] 250;
'instantstat' 'string' {'on','off'} 'off';
'newfig' 'string' {'on','off'} 'on';
'nofig' 'string' {'on','off'} 'off';
'statlim' 'string' {'surrogate','parametric'} 'parametric';
'timesout' 'real' [] []; ...
'filterfunc' 'string' { 'eegfilt' 'iirfilt' 'eegfiltfft' } 'eegfiltfft'; ...
'filterphase' '' {} [];
'filteramp' '' {} [];
'ntimesout' 'integer' [] 200; ...
'tlimits' 'real' [] [0 frame/srate];
'title' 'string' [] '';
'vert' 'real' [] [];
'winsize' 'integer' [0 Inf] max(pow2(nextpow2(frame)-3),4) }, 'pac');
if ischar(g), error(g); end
if ~isempty(g.filterphase)
x_freqphase = feval(g.filterphase, X(:)');
else x_freqphase = feval(g.filterfunc, X(:)', srate, g.freqphase(1), g.freqphase(end));
end
if ~isempty(g.filteramp)
x_freqamp = feval(g.filteramp, Y(:)');
else x_freqamp = feval(g.filterfunc, Y(:)', srate, g.freqamp( 1), g.freqamp( end));
end
z_phasedata = hilbert(x_freqphase);
z_ampdata = hilbert(x_freqamp);
phase = angle(z_phasedata);
amplitude = abs( z_ampdata);
z = amplitude.*exp(i*phase); % this is the pac measure
% get time windows
% ----------------
g.verbose = 'on';
g.causal = 'off';
[ timesout1 indexout ] = gettimes(frame, g.tlimits, g.timesout, g.winsize, g.ntimesout, g.causal, g.verbose);
% scan time windows
% -----------------
if ~isempty(g.alpha)
m_raw = zeros(1,length(indexout));
pvals = zeros(1,length(indexout));
end
fprintf('Computing PAC:\n');
for iWin = 1:length(indexout)
x_phaseEpoch = x_freqphase(indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
x_ampEpoch = x_freqamp( indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
z_phaseEpoch = z_phasedata(indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
z_ampEpoch = z_ampdata( indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
z_epoch = z( indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
numpoints=length(x_phaseEpoch);
if rem(iWin,10) == 0, verboseprintf(g.verbose, ' %d',iWin); end
if rem(iWin,120) == 0, verboseprintf(g.verbose, '\n'); end
% Choose method
% -------------
if strcmpi(g.method, 'modulation')
% Modulation index
m_raw(iWin) = abs(sum(z_epoch))/numpoints;
elseif strcmpi(g.method, 'plv')
if iWin == 145
%dsfsd;
end
%amplitude_filt = sgolayfilt(amplitude, 3, 101);
if ~isempty(g.filterphase)
amplitude_filt = feval(g.filterphase, z_ampEpoch);
else amplitude_filt = feval(g.filterfunc , z_ampEpoch, srate, g.freqphase(1), g.freqphase(end));
end
z_amplitude_filt = hilbert(amplitude_filt);
phase_amp_modulation = angle(z_amplitude_filt);
m_raw(iWin) = abs(sum(exp(i*(x_phaseEpoch - phase_amp_modulation)))/numpoints);
elseif strcmpi(g.method, 'corr')
if iWin == inf %145
figure; plot(abs(z_ampdata))
hold on; plot(x_phasedata/10000000000, 'r')
x = X(indexout(iWin)+[-g.winsize/2+1:g.winsize/2]);
hold on; plot(x, 'g');
dsfsd;
end
[r_ESC pval_corr] = corrcoef(x_phaseEpoch, abs(z_ampEpoch));
m_raw(iWin) = r_ESC(1,2);
pvals(iWin) = pval_corr(1,2);
elseif strcmpi(g.method, 'glm')
[b dev stats] = glmfit(x_phaseEpoch', abs(z_ampEpoch)', 'normal');
GLM_beta = stats.beta(2,1);
pvals(iWin) = stats.p(2,1);
m_raw(iWin) = b(1);
end
%% compute statistics (instantaneous)
% -----------------------------------
if ~isempty(g.alpha) && strcmpi(g.instantstat, 'on') && ~strcmpi(g.method, 'corr') && ~strcmpi(g.method, 'glm')
% compute surrogate values
numsurrogate=g.naccu;
minskip=srate;
maxskip=numpoints-srate; % max variation half a second
if maxskip < 1
error('Window size shorter than 1 second; too short for computing surrogate data');
end
skip=ceil(numpoints.*rand(numsurrogate*4,1));
skip(skip>maxskip)=[];
skip(skip<minskip)=[];
skip=skip(1:numsurrogate,1);
surrogate_m=zeros(numsurrogate,1);
for s=1:numsurrogate
surrogate_amplitude=[amplitude(skip(s):end) amplitude(1:skip(s)-1)]; % consider circular shifts
surrogate_m(s)=abs(mean(surrogate_amplitude.*exp(i*phase)));
%disp(numsurrogate-s)
end
if strcmpi(g.statlim, 'surrogate')
pvals(iWin) = stat_surrogate_pvals(surrogate_m, m_raw(iWin), 'upper');
%fprintf('Raw PAC is %3.2f (p-value=%1.3f)\n', m_raw(iWin), pvals(iWin));
else
% Canolty method below
%% fit gaussian to surrogate data, uses normfit.m from MATLAB Statistics toolbox
[surrogate_mean,surrogate_std]=normfit(surrogate_m);
%% normalize length using surrogate data (z-score)
m_norm_length=(abs(m_raw(iWin))-surrogate_mean)/surrogate_std;
pvals(iWin) = normcdf(0, m_norm_length, 1);
m_norm_phase=angle(m_raw(iWin));
m_norm=m_norm_length*exp(i*m_norm_phase);
% compare parametric and non-parametric methods (return similar
% results)
if iWin == length(indexout)
figure;
plot(-log10(pvals));
hold on; plot(-log10(pvals2), 'r');
end
end
end
end
fprintf('\n');
% Computes alpha
% --------------
if ~isempty(g.alpha) && strcmpi(g.instantstat, 'off')
if isempty(g.baseline)
g.baseline = [ timesout1(1) timesout1(end) ];
end
baselineInd = find(timesout1 >= g.baseline(1) & timesout1 <= g.baseline(end));
m_raw_base = abs(m_raw(baselineInd));
if strcmpi(g.statlim, 'surrogate')
for index = 1:length(m_raw)
pvals(index) = stat_surrogate_pvals(m_raw_base, m_raw(index), 'upper');
end
else
[surrogate_mean,surrogate_std]=normfit(m_raw_base);
m_norm_length=(abs(m_raw)-surrogate_mean)/surrogate_std;
pvals = normcdf(0, m_norm_length, 1);
end
if strcmpi(g.mcorrect, 'fdr')
pvals = fdr(pvals);
end
end
%% plot results
% -------------
if strcmpi(g.nofig, 'on')
return
end
if strcmpi(g.newfig, 'on')
figure;
end
if ~isempty(g.alpha)
plotcurve(timesout1, m_raw, 'maskarray', pvals < g.alpha);
else plotcurve(timesout1, m_raw);
end
xlabel('Time (ms)');
ylabel('PAC (0 to 1)');
title(g.title);
% plot vertical lines
% -------------------
if ~isempty(g.vert)
hold on;
yl = ylim;
for index = 1:length(g.vert)
plot([g.vert(index) g.vert(index)], yl, 'g');
end
end
% -------------
% gettime function identical to timefreq function
% DO NOT MODIFY
% -------------
function [ timevals, timeindices ] = gettimes(frames, tlimits, timevar, winsize, ntimevar, causal, verbose);
timevect = linspace(tlimits(1), tlimits(2), frames);
srate = 1000*(frames-1)/(tlimits(2)-tlimits(1));
if isempty(timevar) % no pre-defined time points
if ntimevar(1) > 0
% generate linearly space vector
% ------------------------------
if (ntimevar > frames-winsize)
ntimevar = frames-winsize;
if ntimevar < 0
error('Not enough data points, reduce the window size or lowest frequency');
end
verboseprintf(verbose, ['Value of ''timesout'' must be <= frame-winsize, ''timesout'' adjusted to ' int2str(ntimevar) '\n']);
end
npoints = ntimevar(1);
wintime = 500*winsize/srate;
if strcmpi(causal, 'on')
timevals = linspace(tlimits(1)+2*wintime, tlimits(2), npoints);
else timevals = linspace(tlimits(1)+wintime, tlimits(2)-wintime, npoints);
end
verboseprintf(verbose, 'Generating %d time points (%1.1f to %1.1f ms)\n', npoints, min(timevals), max(timevals));
else
% subsample data
% --------------
nsub = -ntimevar(1);
if strcmpi(causal, 'on')
timeindices = [ceil(winsize+nsub):nsub:length(timevect)];
else timeindices = [ceil(winsize/2+nsub/2):nsub:length(timevect)-ceil(winsize/2)-1];
end
timevals = timevect( timeindices ); % the conversion at line 741 leaves timeindices unchanged
verboseprintf(verbose, 'Subsampling by %d (%1.1f to %1.1f ms)\n', nsub, min(timevals), max(timevals));
end
else
timevals = timevar;
% check boundaries
% ----------------
wintime = 500*winsize/srate;
if strcmpi(causal, 'on')
tmpind = find( (timevals >= tlimits(1)+2*wintime-0.0001) & (timevals <= tlimits(2)) );
else tmpind = find( (timevals >= tlimits(1)+wintime-0.0001) & (timevals <= tlimits(2)-wintime+0.0001) );
end
% 0.0001 account for numerical inaccuracies on opteron computers
if isempty(tmpind)
error('No time points. Reduce time window or minimum frequency.');
end
if length(timevals) ~= length(tmpind)
verboseprintf(verbose, 'Warning: %d out of %d time values were removed (now %3.2f to %3.2f ms) so the lowest\n', ...
length(timevals)-length(tmpind), length(timevals), timevals(tmpind(1)), timevals(tmpind(end)));
verboseprintf(verbose, ' frequency could be computed with the requested accuracy\n');
end
timevals = timevals(tmpind);
end
% find closet points in data
% --------------------------
timeindices = round(eeg_lat2point(timevals, 1, srate, tlimits, 1E-3));
if length(timeindices) < length(unique(timeindices))
timeindices = unique_bc(timeindices)
verboseprintf(verbose, 'Warning: duplicate times, reduce the number of output times\n');
end
if length(unique(timeindices(2:end)-timeindices(1:end-1))) > 1
verboseprintf(verbose, 'Finding closest points for time variable\n');
verboseprintf(verbose, 'Time values for time/freq decomposition is not perfectly uniformly distributed\n');
else
verboseprintf(verbose, 'Distribution of data point for time/freq decomposition is perfectly uniform\n');
end
timevals = timevect(timeindices);
function verboseprintf(verbose, varargin)
if strcmpi(verbose, 'on')
fprintf(varargin{:});
end
Datasets
Models
