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.