spm_fmri_concatenate.m 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124
  1. function spm_fmri_concatenate(P, scans)
  2. % Adjust an SPM.mat which has concatenated sessions.
  3. % FORMAT spm_post_concatenate(P, scans)
  4. % Session regressors are added and the high-pass filter and non-sphericity
  5. % estimates adjusted as if sessions are separate.
  6. %
  7. % P - filename of the SPM.mat file to adjust
  8. % scans - [1 x n] vector with the original number of scans in each session
  9. %
  10. % The expected workflow is:
  11. %
  12. % 1. Manually specify a GLM with timeseries and onsets concatenated
  13. % 2. Run spm_post_concatenate on the saved SPM.mat.
  14. % 3. Estimate the SPM.mat in the normal way.
  15. %
  16. % Tips:
  17. %
  18. % - The BOLD-response may overhang from one session to the next. To reduce
  19. % this, acquire additional volumes at the end of each session and / or
  20. % add regressors to model the trials at the session borders.
  21. %__________________________________________________________________________
  22. % Copyright (C) 2015-2017 Wellcome Trust Centre for Neuroimaging
  23. % Guillaume Flandin & Peter Zeidman
  24. % $Id: spm_fmri_concatenate.m 7018 2017-02-15 13:36:48Z guillaume $
  25. %-Input parameters
  26. %==========================================================================
  27. %-Get SPM.mat
  28. %--------------------------------------------------------------------------
  29. if ~nargin || isempty(P)
  30. [P, sts] = spm_select(1,'^SPM\.mat$','select SPM.mat');
  31. if ~sts, return; end
  32. end
  33. if iscell(P), P = P{1}; end
  34. %-Get number of scans per session
  35. %--------------------------------------------------------------------------
  36. if nargin < 2
  37. scans = spm_input('scans per session','!+0','r');
  38. end
  39. % Load SPM.mat
  40. %--------------------------------------------------------------------------
  41. SPM = load(P);
  42. % Validate
  43. %--------------------------------------------------------------------------
  44. try
  45. SPM = SPM.SPM;
  46. catch
  47. error('Input file is not a valid SPM.mat.');
  48. end
  49. if ~isfield(SPM,'Sess') || numel(SPM.nscan) ~= 1
  50. error('Input file is not a single session fMRI SPM.mat.');
  51. end
  52. if sum(scans) ~= SPM.nscan
  53. error('Number of scans does not match.');
  54. end
  55. % Create backup
  56. %--------------------------------------------------------------------------
  57. copyfile(P, fullfile(spm_file(P,'fpath'),'SPM_backup.mat'));
  58. %-Session-specific whitening and filtering
  59. %==========================================================================
  60. SPM.nscan = scans;
  61. % Session-specific grand mean scaling (see spm_fmri_spm_ui.m)
  62. %--------------------------------------------------------------------------
  63. % Session effects (see spm_fMRI_design.m)
  64. %--------------------------------------------------------------------------
  65. if numel(SPM.xX.iB) == 1 && SPM.xX.iB == size(SPM.xX.X,2)
  66. Xb = [];
  67. Bn = {};
  68. for s=1:numel(SPM.nscan)
  69. Xb = blkdiag(Xb,ones(SPM.nscan(s),1));
  70. Bn{s} = sprintf('Sn(%i) constant',s);
  71. end
  72. SPM.xX.X = [SPM.xX.X(:,1:end-1) Xb];
  73. SPM.xX.iB = SPM.xX.iB:(SPM.xX.iB+size(Xb,2)-1);
  74. SPM.xX.name = {SPM.xX.name{1:end-1} Bn{:}};
  75. end
  76. % High-pass filter (see spm_spm.m)
  77. %--------------------------------------------------------------------------
  78. s = cumsum([0 SPM.nscan]);
  79. for i=1:numel(SPM.nscan)
  80. K(i) = struct('HParam', SPM.xX.K(1).HParam,...
  81. 'row', s(i) + (1:SPM.nscan(i)),...
  82. 'RT', SPM.xY.RT);
  83. end
  84. SPM.xX.K = spm_filter(K);
  85. % Temporal non-sphericity (see spm_spm.m)
  86. %--------------------------------------------------------------------------
  87. switch lower(SPM.xVi.form)
  88. case {'ar(1)','ar(0.2)'}
  89. SPM.xVi.Vi = spm_Ce('ar',SPM.nscan,0.2);
  90. SPM.xVi.form = 'AR(0.2)';
  91. case 'fast'
  92. SPM.xVi.Vi = spm_Ce('fast',SPM.nscan,SPM.xY.RT);
  93. case {'i.i.d', 'none'}
  94. otherwise
  95. warning('Unhandled temporal non-sphericity.');
  96. end
  97. % Set back nscan as if single session
  98. %--------------------------------------------------------------------------
  99. %SPM.nscan = sum(SPM.nscan);
  100. %-Save
  101. %==========================================================================
  102. save(P, 'SPM', spm_get_defaults('mat.format'));
  103. disp('SPM.mat adjusted for sessions: Please estimate to apply changes.');