function returns estimated AR coefficients using fmristat algorithm FORMAT [rho,Vmhalf,V] = pr_fmristat_ar(res,X,nlags) See http://www.math.mcgill.ca/keith/fmristat/ and fmrilm.m in fmristat package for code, and Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., Morales, F., Evans, A.C. (2002). A general statistical analysis for fMRI data. NeuroImage, 15:1-15 - for description of the algorithm $Id: pr_fmristat_ar.m 301 2004-02-04 03:53:23Z matthewbrett $
0001 function [rho,Vmhalf,V] = pr_fmristat_ar(res,X,nlags) 0002 % function returns estimated AR coefficients using fmristat algorithm 0003 % FORMAT [rho,Vmhalf,V] = pr_fmristat_ar(res,X,nlags) 0004 % 0005 % See http://www.math.mcgill.ca/keith/fmristat/ and 0006 % fmrilm.m in fmristat package for code, and 0007 % Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., Morales, 0008 % F., Evans, A.C. (2002). A general statistical analysis for fMRI 0009 % data. NeuroImage, 15:1-15 - for description of the algorithm 0010 % 0011 % $Id: pr_fmristat_ar.m 301 2004-02-04 03:53:23Z matthewbrett $ 0012 0013 % This is the copyright notice from fmrilm: 0014 %############################################################################ 0015 % COPYRIGHT: Copyright 2002 K.J. Worsley 0016 % Department of Mathematics and Statistics, 0017 % McConnell Brain Imaging Center, 0018 % Montreal Neurological Institute, 0019 % McGill University, Montreal, Quebec, Canada. 0020 % worsley@math.mcgill.ca, liao@math.mcgill.ca 0021 % 0022 % Permission to use, copy, modify, and distribute this 0023 % software and its documentation for any purpose and without 0024 % fee is hereby granted, provided that the above copyright 0025 % notice appear in all copies. The author and McGill University 0026 % make no representations about the suitability of this 0027 % software for any purpose. It is provided "as is" without 0028 % express or implied warranty. 0029 %############################################################################ 0030 0031 if nargin < 2 0032 error('Need covariance and design'); 0033 end 0034 if nargin < 3 0035 nlags = 1; 0036 end 0037 0038 sX = spm_sp('Set', X); 0039 R = spm_sp('r', sX); 0040 0041 nlp1 = nlags+1; 0042 [n nvox] = size(res); 0043 0044 % Bias reduction 0045 M=zeros(nlp1); 0046 for i=1:(nlp1) 0047 for j=1:(nlp1) 0048 Di=(diag(ones(1,n-i+1),i-1)+diag(ones(1,n-i+1),-i+1))/(1+(i==1)); 0049 Dj=(diag(ones(1,n-j+1),j-1)+diag(ones(1,n-j+1),-j+1))/(1+(j==1)); 0050 M(i,j)=trace(R*Di*R*Dj)/(1+(i>1)); 0051 end 0052 end 0053 invM = inv(M); 0054 0055 a = zeros(nlp1,nvox); 0056 for lag = 0:nlags 0057 a(lag+1,:)= sum(res(1:(n-lag),:).*res((lag+1):n,:)); 0058 end 0059 0060 vhat = invM*a; 0061 rho = vhat(2:end,:) ./ (ones(nlags, 1) * vhat(1,:)); 0062 rho = mean(rho,2)'; 0063 0064 if nargout > 1 0065 % Whitening matrix; Appendix A3 Worsley et al (2002) 0066 % Modified according to fmrilm code 0067 [Ainvt posdef] = chol(toeplitz([1 rho])); 0068 nl=size(Ainvt,1); 0069 A=inv(Ainvt'); 0070 Vmhalf = zeros(n,n); 0071 B=ones(n-nl,1)*A(nl,:); 0072 Vmhalf(nl+1:end,:) = spdiags(B,1:nl,n-nl,n); 0073 Vmhalf(1:nl,1:nl) = A; 0074 end 0075 0076 if nargout > 2 0077 % Estimated covariance 0078 V = inv(Vmhalf*Vmhalf'); 0079 end