autocorr

 method to set autocorrelation types for design
 FORMAT D = autocorr(D, autocorr_type, varargin)
 
 D             - design object
 autocorr_type - autocorrelation type specification, one of
                 'SPM'
                 'fmristat'
                 'none'
 varargin      - parameters defining autocorrelation model.  
                 If autocorr type is 'SPM':
                   varargin{1} should be vector with estimated AR
                      coefficients (default is [0.2])
                   varargin{2} is optional, and is flag; non-zero value
                   specifies voxel-wise covariance estimates (default 1)
                 If autocorr type is 'fmristat'
                   varargin{1} is scalar value for order of fmristat
                   model (default is 1);
 
 $Id$

0001 function D = autocorr(D, autocorr_type, varargin)
0002 % method to set autocorrelation types for design
0003 % FORMAT D = autocorr(D, autocorr_type, varargin)
0004 %
0005 % D             - design object
0006 % autocorr_type - autocorrelation type specification, one of
0007 %                 'SPM'
0008 %                 'fmristat'
0009 %                 'none'
0010 % varargin      - parameters defining autocorrelation model.
0011 %                 If autocorr type is 'SPM':
0012 %                   varargin{1} should be vector with estimated AR
0013 %                      coefficients (default is [0.2])
0014 %                   varargin{2} is optional, and is flag; non-zero value
0015 %                   specifies voxel-wise covariance estimates (default 1)
0016 %                 If autocorr type is 'fmristat'
0017 %                   varargin{1} is scalar value for order of fmristat
0018 %                   model (default is 1);
0019 %
0020 % $Id$
0021   
0022 if nargin < 2
0023   error('Need autocorr type');
0024 end
0025 
0026 if ~is_fmri(D)
0027   warning('Can only set autocorrelation for FMRI design types');
0028   return
0029 end
0030 
0031 % Get design, put into some useful variables
0032 v_f = verbose(D);
0033 SPM = des_struct(D);
0034 nscan = SPM.nscan;
0035 
0036 SPM.xVi.cov_calc = 'summary';
0037 
0038 switch lower(autocorr_type)
0039  case 'fmristat'
0040   % Fit fmristat model AR(n)
0041   if nargin < 3, varargin{1} = 1; end
0042   cVi = varargin{1};
0043   if prod(size(cVi)) > 1
0044     error('Expecting scalar for fmristat order');
0045   end
0046   SPM.xVi.Vi = struct('type', 'fmristat', 'order', cVi);
0047   cVi        = sprintf('fmristat AR(%d)',cVi);
0048   f2cl       = 'V'; % Field to CLear
0049   
0050  case 'spm'
0051   % SPM AR coefficient(s) to be specified
0052   if nargin < 3, varargin{1} = 0.2; end
0053   if nargin < 4, varargin{2} = 1; end
0054   cVi = varargin{1};
0055   if any(cVi > 1 | cVi < 0)
0056     error('Rho estimates should be > 0 and < 1');
0057   end
0058   SPM.xVi.Vi = pr_spm_ce(nscan, cVi);
0059   cVi        = sprintf('AR(%0.1f)',cVi(1));
0060   f2cl       = 'V'; 
0061   if varargin{2}
0062     SPM.xVi.cov_calc = 'vox';
0063   end
0064   
0065  case 'none'        
0066   %  xVi.V is i.i.d
0067   %---------------------------------------------------------------
0068   SPM.xVi.V  = speye(sum(nscan));
0069   cVi        = 'i.i.d';
0070   f2cl       = 'Vi'; 
0071   
0072  otherwise        
0073   error(['Eccentric autocorr type ' autocorr_type]);
0074 end
0075 
0076 % If we've set V, need to clear Vi, because the
0077 % estimate method takes the presence of Vi to mean that
0078 % V can be cleared, with 'redo_covar' flag
0079 % Conversely V needs to be cleared if Vi was estimated
0080 if isfield(SPM.xVi, f2cl)
0081   SPM.xVi = rmfield(SPM.xVi, f2cl);
0082   if v_f, fprintf('Clearing previous %s matrix\n', f2cl); end
0083 end
0084 
0085 % Also: remove previous W matrices
0086 % Either will need to be recalculated or won't be used
0087 if isfield(SPM.xX, 'W')
0088   SPM.xX = rmfield(SPM.xX, 'W');
0089   if v_f, fprintf('Clearing previous W matrix\n'); end
0090 end
0091 
0092 % fill into design
0093 SPM.xVi.form = cVi;
0094 xsDes = struct('Serial_correlations', SPM.xVi.form);
0095 SPM.xsDes = mars_struct('ffillmerge', SPM.xsDes, xsDes);
0096 
0097 % put stuff into object
0098 D = des_struct(D,SPM);