Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_eeg_firstlevel.m
8a186cd82567master

Last Change
Raw File

spm_eeg_firstlevel.m
View Options

function D = spm_eeg_firstlevel(S)
% Convolution modelling for M/EEG data
% FORMAT D = spm_eeg_firstlevel(S)
%
% S - job structure generated by spm_cfg_eeg_firstlevel
%__________________________________________________________________________
% Reference: Litvak V, Jha A, Flandin G, Friston K. Convolution models for
% induced electromagnetic responses. Neuroimage. 2013, 64:388-98
%__________________________________________________________________________
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging
% Vladimir Litvak
% $Id: spm_eeg_firstlevel.m 6146 2014-09-02 11:27:43Z vladimir $
SVNrev = '$Rev: 6146 $';
%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG first level');
D = spm_eeg_load(char(S.sess.D));
statdir = D.path;
try
delete(fullfile(statdir, 'SPM.mat'));
end
job{1}.spm.stats.fmri_design = [];
job{1}.spm.stats.fmri_design.dir = {statdir};
job{1}.spm.stats.fmri_design.timing.units = 'secs';
job{1}.spm.stats.fmri_design.timing.fmri_t0 = 1;
job{1}.spm.stats.fmri_design.volt = S.volt;
job{1}.spm.stats.fmri_design.cvi = 'none';
isTF = strncmpi(D.transformtype,'TF',2);
time = D.time;
dt = mean(diff(time));
nt = D.nsamples;
channels = D.selectchannels(S.channels);
nchan = numel(channels);
nf = D.nfrequencies;
job{1}.spm.stats.fmri_design.timing.RT = dt;
job{1}.spm.stats.fmri_design.sess = rmfield(S.sess, {'D', 'cond', 'convregress', 'multi_conv_reg', 'savereg'});
job{1}.spm.stats.fmri_design.sess.hpf = Inf;
job{1}.spm.stats.fmri_design.sess.nscan = nt;
nc = numel(S.sess.cond);
if nc == 0
job{1}.spm.stats.fmri_design.sess.cond = rmfield(S.sess.cond, 'define');
else
for j = 1:nc
conditionlabels = {};
trigsample = {};
duration = {};
durationsample = {};
if isfield(S.sess.cond(j).define, 'event')
S1 = [];
S1.D = D;
S1.reviewtrials = 0;
S1.save = 0;
S1.trialdef = S.sess.cond(j).define.event;
S1.timewin = [0 0];
[S1.trialdef(:).conditionlabel] = deal(S.sess.cond(j).name);
trl = spm_eeg_definetrial(S1);
trigsample{j} = trl(:, 1);
trig{j} = time(trigsample{j});
duration{j} = 0*trig{j};
durationsample{j} = 0*trig{j};
elseif isfield(S.sess.cond(j).define, 'manual')
if isequal(S.timing.units, 'secs')
trig{j} = S.sess.cond(j).define.manual.onset;
trigsample{j} = D.indsample(trig{j});
duration{j} = S.sess.cond(j).define.manual.duration;
if isempty(duration)
duration{j} = 0*trig{j};
end
durationsample{j} = round(duration{j}*D.fsample);
else
trigsample{j} = S.sess.cond(j).define.manual.onset;
trig{j} = time(trigsample{j});
durationsample{j} = S.sess.cond(j).define.manual.duration;
if isempty(durationsample)
durationsample{j} = 0*trig{j};
end
duration{j} = durationsample{j}/D.fsample;
end
conditionlabels{j} = repmat({S.sess.cond(j).name}, length(trig{j}), 1);
else
error('Unsupported option');
end
job{1}.spm.stats.fmri_design.sess.cond(j).name = S.sess.cond(j).name;
job{1}.spm.stats.fmri_design.sess.cond(j).onset = trig{j} - D.timeonset;
job{1}.spm.stats.fmri_design.sess.cond(j).duration = duration{j};
job{1}.spm.stats.fmri_design.sess.cond(j).tmod = S.sess.cond(j).tmod;
job{1}.spm.stats.fmri_design.sess.cond(j).pmod = S.sess.cond(j).pmod;
job{1}.spm.stats.fmri_design.sess.cond(j).orth = S.sess.cond(j).orth;
end
end
eventchanind = D.indchannel('event');
if ~isempty(eventchanind) && nc
sess = job{1}.spm.stats.fmri_design.sess;
job{1}.spm.stats.fmri_design.sess.cond = ...
rmfield(sess.cond, {'tmod', 'pmod', 'orth'});
job{1}.spm.stats.fmri_design.timing.fmri_t = 1;
job{1}.spm.stats.fmri_design.bases.fir.length = dt;
job{1}.spm.stats.fmri_design.bases.fir.order = 1;
job{1}.spm.stats.fmri_design.sess.regress = [];
job{1}.spm.stats.fmri_design.sess =...
rmfield(job{1}.spm.stats.fmri_design.sess, 'regress');
[dummy, job] = spm_jobman('harvest', job);
spm_run_fmri_spec(job{1}.spm.stats.fmri_design);
load(fullfile(statdir, 'SPM.mat'));
X = SPM.xX.X*SPM.xBF.dt;
X = X(:, 1:(end-1));
aX = any(X', 1);
caX = cumsum(aX);
if isTF
Y = squeeze(mean(D(eventchanind, :, :), 2));
else
Y = D(eventchanind, :);
end
Y = Y-min(Y);
Y = Y/max(Y);
[c, lags] = xcorr(aX, Y, 500, 'coeff');
[dummy, ind] = max(c);
shiftcorr = lags(ind);
if ~spm('CmdLine')
spm_figure('GetWin','Timing check');
clf
subplot(2, 3, 1:3);
plot(lags, c);
xlim(100*[-1 1]);
ntoplot = 5;
i1 = find(caX == 1);
i2 = find(caX == ntoplot);
i3= find((caX-caX(floor(nt/2))) == 1);
i4= find((caX-caX(floor(nt/2))) == ntoplot);
i5= find((caX-caX(end)) == -ntoplot);
i6= find((caX-caX(end)) == 0);
subplot(2, 3, 4);
plot(time(i1(1):i2(1)), Y(i1(1):i2(1)), 'k');
hold on
plot(time(i1(1):i2(1)), X(i1(1):i2(1), :));
subplot(2, 3, 5);
plot(time(i3(1):i4(1)), Y(i3(1):i4(1)), 'k');
hold on
plot(time(i3(1):i4(1)), X(i3(1):i4(1), :));
subplot(2, 3, 6);
plot(time(i5(1):i6(1)), Y(i5(1):i6(1)), 'k');
hold on
plot(time(i5(1):i6(1)), X(i5(1):i6(1), :));
end
delete(fullfile(statdir, 'SPM.mat'));
job{1}.spm.stats.fmri_design.bases = [];
job{1}.spm.stats.fmri_design.sess = sess;
else
shiftcorr = 0;
end
for c = 1:nc
job{1}.spm.stats.fmri_design.sess.cond(c).onset = ...
job{1}.spm.stats.fmri_design.sess.cond(c).onset + 1e-3*S.timing.timewin(1) - shiftcorr*dt;
end
U = struct([]);
for i = 1:numel(S.sess.convregress)
U(i).name = S.sess.convregress(i).name;
U(i).u = S.sess.convregress(i).val;
end
ncr = numel(U);
if ~isequal(S.sess.multi_conv_reg, {''})
temp = load(S.sess.multi_conv_reg{1});
for i = 1:size(temp.R, 2)
U(ncr+i).u = temp.R(:, i);
if isfield(temp, 'names')
U(ncr+i).name = temp.names{i};
else
U(ncr+i).name = ['R' num2str(i)];
end
end
end
ncr = numel(U);
% Create dummy conditions to replace with continuous regressors later
for i = 1:ncr
job{1}.spm.stats.fmri_design.sess.cond(nc+i).name = U(i).name;
job{1}.spm.stats.fmri_design.sess.cond(nc+i).onset = D.time(round(D.nsamples/2));
job{1}.spm.stats.fmri_design.sess.cond(nc+i).duration = 0;
job{1}.spm.stats.fmri_design.sess.cond(nc+i).tmod = 0;
job{1}.spm.stats.fmri_design.sess.cond(nc+i).pmod = struct([]);
job{1}.spm.stats.fmri_design.sess.cond(nc+i).orth = 0;
end
job{1}.spm.stats.fmri_design.timing.fmri_t = S.timing.utime;
bname = char(fieldnames(S.bases));
job{1}.spm.stats.fmri_design.bases = S.bases;
job{1}.spm.stats.fmri_design.bases.(bname).length = 1e-3*diff(sort(S.timing.timewin));
[dummy, job] = spm_jobman('harvest', job);
spm_run_fmri_spec(job{1}.spm.stats.fmri_design);
clear job
load(fullfile(statdir, 'SPM.mat'));
X = SPM.xX.X*SPM.xBF.dt;
nb = size(SPM.xBF.bf, 1);
nr = size(SPM.xBF.bf, 2);
T = round(1/(SPM.xBF.dt*D.fsample));
if ~isempty(U)
if S.timing.timewin(1)<0
pad = round(abs(1e-3*S.timing.timewin(1)*D.fsample));
else
pad = 0;
end
for i = 1:ncr
u = U(i).u;
u = u(:);
if length(u)~= D.nsamples
error('Convolution regressor should be the same length as the data.');
end
u = [u(:); zeros(pad, 1)];
if T~=1
u = interp1(1:length(u), u, 1:1/T:length(u));
end
U(i).u = u(:);
U(i).dt = SPM.xBF.dt;
if ~iscell(U(i).name)
U(i).name = {U(i).name};
end
end
%-Convolve stimulus functions with basis functions
%----------------------------------------------------------------------
[Xu] = spm_Volterra(U, SPM.xBF.bf);
%-Resample regressors
%----------------------------------------------------------------------
Xu = Xu((0:(nt + pad - 1))*T+1, :);
Xu = Xu((1+pad):end, :);
X(:, (nc*nr+1):((nc+ncr)*nr)) = Xu*SPM.xBF.dt;
end
SPM.xX.X = X./SPM.xBF.dt;
%-Save SPM.mat
%--------------------------------------------------------------------------
fprintf('%-40s: ','Saving SPM configuration') %-#
save(fullfile(statdir, 'SPM.mat'), 'SPM', spm_get_defaults('mat.format'));
fprintf('%30s\n','...SPM.mat saved') %-#
if S.sess.hpf
K = [];
K.RT = dt;
K.row = SPM.Sess.row;
K.HParam = S.sess.hpf;
else
K = 1;
end
xX = spm_sp('Set', spm_filter(K, X)); % filter design
if ~isfield(xX,'pX')
xX.pX = spm_sp('x-',xX);
end
if ~isempty(SPM.Sess.U)
label = {};
for i = 1:numel(SPM.Sess.U)
label = [label SPM.Sess.U(i).name];
end
ne = numel(label);
if isTF
Dout = clone(D, spm_file(D.fname, 'prefix', S.prefix), [nchan nf nb ne]);
else
Dout = clone(D, spm_file(D.fname, 'prefix', S.prefix), [nchan nb ne]);
end
Dout = fsample(Dout, 1/SPM.xBF.dt);
Dout = timeonset(Dout, 1e-3*S.timing.timewin(1));
Dout = chanlabels(Dout, ':', D.chanlabels(channels));
Dout = conditions(Dout, ':', label);
Dout = type(Dout, 'evoked');
spm_progress_bar('Init', nchan, 'channels done');
if nchan > 100, Ibar = floor(linspace(1, nchan, 100));
else Ibar = 1:nchan; end
for c = 1:nchan
W = ones(nt, 1);
W(D.badsamples(channels(c), ':', 1)) = exp(-256);
W = spdiags(W, 0, nt, nt);
axX = spm_sp('Set', spm_filter(K, W*X));
if ~isfield(axX,'pX')
axX.pX = spm_sp('x-',axX);
end
Y = reshape(D(channels(c), :, :, :), nf, nt);
Y = spm_filter(K, W*Y');
B = axX.pX*Y;
for i = 1:ne
xY = SPM.xBF.bf*B((i-1)*nr+(1:nr), :);
if isTF
Dout(c, :, :, i) = shiftdim(xY', -1);
else
Dout(c, :, i) = xY';
end
end
if ismember(c, Ibar), spm_progress_bar('Set', c); end
end
spm_progress_bar('Clear');
%-Save
%--------------------------------------------------------------------------
save(Dout);
else
Dout = [];
nr = 0;
ne = 0;
end
if ~isempty(SPM.Sess.C.C) && S.sess.savereg
rlabel = {};
for i = 1:numel(SPM.Sess.C)
rlabel = [rlabel SPM.Sess.C(i).name];
end
ng = numel(rlabel);
if ~isempty(Dout)
preprefix = 'R';
else
preprefix = '';
end
if isTF
Dr = clone(D, spm_file(D.fname, 'prefix', [S.prefix preprefix]), [nchan nf 1 ng]);
else
Dr = clone(D, spm_file(D.fname, 'prefix', [S.prefix preprefix]), [nchan 1 ng]);
end
Dr = fsample(Dr, 0);
Dr = timeonset(Dr, 0);
Dr = chanlabels(Dr, ':', D.chanlabels(channels));
Dr = conditions(Dr, ':', rlabel);
Dr = type(Dr, 'evoked');
spm_progress_bar('Init', nchan, 'channels done');
if nchan > 100, Ibar = floor(linspace(1, nchan, 100));
else Ibar = 1:nchan; end
for c = 1:nchan
W = ones(nt, 1);
W(D.badsamples(channels(c), ':', 1)) = exp(-256);
W = spdiags(W, 0, nt, nt);
axX = spm_sp('Set', spm_filter(K, W*X));
if ~isfield(axX,'pX')
axX.pX = spm_sp('x-',axX);
end
Y = reshape(D(channels(c), :, :, :), nf, nt);
Y = spm_filter(K, W*Y');
B = axX.pX*Y;
for i = 1:ng
xY = B(ne*nr+i, :);
if isTF
Dr(c, :, :, i) = xY;
else
Dr(c, :, i) = xY;
end
end
if ismember(c, Ibar), spm_progress_bar('Set', c); end
end
spm_progress_bar('Clear');
%-Save
%--------------------------------------------------------------------------
save(Dr);
else
Dr = [];
end
if ~isempty(Dout)
D = Dout;
else
D = Dr;
end
c4science · Help