% CORRIMAGE - compute correlation image between an event and amplitude
% and phase in the time-frequency domain.
%
% Usage:
% corrimage(data, sortvar, times);
% [times, freqs, alpha, sigout, limits ] = corrimage(data, sortvar, ...
% times, 'key1', val1, 'key2', val2, ...)
%
% Inputs:
% data - [float array] data of size (times,trials)
% sortvar - [float array] sorting variable of size (trials)
% times - [float array] time indices for every data point. Size of
% (times). Note: same input as the times vector for erpimage.
%
% Optional input:
% 'mode' - ['phase'|'amp'] compute correlation of event values with
% phase ('phase') or amplitude ('amp') of signal at the
% given time-frequency points. Default is 'amp'.
% 'freqs' - [min nfreqs max] build a frequency vector of size nfreqs
% with frequency spaced in a log scale. Then compute
% correlation at these frequency value.
% Default is 50 points in a log scale between 2.5Hz to 50Hz.
% 'times' - [float vector] compute correlation at these time value (ms).
% (uses closest times points in 'times' input and return
% them in the times output).
% Default is 100 steps between min time+5% of time range and
% max time-5%. Enter only 2 values [N X] to generate N
% time points and trim by X %. If N is negative, uses it as
% a subsampling factor [-3 5] trims times by 5% and subsample
% by 3 (by subsampling one obtains a regularly spaced times).
% 'trim' - [low high] low and high percentile of sorted sortvar values
% to retain. i.e. [5 95] remove the 5 lowest and highest
% percentile of sortvar values (and associated data) before
% computing statistics. Default is [0 100].
% 'align' - [float] same as 'align' parameter of ERPIMAGE. This
% parameter is used to constrain the 'times' parameter so
% correlation with data trials containing 0-values (as a
% result of data alignment) are avoided: computing these
% correlations would produce spurious significant results.
% Default is no alignment.
% 'method' - ['erpimage'|timefreq'] use either the ERPIMAGE function
% of the TIMEFREQ function to compute spectral decomposition.
% Default is 'timefreq' for speed reasons (note that both
% methods should return the same results).
% 'erpout' - [min max] regress out ERP using the selected time-window [min
% max] in ms (the ERP is subtracted from the whole time period
% but only regressed out in the selected time window).
% 'triallimit' - [integer array] specify trial boundaries for subjects. For
% instance [1 200 400] indicates 2 subjects, trials 1 to 199
% for subject 1 and trials 200 to 399 for subject 2. This is
% currently only used for regressing erp out.
%
% Processing options:
% 'erpimopt' - [cell array] erpimage additional options (number of cycle ...).
% 'tfopt' - [cell array] timefreq additional options (number of cycle ...).
%
% Plotting options:
% 'plotvals' - [cell array of output values] enter output values here
% { times freqs alpha sigout} to replot them.
% 'nofig' - ['on'|'off'] do not create figure.
% 'cbar' - ['on'|'off'] plot color bar. Default: 'on'.
% 'smooth' - ['on'|'off'] smooth significance array. Default: 'on'.
% 'plot' - ['no'|'alpha'|'sigout'|'sigoutm'|'sigoutm2'] plot -10*log(alpha)
% values ('alpha'); output signal (slope or ITC) ('sigout'),
% output signal masked by significance ('sigoutm') or the last
% 2 option combined ('sigoutm2'). 'no' prevent the function
% from plotting. In addition, see pmask. Default is 'sigoutmasked'.
% 'pmask' - [real] maximum p value to show in plot. Default is 0.00001
% (0.001 taking into account multiple comparisons (100)). Enter
% 0.9XXX to get the higher tail or a negative value (e.e., -0.001
% to get both tails).
% 'vert' - [real array] time vector for vertivcal lines.
% 'limits' - [min max] plotting limits.
%
% Outputs:
% times - [float vector] vector of times (ms)
% freqs - [float vector] vector of frequencies (Hz)
% alpha - [float array] array (freqs,times) of p-values.
% sigout - [float array] array (freqs,times) of signal out (coherence
% for phase and slope for amplitude).
%
% Important note: the 'timefreq' method may truncate the time window to
% compute the spectral decomposition at the lowest freq.
%
% Author: Arnaud Delorme & Scott Makeig, SCCN UCSD,
% and CNL Salk Institute, 18 April 2003
%
% See also: ERPIMAGE
%123456789012345678901234567890123456789012345678901234567890123456789012
% Copyright (C) 2002 Arnaud Delorme
%
% 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 [time, freq, alpha, sigout, limits, tf, sortvar] = corrimage(data, sortvar, timevect, varargin)
if nargin < 3
help corrimage;
return;
end
if isempty(timevect), timevect = [1 2]; end
% check inputs
% ------------
g = finputcheck(varargin, { ...
'freqs' 'real' [0 Inf] [2.5 50 50];
'times', 'real', [], [100 5];
'mode', 'string', { 'phase', 'amp' }, 'amp';
'vert' 'real' [] [];
'align' { 'real','cell' } [] [];
'plotvals' 'cell' [] {};
'pmask' 'real' [] 0.00001;
'triallimit' 'integer' [] [];
'trim' 'real' [0 100] [0 100];
'limits' 'real' [] [];
'method' 'string' { 'erpimage','timefreq' } 'timefreq';
'plot' 'string' { 'no','alpha','sigout','sigoutm','sigoutp','sigoutm2' } 'sigoutm';
'nofig' 'string' { 'on','off' } 'off';
'cbar' 'string' { 'on','off' } 'on';
'smooth' 'string' { 'on','off' } 'off';
'erpout' 'real' [] [];
'tfopt' 'cell' [] {};
'erpimopt' 'cell' [] {} });
if ischar(g), error(g); end
fprintf('Generating %d frequencies in log scale (ignore message on linear scale)\n', g.freqs(2));
g.freqs = logscale(g.freqs(1), g.freqs(3), g.freqs(2));
frames = length(timevect);
if size(data,1) == 1
data = reshape(data, frames, size(data,2)*size(data,3)/frames);
end
% trim sortvar values
% -------------------
[sortvar sortorder] = sort(sortvar);
data = data(:, sortorder);
len = length(sortvar);
lowindex = round(len*g.trim(1)/100)+1;
highindex = round(len*g.trim(2)/100);
sortvar = sortvar(lowindex:highindex);
data = data(:, lowindex:highindex);
if lowindex ~=1 || highindex ~= length(sortorder)
fprintf('Actual percentiles %1.2f-%1.2f (indices 1-%d -> %d-%d): event vals min %3.2f; max %3.2f\n', ...
100*(lowindex-1)/len, 100*highindex/len, len, lowindex, highindex, min(sortvar), max(sortvar));
end
% assign subject number for each trial
% ------------------------------------
if ~isempty(g.triallimit)
alltrials = zeros(1,len);
for index = 1:length(g.triallimit)-1
alltrials([g.triallimit(index):g.triallimit(index+1)-1]) = index;
end
alltrials = alltrials(sortorder);
alltrials = alltrials(lowindex:highindex);
end
%figure; hist(sortvar)
% constraining time values depending on data alignment
% ----------------------------------------------------
srate = 1000*(length(timevect)-1)/(timevect(end)-timevect(1));
if ~isempty(g.align)
if iscell(g.align)
fprintf('Realigned sortvar plotted at %g ms.\n',g.align{1});
shifts = round((g.align{1}-g.align{2})*srate/1000); % shifts can be positive or negative
g.align = g.align{1};
else
if isinf(g.align)
g.align = median(sortvar);
end
fprintf('Realigned sortvar plotted at %g ms.\n',g.align);
% compute shifts
% --------------
shifts = round((g.align-sortvar)*srate/1000); % shifts can be positive or negative
end
%figure; hist(shifts)
minshift = min(shifts); % negative
maxshift = max(shifts); % positive
if minshift > 0, error('minshift has to be negative'); end
if maxshift < 0, error('maxshift has to be positive'); end
% realign data for all trials
% ---------------------------
aligndata=zeros(size(data))+NaN; % begin with matrix of zeros()
for t=1:size(data,2), %%%%%%%%% foreach trial %%%%%%%%%
shft = shifts(t);
if shft>0, % shift right
aligndata(shft+1:frames,t)=data(1:frames-shft,t);
elseif shft < 0 % shift left
aligndata(1:frames+shft,t)=data(1-shft:frames,t);
else % shft == 0
aligndata(:,t) = data(:,t);
end
end % end trial
data = aligndata(maxshift+1:frames+minshift,:);
if any(any(isnan(data))), error('NaNs remaining after data alignment'); end
timevect = timevect(maxshift+1:frames+minshift);
% take the time vector subset
% ---------------------------
if isempty(timevect), error('Shift too big, empty time vector');
else fprintf('Time vector truncated for data alignment between %1.0f and %1.0f\n', ...
min(timevect), max(timevect));
end;
end
% regress out the ERP from the data (4 times to have residual ERP close to 0)
% ---------------------------------------------------------------------------
if ~isempty(g.erpout)
%data = erpregout(data);
%data = erpregout(data);
%data = erpregout(data);
disp('Regressing out ERP');
erpbeg = mean(data,2);
if ~isempty(g.triallimit)
for index = 1:length(g.triallimit)-1
trials = find(alltrials == index);
%[data(:,trials) erp factors]= erpregout(data(:,trials), [timevect(1) timevect(end)], [300 400]);
erpstart = mean(data(:,trials),2);
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
erpend = mean(data(:,trials),2);
fprintf([ '***********************************************\n' ...
'Ratio in ERP after regression (compare to before) is %3.2f\n' ...
'***********************************************\n' ], mean(erpend./erpstart));
end
end
erpend = mean(data,2);
fprintf([ '***********************************************\n' ...
'Ratio in grand ERP after regression (compare to before) is %3.2f\n' ...
'***********************************************\n' ], mean(erpend./erpbeg));
if 0
%trials = find(alltrials == 1);
data2 = data;
dasf
[data(:,trials) erp factors]= erpregout(data(:,trials), [timevect(1) timevect(end)], [300 400]);
figure; subplot(1,2,1); erpimage(data, sortvar, timevect, '', 300, 500, 'erp');
figure;
subplot(1,2,1); erpimage(data2, sortvar, timevect, '', 300, 500, 'erp');
subplot(1,2,2); erpimage(data , sortvar, timevect, '', 300, 500, 'erp');
figure;
subplot(1,2,1); erpimage(data2(:,trials), sortvar(trials), timevect, '', 0, 0, 'erp');
subplot(1,2,2); erpimage(data (:,trials), sortvar(trials), timevect, '', 0, 0, 'erp');
figure;
for index = 1:length(g.triallimit)-1
subplot(3,5,index);
trials = find(alltrials == index);
erpimage(data(:,trials), sortvar(trials), timevect, '', 50, 100, 'erp');
end
disp('Regressing out ERP');
if ~isempty(g.triallimit)
for index = 1:length(g.triallimit)-1
trials = find(alltrials == index);
[data(:,trials) erp factors]= erpregout(data(:,trials), [timevect(1) timevect(end)], [300 400]);
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
data(:,trials) = erpregout(data(:,trials));
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
%data(:,trials) = erpregout(data(:,trials), [timevect(1) timevect(end)], g.erpout);
end
else
%data = erpregout(data, [timevect(1) timevect(end)], g.erpout);
end
end
end
% generate time points
% --------------------
g.times = gettimes(timevect, g.times);
data = reshape(data, 1, size(data,2)*size(data,1));
% time frequency decomposition
% ----------------------------
if strcmpi(g.method, 'timefreq') && isempty(g.plotvals)
data = reshape(data, length(data)/length(sortvar), length(sortvar));
[tf, g.freqs, g.times] = timefreq(data, srate, 'freqs', g.freqs, 'timesout', g.times, ...
'tlimits', [min(timevect) max(timevect)], 'wavelet', 3);
outvar = sortvar;
end
% compute correlation
% -------------------
if strcmpi(g.mode, 'phase') && isempty(g.plotvals)
for freq = 1:length(g.freqs)
fprintf('Computing correlation with phase %3.2f Hz ----------------------\n', g.freqs(freq));
for time = 1:length(g.times)
if strcmpi(g.method, 'erpimage')
[outdata,outvar,outtrials,limits,axhndls,erp, ...
amps,cohers,cohsig,ampsig,outamps,phsangls,phsamp,sortidx] =erpimage(data,sortvar,timevect, ...
'', 0,0,g.erpimopt{:}, 'phasesort', [g.times(time) 0 g.freqs(freq)], 'noshow', 'yes');
else
phsangls = angle( squeeze(tf(freq, time, :))' );
end
% computing ITCs
[ ITC(freq, time) alpha(freq, time) ] = ...
bootcircle(outvar/mean(outvar), exp(j*phsangls), 'naccu', 250);
% accumulate 200 values, fitted with a normal distribution
% --------------------------------------------------------
%cmplx = outvar .* exp(j*phsangls)/mean(outvar);
%ITC(freq, time) = mean(cmplx);
%alpha(freq,time) = bootstat(outvar/mean(outvar), exp(j*phsangls), 'res = mean(arg1 .* arg2)', ...
% 'naccu', 100, 'vals', abs(ITC(freq, time)), 'distfit', 'norm');
end
end
sigout = ITC;
elseif isempty(g.plotvals)
for freq = 1:length(g.freqs)
fprintf('Computing correlation with amplitude %3.2f Hz ----------------------\n', g.freqs(freq));
for time = 1:length(g.times)
if strcmpi(g.method, 'erpimage')
[outdata,outvar,outtrials,limits,axhndls,erp, ...
amps,cohers,cohsig,ampsig,outamps,phsangls,phsamp,sortidx] =erpimage(data,sortvar,timevect, ...
'', 0,0,g.erpimopt{:}, 'ampsort', [g.times(time) 0 g.freqs(freq)], 'noshow', 'yes');
else
phsamp = abs(squeeze(tf(freq, time, :)))';
end
% computing ITCs
[ypred alpha(freq, time) Rsq slope(freq, time)] = fastregress(outvar, 20*log10(phsamp));
end
end
sigout = slope;
else
g.times = g.plotvals{1};
g.freqs = g.plotvals{2};
alpha = g.plotvals{3};
sigout = g.plotvals{4};
end
% plot correlation
% ----------------
if ~strcmp('plot', 'no')
if ~isreal(sigout),
if strcmpi(g.plot, 'sigoutp')
sigoutplot = angle(sigout);
else
sigoutplot = abs(sigout);
end
else sigoutplot = sigout;
end
sigouttmp = sigoutplot;
if strcmpi(g.smooth, 'on')
tmpfilt = gauss2d(3,3,.3,.3);
tmpfilt = tmpfilt/sum(sum(tmpfilt));
alpha = conv2(alpha, tmpfilt, 'same');
end
% mask signal out
if g.pmask > 0.5
indices = find( alpha(:) < g.pmask);
sigouttmp = sigoutplot;
sigouttmp(indices) = 0;
elseif g.pmask < 0 % both sides, i.e. 0.01
sigouttmp = sigoutplot;
indices = intersect_bc(find( alpha(:) > -g.pmask), find( alpha(:) < 1+g.pmask));
sigouttmp(indices) = 0;
else
sigouttmp = sigoutplot;
indices = find( alpha(:) > g.pmask);
sigouttmp(indices) = 0;
end
if strcmpi(g.nofig, 'off'), figure; end
switch g.plot
case 'alpha', limits = plotfig(g.times, g.freqs, -log10(alpha), g);
case 'sigout', limits = plotfigsim(g.times, g.freqs, sigoutplot, g);
case { 'sigoutm' 'sigoutp' }, limits = plotfigsim(g.times, g.freqs, sigouttmp, g);
case 'sigoutm2', limits = subplot(1,2,1); plotfigsim(g.times, g.freqs, sigoutplot, g);
limits = subplot(1,2,2); plotfigsim(g.times, g.freqs, sigouttmp, g);
end
end
time = g.times;
freq = g.freqs;
return;
% formula for the normal distribution
% -----------------------------------
function y = norm(mu, sigma, x);
y = 1/sqrt(2) * exp( -(x-mu).^2/(sigma*sigma*2) ) / (sqrt(pi)*sigma);
% get time points
% ---------------
function val = gettimes(timevect, timevar);
if length(timevar) == 2
if timevar(1) > 0
% generate linearly space vector
% ------------------------------
npoints = timevar(1);
trim = timevar(2);
if length(timevect)-2*round(trim/100*length(timevect)) < npoints
npoints = length(timevect)-round(trim/100*length(timevect));
end
fprintf('Generating %d times points trimmed by %1.1f percent\n', npoints, trim);
timer = max(timevect) - min(timevect);
maxt = max(timevect)-timer*trim/100;
mint = min(timevect)+timer*trim/100;
val = linspace(mint,maxt, npoints);
else
% subsample data
% --------------
nsub = -timevar(1);
trim = timevar(2);
len = length(timevect);
trimtimevect = timevect(round(trim/100*len)+1:len-round(trim/100*len));
fprintf('Subsampling by %d and trimming data by %1.1f percent (%d points)\n', nsub, trim, round(trim/100*len));
val = trimtimevect(1:nsub:end);
end
else
val = timevar;
end;
% find closet points in data
oldval = val;
for index = 1:length(val)
[dum ind] = min(abs(timevect-val(index)));
val(index) = timevect(ind);
end
if length(val) < length(unique(val))
disp('Warning: duplicate times, reduce the number of output times');
end
if all(oldval == val)
disp('Debug msg: Time value unchanged by finding closest in data');
end;
% get log scale (for frequency)
% -------------
function val = logscale(a,b,n);
%val = [5 7 9]; return;
val = linspace(log(a), log(b), n);
val = exp(val);
% plot figure
% -----------
function limits = plotfig(times, freqs, vals, g)
icadefs;
imagesc(times, [1:size(vals,1)], vals);
colormap(DEFAULT_COLORMAP);
ticks = linspace(1, size(vals,1), length(freqs));
ticks = ticks(1:4:end);
set(gca, 'ytick', ticks);
set(gca, 'yticklabel', num2str(freqs(1:4:end)', 3))
xlabel('Time (ms)'); ylabel('Freq (Hz)');
if ~isempty(g.limits), caxis(g.limits); end
limits = caxis;
if ~isempty(g.vert)
for vert = g.vert(:)'
hold on; plot([vert vert], [0.001 500], 'k', 'linewidth', 2);
end
end
if strcmpi(g.cbar,'on')
cbar;
end
% plot figure with symmetrical phase
% ---------------------------------
function limits = plotfigsim(times, freqs, vals, g)
icadefs;
imagesc(times, [1:size(vals,1)], vals);
colormap(DEFAULT_COLORMAP);
ticks = linspace(1, size(vals,1), length(freqs));
ticks = ticks(1:4:end);
set(gca, 'ytick', ticks);
set(gca, 'yticklabel', num2str(freqs(1:4:end)', 3))
if ~isempty(g.limits)
caxis(g.limits);
limits = g.limits;
else
clim = max(abs(caxis));
limits = [-clim clim];
caxis(limits);
end
xlabel('Time (ms)'); ylabel('Freq (Hz)');
if ~isempty(g.vert)
for vert = g.vert(:)'
hold on; plot([vert vert], [0.01 500], 'k', 'linewidth', 2);
end
end
if strcmpi(g.cbar,'on')
cbar;
end
return;
% -----------------------------------------
% plot time-freq point (when clicking on it
% -----------------------------------------
function plotpoint(data, sortvar, timevect, freq, timepnts);
figure;
subplot(2,2,1);
erpimage(act_all(:,:),sortvar_all,timepnts, ...
'', 300,10,'erp', 'erpstd', 'caxis',[-1.0 1.0], 'srate',256,'align',352, 'yerplabel', '', erpimopt{:}, ...
'phasesort', [500 0 5]); % aligned to median rt
% plot in polar coordinates phase versus RT
% -----------------------------------------
phaseang2 = [phsangls phsangls-2*pi]; phaseang2 = movav(phaseang2,[], 100);
outvar2 = [outvar outvar]; outvar2 = movav(outvar2,[], 100);
phaseang2 = phaseang2(length(phsangls)/2-50:length(phsangls)+length(phsangls)/2-50);
outvar2 = outvar2 (length(phsangls)/2-50:length(phsangls)+length(phsangls)/2-50);
subplot(2,2,3);
polar(phsangls, outvar, '.');
hold on; polar(phaseang2, outvar2, 'r');
% computing ITC
% -------------
cmplx = outvar .* exp(j*phsangls);
ITCval = mean(cmplx);
angle(ITCval)/pi*180;
abs(ITCval)/mean(outvar);
text(-1300,-1400,sprintf('ITC value: amplitude %3.4f, angle %3.1f\n', abs(ITCval)/mean(outvar), angle(ITCval)/pi*180));
% accumulate 200 values
% ---------------------
alpha = 0.05;
if exist('accarray') ~= 1, accarray = NaN; end
[sigval accarray] = bootstat(outvar/mean(outvar), phsangls, 'res = mean(arg1 .* exp(arg2*complex(0,1)))', ...
'accarray', accarray, 'bootside', 'upper', 'alpha', alpha);
text(-1300,-1600,sprintf('Threshold for significance at %d percent(naccu=200): %3.4f\n', alpha*100, sigval));
title(sprintf('Clust %d corr. theta phase at 500 ms and RT', clust));
% amplitude sorting
% -----------------
figure;
[outdata,outvar,outtrials,limits,axhndls,erp, ...
amps,cohers,cohsig,ampsig,outamps,phsangls,phsamp,sortidx]=erpimage(act_all(:,:),sortvar_all,timepnts, ...
'', 0,0,'erp', 'erpstd', 'caxis', 0.5, 'cbar', ...
'srate',256,'align',352, 'yerplabel', '', erpimopt{:}, 'ampsort', [500 0 5]); % aligned to median rt
close;
subplot(2,2,2);
erpimage(act_all(:,:),sortvar_all,timevect, ...
'', 300,10,'erp', 'erpstd', 'caxis', 0.5, 'cbar', ...
'srate',256,'align',352, 'yerplabel', '', erpimopt{:}, 'ampsort', [500 0 5]); % aligned to median rt
% compute correlation
% -------------------
[ypred alpha Rsq] = myregress(outvar, phsamp);
subplot(2,2,4);
plot(outvar, phsamp, '.');
hold on;
plot(outvar, ypred, 'r');
xlabel('Reaction time');
ylabel('Amplitude');
title(sprintf('Theta amp. at 500 ms vs RT (p=%1.8f, R^2=%3.4f)', alpha, Rsq));
set(gcf, 'position', [336 485 730 540], 'paperpositionmode', 'auto');
setfont(gcf, 'fontsize', 10);
eval(['print -djpeg clust' int2str(clust) 'corrthetart.jpg']);
% set up for command line call
% copy text from plotcorrthetaart
%timevect = times;
sortvar = sortvar_all;
data = act_all;
time = 1
freq = 1
g.times = 0;
g.freqs = 5;
g.erpimopt = {};
Datasets
Models
