Attachment 'Session2-script.m'
Download 1 % Batch script (linear pipeline) for MEG analysis in SPM5
2 % Rik Henson MRC CBU Sep 08
3 %
4 % Note that would benefit from modularisation (eg for stop/start
5 % running), but wanted all steps in one file to ease exposition
6 %
7 % First start SPM with (at the CBU) "spm 5 eeg" from Linux prompt
8
9 clear
10
11 wd = '/imaging/ms02/MasterclassDemo/Again/'; % !!!REPLACE WITH YOUR DIRECTORY!!!
12 cd(wd)
13
14 % Extra bit for MEG masterclass as maxfiltered files not stored in own
15 % working directory.
16 fifd = '/imaging/ms02/MasterclassDemo/'; % Location of the maxfiltered files
17
18 addpath /imaging/local/spm/spm5/cbu_updates/
19 % Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/spm5_cbu_updates/devel/
20 addpath /imaging/local/meg_misc/
21 % Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/meg_misc/devel/
22
23 %%%%%%%%%% Different types of sensor
24 typ{1} = 'mags'; typ{2} = 'grds'; typ{3} = 'eeg';
25 typ{4} = 'grms'; % RMS of grads; just for sensor-level analyses
26 Ntyp = length(typ);
27 eegflag = [0 0 1 0]; % Binary flag for which typ is EEG
28
29 % Num channels per type (if differs across subjects, can re-specify below)
30 Nchan(1) = 102; Nchan(2) = 204; Nchan(3) = 70;
31 Nchan(4) = 102;
32
33 % Any user-thresholds for rejection for each of above (Inf=none)
34 thr(1) = Inf; thr(2) = Inf; thr(3) = 120;
35 thr(4) = Inf;
36 EOGthr = 120; % EOG threshold (EOG is duplicated in each file)
37
38 %%%%%%%%%% Experiment-specific parameters
39 expwin = [-100 400]; % Epoch window (in ms)
40 expcodes = [1:8]; % Codes in FIF file to extract
41 offset = [ones(1,length(expcodes))*34]; % t=0 offset (eg visual delay)
42 newcodes = [1 1 1 1 2 2 2 2]; % Any new (recoding) of above
43 expcons = [1 0; 0 1; 1 -1; 1/2 1/2]; % Contrasts (on NEW codes)
44 Ncons = size(expcons,1);
45
46 %%%%%%%%%% Raw FIF files on network
47 % (cell of cells for each session of each sub (only one session here))
48 % Eg if instead two sessions per subject:
49 % Fifdata = { { '/megdata/cbu/lfp/meg08_0014/080117/block1_raw.fif';
50 % '/megdata/cbu/lfp/meg08_0014/080117/block2_raw.fif'} };
51 fifdata = {
52 {'/megdata/cbu/lfp/meg08_0014/080117/block4_raw.fif'};
53 {'/megdata/cbu/lfp/meg08_0015/080118/block4_raw.fif'};
54 {'/megdata/cbu/lfp/meg08_0016/080118/block4_raw.fif'};
55 {'/megdata/cbu/lfp/meg08_0038/080131/block4_raw.fif'};
56 {'/megdata/cbu/lfp/meg08_0064/080212/block4_raw.fif'};
57 {'/megdata/cbu/lfp/meg08_0111/080303/block4_raw.fif'};
58 {'/megdata/cbu/lfp/meg08_0112/080303/block4_raw.fif'};
59 {'/megdata/cbu/lfp/meg08_0113/080303/block4_raw.fif'};
60 {'/megdata/cbu/lfp/meg08_0127/080310/block4_raw.fif'};
61 }
62
63
64 %%%%%%%%%% Raw DICOM files on network (one per subject)
65 mridata = {
66 '/mridata/cbu/CBU080769_CBU080769/20080916_095501/Series_002_CBU_MPRAGE';
67 '/mridata/cbu/CBU071255_CBU071255/20071218_150851/Series_002_CBU_MPRAGE';
68 '/mridata/cbu/CBU071151_CBU071151/20071116_140044/Series_002_CBU_MPRAGE';
69 '/mridata/cbu/CBU071229_CBU071229/20071211_112249/Series_002_CBU_MPRAGE';
70 '/mridata/cbu/CBU080034_RIS13667/20080111_102845/Series_002_CBU_MPRAGE';
71 '/mridata/cbu/CBU080292_CBU080292/20080418_161843/Series_002_CBU_MPRAGE';
72 '/mridata/cbu/CBU080358_CBU080358/20080506_151233/Series_003_CBU_MPRAGE';
73 '/mridata/cbu/CBU080413_CBU080413/20080522_124313/Series_002_CBU_MPRAGE';
74 '/mridata/cbu/CBU080802_CBU080802/20080930_090538/Series_002_CBU_MPRAGE';
75 };
76
77
78 %%%%%%%%%% Bad channel names noticed by Operator or post hoc (eeg/meg,sub,ses)
79 badchans{1,1,1} = [1143];
80 badchans{1,2,1} = [1143 0913 1013 0313 0933 2222 0642 1232 1913];
81 badchans{1,3,1} = [1143];
82 badchans{1,4,1} = [1143 2223 0432];
83 badchans{1,5,1} = [1143];
84 badchans{1,6,1} = [1143];
85 badchans{1,7,1} = [1143 2223];
86 badchans{1,8,1} = [1143 1412];
87 badchans{1,9,1} = [1143 2223];
88
89 %%%%%%%%%% Bad EEG channels noticed by Operator or post hoc (in terms of NAME)
90 badchans{2,1,1} = [48 50 73]; % Poor contact
91 badchans{2,2,1} = [7 13 16 24 46 48 57]; % Drifting
92 badchans{2,3,1} = [48]; % Drifting
93 badchans{2,4,1} = [48]; % Drifting
94 badchans{2,5,1} = [15 40 49]; % Poor contact
95 badchans{2,6,1} = [48]; % Drifting
96 badchans{2,7,1} = [48 72]; % No deflection (odd N170 topo)?
97 badchans{2,8,1} = [65 48]; % odd in final topo of N170?
98 badchans{2,9,1} = [20 31 45]; % 45 poor contact; 20+31 highfreq noise
99
100 % (Just happens that any EEG Channels 61 onwards are numbered 65 onwards in our lab!)
101 for sub=1:size(badchans,2)
102 for ses=1:size(badchans,3)
103 f = find(badchans{2,sub,ses}>60);
104 badchans{2,sub,ses}(f) = badchans{2,sub,ses}(f)-4;
105 end
106 end
107
108 %%%%%%%%%%% subject-specific parameters
109 %dosubs = [1:9]; % Subjects to run, eg for subset: dosubs = [5:8];
110 dosubs = [1:8]; % Exclude sub9 because MRI segmentation requires manual re-
111 % positioning first (can try yourself; see MRI section below!)
112 % (in fact for many, eg sub5 too, normalisation is much
113 % better after manual repositioning)
114 Nsub = length(dosubs);
115
116 % Below assumes that two EOG channels (VEOG and HEOG), which will be
117 % appended at end of every data-file. Note that this script does not yet
118 % handle concurrent ECG (see Jason Taylor for how to do this)
119 for s = 1:Nsub
120 subnum = dosubs(s);
121 for m = 1:Ntyp;
122 Nsubchan(s,m) = Nchan(m); % In case different number of channels (eg EEG) for some subjects
123 chan_thr{s,m} = [ones(1,Nsubchan(s,m))*thr(m) EOGthr EOGthr];
124 end
125 end
126
127 VitECapsules = ones(1,max(dosubs)); % Whether MRIs have Vitamin E capsule (all here)
128
129 %%%%%%%%%%% control flags (1=do step; 0=omit step)
130 maxfiltflag = 0; maxmvcompflag = 0; maxtransflag = 0;
131 usemaxmvcompflag = 0; usemaxtransflag = 1;
132 chancheckflag = 0; respmflag = 1; reprocessflag = 1;
133 getmriflag = 0;
134 write2Dflag = [1 0 1 1]; write3Dflag = 0;
135 grpinvertflag = 0; % Whether invert using group priors (after individual inversions)
136 % (Litvak & Friston, 2008, Neuroimage)
137 fusionflag = 0; % Whether want to simultaneously invert (fuse) modalities (typs)
138 % (Henson et al, submitted)
139
140 %%%%%%%%%% General parameters
141 % Maxfilter
142 downsample_maxf = '-ds 4'; % maxfilter downsampling factor (eg 4 to save space)
143 %downsample_maxf = ''; % (if none)
144 format_maxf = '-format float'; % maxfilter data format (native = float32)
145 %format_maxf = '-format short'; % (if to save space)
146 trans_offset_maxf = [0 -13 +6]; % New centre for trans default
147 % (see our maxfilter WIKI)
148
149 % SPM
150 downsample_spm = 100; % spm new sample rate (Hz; should be at
151 % least twice filt_cut_spm below)
152 filt_cut_spm = 40; % Lowpass filter cutoff for SPM (Hz)
153 dformat_spm = 'int16'; % data format for spm (eg to save space)
154 sensor_smooth_spm = [5 5 20]; % Smoothing in SensorSpace-Time [mm mm ms];
155 source_smooth_spm = [12 12 12]; % Smoothing in Source space [mm mm mm]
156
157 %%%%%%%%%%% Initialisation
158 subsphere = [];
159 allsssbad = cell(Nsub,1);
160 nbadchan = cell(Ntyp,Nsub);
161 nrejects = cell(Ntyp,Nsub);
162 nevents = cell(Ntyp,Nsub);
163
164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
165 %%%%%%%%%%% Preprocess MEG/EEG data
166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167
168 for sub = 1:Nsub
169
170 subnum = dosubs(sub);
171 Nses(sub) = size(fifdata{subnum},1);
172
173 cd(wd)
174 swd = sprintf('sub%d',subnum);
175 try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end
176
177 for ses = 1:Nses(sub)
178
179 allsssbad{sub} = cell(1,Nses);
180
181 fnm_stem = sprintf('sub%d_ses%d',subnum,ses);
182
183 if maxfiltflag
184
185 %%%%% Use autobad to detect bad channels from first 20s, and select from output
186 rawfile = fifdata{subnum}{ses};
187 sssfile = sprintf('%s_sss.fif',fnm_stem)
188 logfile = sprintf('%s_autobad.log',fnm_stem)
189 headptsfile = sprintf('%s_headpoints.txt',fnm_stem)
190
191 if exist(sssfile), delete(sssfile), end
192 if exist(logfile), delete(logfile), end
193 [Centre,Radius] = meg_fit_sphere(rawfile,pwd,headptsfile);
194 disp(sprintf('Subject %d (%d): Session %d: Sphere centre [%d %d %d] and radius %d mm',subnum,sub,ses,Centre(1),Centre(2),Centre(3),Radius))
195 subsphere(sub,:) = [Centre Radius];
196
197 disp('Checking first 20s for bad channels...')
198 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d -autobad 20 -skip 21 9999999 -v | tee %s',rawfile,sssfile,Centre(1),Centre(2),Centre(3),logfile))
199
200 badfile = sprintf('%s_badchan.txt',fnm_stem)
201 if exist(badfile), delete(badfile), end
202 eval(sprintf('!cat %s | sed -n ''/Static/p'' | cut -f 5- -d '' '' > %s',logfile,badfile));
203 sssbadchans{sub}{ses} = unique([textread(badfile,'%d')' badchans{1,subnum,ses}]);
204 for bc=1:length(sssbadchans{sub}{ses}); allsssbad{sub}{ses} = [allsssbad{sub}{ses} sprintf(' %04d',sssbadchans{sub}{ses}(bc))]; end
205
206 %%%%% Now mark bad channels and run ST
207 logfile = sprintf('%s_ssst.log',fnm_stem)
208 if exist(sssfile), delete(sssfile), end % Don't need previous file from autobad
209 sssfile = sprintf('%s_ssst.fif',fnm_stem)
210 if exist(logfile), delete(logfile), end
211 posfile = sprintf('%s_mvpos.txt',fnm_stem)
212
213 disp(['Now maxfiltering with bad channels (and tSSS): ' allsssbad{sub}{ses}])
214 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -headpos -hpistep 1000 -hpisubt amp -hp %s %s -v | tee %s',rawfile,sssfile,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},posfile,downsample_maxf,logfile))
215
216 S.logfile = logfile;
217 [max_d(sub),max_r(sub)] = meg_plot_maxfilter_move(S);
218
219 infile = sssfile;
220
221 %%%%% Now mark bad channels and run MVCOMP
222 if maxmvcompflag
223
224 sss2file = sprintf('%s_mvcomp.fif',fnm_stem)
225 logfile = sprintf('%s_mvcomp.log',fnm_stem)
226 if exist(sss2file), delete(sss2file), end
227 if exist(logfile), delete(logfile), end
228
229 disp(['Now compensating movement... (every 200ms)'])
230 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -movecomp -hpistep 200 -hpisubt amp %s -v | tee %s',rawfile,sss2file,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},downsample_maxf,logfile))
231
232 infile = sss2file;
233 end
234
235 %%%%% Now Trans Default (if requested)
236 if maxtransflag
237 trans_origin = Centre(1:3) + trans_offset_maxf;
238 logfile = sprintf('%s_trans.log',fnm_stem)
239 sss3file = sprintf('%s_trans.fif',fnm_stem)
240 if exist(sss3file), delete(sss3file), end
241 if exist(logfile), delete(logfile), end
242
243 disp(['Now trans-ing to DEFAULT ' num2str(round(trans_origin))])
244 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -autobad off -trans default -frame head -origin %d %d %d %s -v -force | tee %s',infile,sss3file,trans_origin(1),trans_origin(2),trans_origin(3),format_maxf,logfile))
245 end
246 end
247
248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
249 %%%%% Olaf's channel checking (have not used for a while; should put in SVN!)
250 if chancheckflag
251 addpath /imaging/olaf/MEG/artefact_scan_tool/
252 global fiffs tasks values criteria sensors datapara figs;
253 ini_check4artefacts; % Initialisation of variables
254 epochlength = 5;
255 amplitude_threshold_mag = 5000;
256 amplitude_threshold_gra = 10000;
257 out_dir = fullfile(wd,'ChanCheck');
258 addfiff(infile,fnm_stem);
259 addtask('maxminamp'); % maximum minus minimum amplitude within epoch
260 addtask('maxmingrad'); % maximum signal gradient (amplitude difference between successive data points in time) with epoch
261 addtask('threshold'); % number of epochs with max-min amplitudes above specified threshold
262 addtask('crosscorr'); % intercorrelation matrix for magnetometers and gradiometers separately
263 addtask('trendamp'); % linear trend in max-min differences across all epochs
264 addtask('trendgrad'); % linear trend in maximum signal gradients across all epochs
265 check4artefacts; % That's where it all happens... lean back and enjoy!
266 end
267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
268
269 %%%%% Read into SPM (writes separate files for MEG and EEG)
270
271 clear S
272 if usemaxtransflag
273 S.Fdata = sprintf('%s_trans.fif',fnm_stem);
274 if maxmvcompflag
275 S.HPIfile = sprintf('%s_mvcomp.fif',fnm_stem); % If trans'ed (hpi lost)
276 else
277 S.HPIfile = sprintf('%s_ssst.fif',fnm_stem); % If trans'ed (hpi lost)
278 end
279 elseif usemaxmvcompflag
280 S.Fdata = sprintf('%s_mvcomp.fif',fnm_stem);
281 S.HPIfile = [];
282 else
283 S.Fdata = sprintf('%s_ssst.fif',fnm_stem);
284 S.HPIfile = [];
285 end
286
287 % Extra bits for MEG masterclass as maxfiltered files not stored in own working directory.
288 S.Fdata = fullfile(fifd, swd, S.Fdata);
289 if ~isempty(S.HPIfile)
290 S.HPIfile = fullfile(fifd, swd, S.HPIfile);
291 end
292 [path, name, ext] = fileparts(S.Fdata);
293 S.Pout = fullfile(wd, swd, sprintf('%s.mat', name)); % New output location
294 % Extra bit ends
295
296 % meg_plot_fifpoints(S); % If want to check points
297 S.Fchannels = 'FIF306_setup.mat';
298 S.veogchan = [62 63]; S.veog_unipolar = 0;
299 S.heogchan = [61 64]; S.heog_unipolar = 0;
300 S.trig_chan = 'STI101';
301 S.twin = [0 Inf];
302 S.eeg_ref = 'average';
303 S.trig_type = 'positive_from_zero';
304 % S.Dtype = 'float32'; % make int16 if files too big
305 S.Fchannels_eeg = fullfile(spm('dir'),'EEGtemplates', sprintf('cbu_meg_%deeg_montage_old.mat',Nsubchan(sub,find(eegflag)))); % use default montage (which simply explains how to display in 2D), prior to new CBU caps in mid-2008
306 % S.Fchannels_eeg = fullfile(wd,swd,sprintf('eeg_montage_sub%02d.mat',subnum)); % use subject-specific (care labels may not be standard)
307
308 if respmflag
309 D = spm_eeg_rdata_FIF(S);
310 basefile = D.fname(1:(end-4));
311
312 %%%%% Downsample, then split Mags and Grads (file too big otherwise)
313
314 clear S; S.D = sprintf('%s-eeg.mat',basefile); disp(S.D)
315 S.Radc_new = downsample_spm;
316 D = spm_eeg_downsample(S);
317
318 clear S; S.D = basefile;
319 S.Radc_new = downsample_spm;
320 D = spm_eeg_downsample(S);
321 basefile = D.fname(1:(end-4));
322
323 clear S; S.D = D.fname
324 D = spm_eeg_splitFIF(S);
325 else
326 basefile = ['d' S.Fdata(1:(end-4))];
327 end
328
329 %%%%% If want to view rawdata (with Danny's tool):
330 % meg_viewdata(D.fname);
331
332 %%%%% Preprocess each data-type
333
334 for typn = 1:Ntyp
335
336 if reprocessflag
337 clear S
338 S.D = sprintf('%s-%s.mat',basefile,typ{typn}); disp(S.D)
339 S.filter.type = 'butterworth';
340 S.filter.band = 'low'; S.filter.PHz = filt_cut_spm;
341 % S.filter.band = 'stop'; S.filter.PHz = [49 51];
342 S.Dtype = dformat_spm;
343 D = spm_eeg_filter(S);
344
345 %%%%%% !!!Insert ICA here if want to remove EOG/ECG artifacts!!!
346 %%%%%% (See Jason Taylor; assume filtered, downsampled data appropriate)
347
348 clear S; S.D = D.fname;
349 S.events.start = expwin(1);
350 S.events.stop = expwin(2);
351 S.events.types = expcodes;
352 S.events.Inewlist = 0;
353 S.events.resynch = offset;
354 D = spm_eeg_epochs(S);
355
356 % re-classify event-codes
357 for e = 1:length(expcodes);
358 fi = find(D.events.code==expcodes(e));
359 if ~isempty(fi)
360 D.events.code(fi) = newcodes(e);
361 end
362 end
363 D.events.types = unique(D.events.code);
364 D.events.Ntypes = length(D.events.types);
365 save(D.fname,'D')
366 clear S; S.D = D.fname;
367
368 else
369 S.D = sprintf('e_fd%s-%s.mat',basefile,typ{typn}); disp(S.D)
370 D = spm_eeg_ldata(S.D);
371 end
372
373 %%%%% Threshold channels (eg EOG)
374
375 nc = D.Nchannels;
376 S.thresholds.External_list = 0;
377 S.thresholds.Check_Threshold = 1;
378 S.thresholds.threshold = chan_thr{sub,typn};
379 S.artefact.weighted = 0;
380 if eegflag(typn)==0
381 % S.channels.Bad = unique([D.channels.Bad strcmpi(D.channels.name....sssbadchans{sub}{ses}]); % If dont trust MaxFilter channel reconstruction
382 else
383 S.channels.Bad = unique([D.channels.Bad badchans{2,subnum,ses}]);
384 end
385 D = spm_eeg_artefact(S);
386
387 nbadchan{typn,sub} = [nbadchan{typn,sub} length(D.channels.Bad)];
388 nrejects{typn,sub} = [nrejects{typn,sub} length(find(D.events.reject))];
389
390 %%%%% Re-reference EEG if removed channels
391 if eegflag(typn) & ~isempty(D.channels.Bad)
392 D = spm_eeg_ldata(D.fname);
393 meg_reavg_ref(D);
394 save(D.fname,'D')
395 end
396 sesnam{typn,ses} = D.fname;
397
398 %%%%% Uncomment if want to average and contrast each session separately
399 %%%%% (as well as after merging below)
400 %
401 % clear S; S.D = D.fname;
402 % D = spm_eeg_average(S);
403 %
404 % clear S; S.D = D.fname;
405 % S.c = expcons;
406 % S.WeightAve = 1;
407 % D = spm_eeg_weight_epochs(S);
408 % nevents{typn,sub} = D.events.repl;
409
410 end % end of typ loop
411 end % end of session loop
412
413 %%%%% Merge sessions, average, contrast
414
415 for typn = 1:Ntyp
416
417 if Nses>1
418 if ~maxtransflag
419 disp(['Concatenating sessions... (NB: merged files will not have correct channel Loc/Orient info unless trans-ed to first or to default!!)']);
420 end
421 clear S; S.D = strvcat(sesnam{typn,:});
422 S.recode = cell(1,Nses);
423 D = spm_eeg_merge(S);
424 clear S; S.D = D.fname;
425 else
426 clear S; S.D = sesnam{typn,1};
427 end
428
429 D = spm_eeg_average(S);
430
431 clear S; S.D = D.fname;
432 S.c = expcons;
433 S.WeightAve = 0;
434 D = spm_eeg_weight_epochs(S);
435 nevents{typn,sub} = D.events.repl;
436 finalnam{typn,sub} = fullfile(D.path,D.fname);
437 end
438
439 if any(write2Dflag)
440 for typn = find(write2Dflag)
441 [pth,nam,ext] = fileparts(finalnam{typn,sub});
442 clear S; S.Fname = [nam ext];
443 S.interpolate_bad = 0;
444 S.n = 32;
445 S.pixsize = 3;
446 S.trialtypes = [1:Ncons];
447 P = spm_eeg_convertmat2ana3D(S);
448 if any(sensor_smooth_spm)
449 if length(sensor_smooth_spm)==1; sensor_smooth_spm = ones(1,3)*sensor_smooth_spm; end
450 for con = 1:size(P,1);
451 [pth,nam,ext] = fileparts(P(con,:));
452 Pout = fullfile(pth,['s' nam ext]);
453 spm_smooth(spm_vol(P(con,:)),Pout,sensor_smooth_spm);
454 Pin = strvcat(P(con,:),Pout);
455 spm_imcalc_ui(Pin,Pout,'((i1+eps).*i2)./(i1+eps)',{[],[],'float32',0}); %Reinsert NaNs
456 SensorTimeImgs{typn}{sub}(con,:) = fullfile(wd,swd,Pout);
457 end
458 else
459 for con = 1:size(P,1);
460 SensorTimeImgs{typn}{sub}(con,:) = fullfile(wd,swd,P(con,:));
461 end
462 end
463 end
464 end
465 end
466
467 cd(wd)
468
469 for typn = 1:Ntyp
470 disp(sprintf('%s...',typ{typn}))
471 for sub = 1:Nsub
472 subnum = dosubs(sub);
473 disp(sprintf('\tSub %d (%d)...',subnum,sub))
474 disp(sprintf('\t\tNevents (per condition): %s',mat2str(nevents{typn,sub})))
475 disp(sprintf('\t\tRejects (per session): %s',mat2str(nrejects{typn,sub})))
476 disp(sprintf('\t\tBadchans (per session): %s',mat2str(nbadchan{typn,sub})))
477 end
478 %%%% Create Grand Mean across subjects
479 clear S; S.P = strvcat(finalnam{typn,:});
480 S.Pout = sprintf('G%d-%s.mat',Nsub,typ{typn})
481 D = spm_eeg_grandmean(S);
482 end
483
484 %return
485
486
487 eval(sprintf('save preproc%d.mat subsphere nrejects nbadchan nevents allsssbad dosubs',length(dosubs)))
488
489 disp('!!!NOW USE SPM-Display-M/EEG button to examine results, eg G8-*.mat files !!!')
490
491 % Once displaying one of the modalities, you can hold shift and select both
492 % conditions 1+2, and click on a posterior right channel to see the M170
493 % difference between faces (1) and scrambled faces (2). You can also select
494 % condition 3, which is the differential ERF/ERP between faces and scrambled,
495 % then press "topography", enter -100ms for initial time, then press and hold
496 % the time slider to see topography of difference at every time point (displayed
497 % in Matlab window) - random noise until ~100ms, when strong differences emerge
498
499
500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
501 %%%%%%%%%%%%%%%%% 2D sensor x 1D Time SPMs
502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
503
504 stdir = fullfile(wd,'SensorTimeSPMs');
505 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
506
507 grpsubs{1} = [1:Nsub];
508 % If want to split subjects, eg into two groups, eg:
509 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
510
511 grpcons = [1:2];
512 Ngrp = length(grpsubs);
513 for typn = find(write2Dflag)
514 outdir = fullfile(stdir,typ{typn})
515 for grp = 1:Ngrp
516 for sub = 1:length(grpsubs{grp});
517 subnum = grpsubs{grp}(sub);
518 imgfiles{grp}{sub} = SensorTimeImgs{typn}{subnum}(grpcons,:);
519 end
520 end
521 meg_batch_anova(imgfiles,outdir);
522 end
523
524 disp('!!!NOW press Results button in SPM!!!')
525
526 % You need to know a bit about the SPM Results interface, but press the Results
527 % button (make sure SPM is in "EEG" mode) (or better still, use the "ns" toolbox
528 % under "Toolbox" button - see WIKI page below) and select the SPM.mat file for
529 % one of the modalities, select the default effects of interest F-contrast
530 % (because we don't care about polarity, at least for Mags and EEG), choose
531 % 0.001 uncorrected. Then press "Volume" to get table of stats (ignore brain image).
532 % You should see, at least in Mags and EEG, frontal and posterior clusters
533 % (simply of different polarity) maximal around 170ms. Not much survives whole-
534 % image correction because only 8 subjects, so few df's ie low power and RFT
535 % conservative) - though effects around 170ms do for corrected extent-level
536 % for Mags and EEG (using "ns" toolbox); GRMS results are very poor though.
537 % Correction for height would also be significant if you used SVC because
538 % of an a priori timewindow of interest around 170ms (again see WIKI page below).
539 %
540 % For more info on interpreting these SPMs, see
541 % http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm
542 %
543 % Also, for an example use with EEG data, see Henson et al (2008), Neuroimage.
544 %
545 % Note that you can do SVC for a priori timewindows of interest (see end of above WIKI)
546 % But one use of this SPMs is to define timewindows of interest (in the absence of
547 % a priori ones) for the subsequent source localisation below...
548 %
549 % (Note that these SPMs can also handle different numbers and locations of channels
550 % across subjects, provided their locations are digitised in same space, because the
551 % data are interpolated to same grid. For MEG data, particularly gradiometers, the
552 % results benefit from using 'trans -default' above, to realign across different
553 % headpositions for different subjects - see Smith et al (in press), HBM
554 % Abstract.
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.