% 1st-level SPM5 batch stats script SPM5	Rik Henson, Jan 06
% No contrasts, but see SPM2 script on SPM webpage for examples (has not changed)
% (or example SPM5 batch for 2nd-level ANOVAs)

spm('defaults','FMRI')
global defaults
global UFp; UFp = 0.001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nses = 2;		  % Number of sessions
ncon = [4 4];		  % Number of conditions per session
sesdir = {'Ses1';'Ses2'}; % Session directory names

nscan = [740 740];	  % Vector of scans per session

owd = '/myhome';

subnum = [1:18];
subnam = [050045 060001 060002 060003 060004 060005 060013 060014 060015 060016 060017 060018 060019 060023 060024 060025 060035 060036];
dosubs = [1:2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for s = 1:length(dosubs)

    clear SPM

    sub = dosubs(s)
    f = find(subnum==sub);
    cbusub = subnam(f);
	
    cd(owd)

% Replace the next line with a function call to return cell array of onsets,
% durations (could be all 0 for events), and possibly any parametric modulators
% (or to load these from a mat file)
    [sots,durs,pars] = get_sots(sub);

    for ses=1:nses
      for c=1:ncon(ses)
	   onsets{ses,c} = sots{ses,c};
	   if(isempty(onsets{ses,c}))	% Hack to handle any missing events
		onsets{ses,c}=nscan(ses)-3;
  	   end
      end
    end

    swd = sprintf('%s/sub%06d',owd,cbusub);
    cwd = sprintf('%s/%s/Stats',swd,sesdir{1});
    eval(sprintf('!mkdir %s',cwd));
    cd(cwd)

% Specify design 
%==========================================================

    SPM.nscan          = nscan;
    SPM.xBF.name       = 'hrf';
    SPM.xBF.order      = 1;                 
    SPM.xBF.length     = 32; 
    SPM.xBF.T          = 16;     
    SPM.xBF.T0         = 8;     
    SPM.xBF.UNITS      = 'scans';
    SPM.xBF.Volterra   = 1;      

% Trial specification: Onsets, duration (UNITS) and parameters for modulation
%-------------------------------------------------------------
 
    npar = zeros(1,nses);
    for ses=1:nses
	  for c=1:ncon(ses)

	    SPM.Sess(ses).U(c).name      = {sprintf('Con%d',c)};
	    SPM.Sess(ses).U(c).ons       = onsets{c};
	    SPM.Sess(ses).U(c).dur       = durs{c};

	    if isempty(pars{c})
		   SPM.Sess(ses).U(c).P(1).name = 'none';
	    else
                for p=1:size(pars{c},2)
		          SPM.Sess(ses).U(c).P(p).name = sprintf('par %d',p);
 	              SPM.Sess(ses).U(c).P(p).P    = pars{c}(:,p);
	              SPM.Sess(ses).U(c).P(p).h    = 1;	% order (1=linear)
	              SPM.Sess(ses).U(c).P(p).i    = [1 p+1];
                  npar(ses) = npar(ses)+1;
                end
        end
	  end
    end


% design (user specified covariates, eg movement parameters)
%-------------------------------------------------------------


   for ses = 1:nses

	 nusr(ses) = 0;

% Eg if want to model out bad scans (eg spikes...)
%  	spks = load(sprintf('%s/spikes_sub%d_ses%d',swd,sub,wses(ses)));
%	if(~isempty(spks.adout))
%    	 spks = unique(spks.adout(:,1));
% 	 nusr(ses) = length(spks);
%    	 usr = zeros(nscan,nusr(ses));
%    	 for spk=1:nusr(ses)
%		usr(spks(spk),spk)=1;
%    	 end
%	 SPM.Sess(ses).C.C = usr;
%	 tmp=[];
%	 for r=1:nusr(ses)
%	 	tmp{r} = sprintf('spike %d',r);
%	 end
%	 SPM.Sess(ses).C.name = tmp;
%	else
%	 SPM.Sess(ses).C.C = [];
%	 tmp=[];
%	end

	  mname = spm_select('List', fullfile(swd,sesdir{ses}), '^rp_.*\.txt$');
      mpars = load(fullfile(swd,sesdir{ses},mname));

	  SPM.Sess(ses).C.C    = [SPM.Sess(ses).C.C mpars];
	  SPM.Sess(ses).C.name = [tmp {'x' 'y' 'z' 'p' 'r' 'y'}];

	  nusr(ses) = nusr(ses)+6;
     end

% global normalization: OPTIONS:'Scaling'|'None'
%-------------------------------------------------------------
     SPM.xGX.iGXcalc    = 'None';

% low frequency confound: high-pass cutoff (secs) [Inf = no filtering]
%-------------------------------------------------------------
     for ses=1:nses
	   SPM.xX.K(ses).HParam = 128;
     end

% intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
%-------------------------------------------------------------
     SPM.xVi.form       = 'AR(1)';


% specify data: matrix of filenames and TR
%=============================================================

    P=[]; files=cell(1,nses);
    for ses=1:nses
	  files{ses} = spm_select('List', fullfile(swd,sesdir{ses}), '^sw.*\.nii$');
    	for f=1:length(files{ses})    
		P = [P; fullfile(swd,sesdir{ses},files{ses}(f,:))];
	  end
    end
    SPM.xY.P = P;

    SPM.xY.RT          = 1.48;

% Configure design matrix
%=============================================================

    SPM = spm_fmri_spm_ui(SPM);

% Update real/nuisance in design matrix (for nonsphericity correction)
%=============================================================

    if sum(nusr)>0
 	  for ses=1:nses
   		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)))];
      end
 	  SPM.xX.iC = setdiff(SPM.xX.iC,SPM.xX.iG);
    end

% Estimate parameters
%==============================================================
    SPM = spm_spm(SPM);

end