1 % 1st-level SPM5 batch stats script SPM5	Rik Henson, Jan 06
   2 % No contrasts, but see SPM2 script on SPM webpage for examples (has not changed)
   3 % (or example SPM5 batch for 2nd-level ANOVAs)
   4 
   5 spm('defaults','FMRI')
   6 global defaults
   7 global UFp; UFp = 0.001;
   8 
   9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  10 
  11 nses = 2;		  % Number of sessions
  12 ncon = [4 4];		  % Number of conditions per session
  13 sesdir = {'Ses1';'Ses2'}; % Session directory names
  14 
  15 nscan = [740 740];	  % Vector of scans per session
  16 
  17 owd = '/myhome';
  18 
  19 subnum = [1:18];
  20 subnam = [050045 060001 060002 060003 060004 060005 060013 060014 060015 060016 060017 060018 060019 060023 060024 060025 060035 060036];
  21 dosubs = [1:2];
  22 
  23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  24 
  25 for s = 1:length(dosubs)
  26 
  27     clear SPM
  28 
  29     sub = dosubs(s)
  30     f = find(subnum==sub);
  31     cbusub = subnam(f);
  32 	
  33     cd(owd)
  34 
  35 % Replace the next line with a function call to return cell array of onsets,
  36 % durations (could be all 0 for events), and possibly any parametric modulators
  37 % (or to load these from a mat file)
  38     [sots,durs,pars] = get_sots(sub);
  39 
  40     for ses=1:nses
  41       for c=1:ncon(ses)
  42 	   onsets{ses,c} = sots{ses,c};
  43 	   if(isempty(onsets{ses,c}))	% Hack to handle any missing events
  44 		onsets{ses,c}=nscan(ses)-3;
  45   	   end
  46       end
  47     end
  48 
  49     swd = sprintf('%s/sub%06d',owd,cbusub);
  50     cwd = sprintf('%s/%s/Stats',swd,sesdir{1});
  51     eval(sprintf('!mkdir %s',cwd));
  52     cd(cwd)
  53 
  54 % Specify design 
  55 %==========================================================
  56 
  57     SPM.nscan          = nscan;
  58     SPM.xBF.name       = 'hrf';
  59     SPM.xBF.order      = 1;                 
  60     SPM.xBF.length     = 32; 
  61     SPM.xBF.T          = 16;     
  62     SPM.xBF.T0         = 8;     
  63     SPM.xBF.UNITS      = 'scans';
  64     SPM.xBF.Volterra   = 1;      
  65 
  66 % Trial specification: Onsets, duration (UNITS) and parameters for modulation
  67 %-------------------------------------------------------------
  68  
  69     npar = zeros(1,nses);
  70     for ses=1:nses
  71 	  for c=1:ncon(ses)
  72 
  73 	    SPM.Sess(ses).U(c).name      = {sprintf('Con%d',c)};
  74 	    SPM.Sess(ses).U(c).ons       = onsets{c};
  75 	    SPM.Sess(ses).U(c).dur       = durs{c};
  76 
  77 	    if isempty(pars{c})
  78 		   SPM.Sess(ses).U(c).P(1).name = 'none';
  79 	    else
  80                 for p=1:size(pars{c},2)
  81 		          SPM.Sess(ses).U(c).P(p).name = sprintf('par %d',p);
  82  	              SPM.Sess(ses).U(c).P(p).P    = pars{c}(:,p);
  83 	              SPM.Sess(ses).U(c).P(p).h    = 1;	% order (1=linear)
  84 	              SPM.Sess(ses).U(c).P(p).i    = [1 p+1];
  85                   npar(ses) = npar(ses)+1;
  86                 end
  87         end
  88 	  end
  89     end
  90 
  91 
  92 % design (user specified covariates, eg movement parameters)
  93 %-------------------------------------------------------------
  94 
  95 
  96    for ses = 1:nses
  97 
  98 	 nusr(ses) = 0;
  99 
 100 % Eg if want to model out bad scans (eg spikes...)
 101 %  	spks = load(sprintf('%s/spikes_sub%d_ses%d',swd,sub,wses(ses)));
 102 %	if(~isempty(spks.adout))
 103 %    	 spks = unique(spks.adout(:,1));
 104 % 	 nusr(ses) = length(spks);
 105 %    	 usr = zeros(nscan,nusr(ses));
 106 %    	 for spk=1:nusr(ses)
 107 %		usr(spks(spk),spk)=1;
 108 %    	 end
 109 %	 SPM.Sess(ses).C.C = usr;
 110 %	 tmp=[];
 111 %	 for r=1:nusr(ses)
 112 %	 	tmp{r} = sprintf('spike %d',r);
 113 %	 end
 114 %	 SPM.Sess(ses).C.name = tmp;
 115 %	else
 116 %	 SPM.Sess(ses).C.C = [];
 117 %	 tmp=[];
 118 %	end
 119 
 120 	  mname = spm_select('List', fullfile(swd,sesdir{ses}), '^rp_.*\.txt$');
 121       mpars = load(fullfile(swd,sesdir{ses},mname));
 122 
 123 	  SPM.Sess(ses).C.C    = [SPM.Sess(ses).C.C mpars];
 124 	  SPM.Sess(ses).C.name = [tmp {'x' 'y' 'z' 'p' 'r' 'y'}];
 125 
 126 	  nusr(ses) = nusr(ses)+6;
 127      end
 128 
 129 % global normalization: OPTIONS:'Scaling'|'None'
 130 %-------------------------------------------------------------
 131      SPM.xGX.iGXcalc    = 'None';
 132 
 133 % low frequency confound: high-pass cutoff (secs) [Inf = no filtering]
 134 %-------------------------------------------------------------
 135      for ses=1:nses
 136 	   SPM.xX.K(ses).HParam = 128;
 137      end
 138 
 139 % intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
 140 %-------------------------------------------------------------
 141      SPM.xVi.form       = 'AR(1)';
 142 
 143 
 144 % specify data: matrix of filenames and TR
 145 %=============================================================
 146 
 147     P=[]; files=cell(1,nses);
 148     for ses=1:nses
 149 	  files{ses} = spm_select('List', fullfile(swd,sesdir{ses}), '^sw.*\.nii$');
 150     	for f=1:length(files{ses})    
 151 		P = [P; fullfile(swd,sesdir{ses},files{ses}(f,:))];
 152 	  end
 153     end
 154     SPM.xY.P = P;
 155 
 156     SPM.xY.RT          = 1.48;
 157 
 158 % Configure design matrix
 159 %=============================================================
 160 
 161     SPM = spm_fmri_spm_ui(SPM);
 162 
 163 % Update real/nuisance in design matrix (for nonsphericity correction)
 164 %=============================================================
 165 
 166     if sum(nusr)>0
 167  	  for ses=1:nses
 168    		SPM.xX.iG = [SPM.xX.iG  [1:nusr(ses)] + (sum(ncon(1:ses))+sum(npar(1:ses)))*SPM.xBF.order + sum(nusr(1:(ses-1)))];
 169       end
 170  	  SPM.xX.iC = setdiff(SPM.xX.iC,SPM.xX.iG);
 171     end
 172 
 173 % Estimate parameters
 174 %==============================================================
 175     SPM = spm_spm(SPM);
 176 
 177 end