function responses = GLMpredictresponses(model,design,tr,numtimepoints,dimdata)
<model> is one of the following:
A where A is X x Y x Z x conditions x time with the timecourse of the
response of each voxel to each condition. XYZ can be collapsed.
{B C} where B is time x 1 with the HRF that is common to all voxels and
all conditions and C is X x Y x Z x conditions with the amplitude of the
response of each voxel to each condition
Note that in both of these cases, the first time point is assumed to be
coincident with condition onset.
<design> is the experimental design. There are three possible cases:
1. A where A is a matrix with dimensions time x conditions.
Each column should be zeros except for ones indicating condition onsets.
(Fractional values in the design matrix are also allowed.)
2. {A1 A2 A3 ...} where each of the A's are like the previous case.
The different A's correspond to different runs, and different runs
can have different numbers of time points.
3. {{C1_1 C2_1 C3_1 ...} {C1_2 C2_2 C3_2 ...} ...} where Ca_b
is a vector of onset times for condition a in run b. Time starts at 0
and is coincident with the acquisition of the first volume. This case
is compatible only with the common-HRF <model>.
<tr> is the sampling rate in seconds
<numtimepoints> is a vector with the number of time points in each run
<dimdata> indicates the dimensionality of the voxels.
A value of 3 indicates X x Y x Z, and a value of 1 indicates XYZ.
Given various inputs, compute the predicted time-series response.
Return:
<responses> as X x Y x Z x time or a cell vector of elements that are
each X x Y x Z x time. The format of <responses> will be a matrix in the
case that <design> is a matrix (case 1) and will be a cell vector in
the other cases (cases 2 and 3).
History:
- 2013/05/12: allow <design> to specify onset times; add <tr>,<numtimepoints> as inputs
- 2013/05/12: update to indicate fractional values in design matrix are allowed.
- 2012/12/03: *** Tag: Version 1.02 ***
- 2012/11/2 - Initial version.0001 function responses = GLMpredictresponses(model,design,tr,numtimepoints,dimdata) 0002 0003 % function responses = GLMpredictresponses(model,design,tr,numtimepoints,dimdata) 0004 % 0005 % <model> is one of the following: 0006 % A where A is X x Y x Z x conditions x time with the timecourse of the 0007 % response of each voxel to each condition. XYZ can be collapsed. 0008 % {B C} where B is time x 1 with the HRF that is common to all voxels and 0009 % all conditions and C is X x Y x Z x conditions with the amplitude of the 0010 % response of each voxel to each condition 0011 % Note that in both of these cases, the first time point is assumed to be 0012 % coincident with condition onset. 0013 % <design> is the experimental design. There are three possible cases: 0014 % 1. A where A is a matrix with dimensions time x conditions. 0015 % Each column should be zeros except for ones indicating condition onsets. 0016 % (Fractional values in the design matrix are also allowed.) 0017 % 2. {A1 A2 A3 ...} where each of the A's are like the previous case. 0018 % The different A's correspond to different runs, and different runs 0019 % can have different numbers of time points. 0020 % 3. {{C1_1 C2_1 C3_1 ...} {C1_2 C2_2 C3_2 ...} ...} where Ca_b 0021 % is a vector of onset times for condition a in run b. Time starts at 0 0022 % and is coincident with the acquisition of the first volume. This case 0023 % is compatible only with the common-HRF <model>. 0024 % <tr> is the sampling rate in seconds 0025 % <numtimepoints> is a vector with the number of time points in each run 0026 % <dimdata> indicates the dimensionality of the voxels. 0027 % A value of 3 indicates X x Y x Z, and a value of 1 indicates XYZ. 0028 % 0029 % Given various inputs, compute the predicted time-series response. 0030 % 0031 % Return: 0032 % <responses> as X x Y x Z x time or a cell vector of elements that are 0033 % each X x Y x Z x time. The format of <responses> will be a matrix in the 0034 % case that <design> is a matrix (case 1) and will be a cell vector in 0035 % the other cases (cases 2 and 3). 0036 % 0037 % History: 0038 % - 2013/05/12: allow <design> to specify onset times; add <tr>,<numtimepoints> as inputs 0039 % - 2013/05/12: update to indicate fractional values in design matrix are allowed. 0040 % - 2012/12/03: *** Tag: Version 1.02 *** 0041 % - 2012/11/2 - Initial version. 0042 0043 % calc 0044 ismatrixcase = ~iscell(design); 0045 dimtime = dimdata + 2; 0046 if iscell(model) 0047 xyzsize = sizefull(model{2},dimdata); 0048 else 0049 xyzsize = sizefull(model,dimdata); 0050 end 0051 0052 % make cell 0053 if ~iscell(design) 0054 design = {design}; 0055 end 0056 0057 % loop over runs 0058 responses = {}; 0059 for p=1:length(design) 0060 0061 % if onset-time case 0062 if iscell(design{p}) 0063 0064 % check that we have the case of common-HRF model 0065 assert(iscell(model)); 0066 0067 % calc 0068 alltimes = linspacefixeddiff(0,tr,numtimepoints(p)); 0069 hrftimes = linspacefixeddiff(0,tr,length(model{1})); 0070 0071 % loop over conditions 0072 temp = zeros(numtimepoints(p),length(design{p})); % this will be time x conditions 0073 for q=1:length(design{p}) 0074 0075 % onset times for qth condition in run p 0076 otimes = design{p}{q}; 0077 0078 % intialize 0079 yvals = 0; 0080 0081 % loop over onset times 0082 for r=1:length(otimes) 0083 0084 % interpolate to find values at the data sampling time points 0085 yvals = yvals + interp1(otimes(r) + hrftimes,model{1}',alltimes,'pchip',0); 0086 0087 end 0088 0089 % record 0090 temp(:,q) = yvals; 0091 0092 end 0093 0094 % weight by the amplitudes 0095 responses{p} = reshape((temp * squish(model{2},dimdata)')',[xyzsize numtimepoints(p)]); % X x Y x Z x time 0096 0097 % if regular matrix case 0098 else 0099 0100 % case of shared HRF 0101 if iscell(model) 0102 0103 % convolve with HRF 0104 temp = conv2(full(design{p}),model{1}); % make full just in case design is sparse 0105 0106 % extract desired subset of time-series 0107 temp = temp(1:numtimepoints(p),:); % time x conditions 0108 0109 % weight by the amplitudes 0110 responses{p} = reshape((temp * squish(model{2},dimdata)')',[xyzsize numtimepoints(p)]); % X x Y x Z x time 0111 0112 % case of individual timecourses 0113 else 0114 0115 % length of each timecourse (L) 0116 len = size(model,dimtime); 0117 0118 % expand design matrix using delta functions 0119 temp = constructstimulusmatrices(design{p}',0,len-1,0); % time x L*conditions 0120 0121 % weight design matrix by the timecourses 0122 responses{p} = reshape((temp * squish(permute(squish(model,dimdata),[3 2 1]),2))',[xyzsize numtimepoints(p)]); % X x Y x Z x time 0123 0124 end 0125 0126 end 0127 0128 end 0129 0130 % undo cell if necessary 0131 if ismatrixcase 0132 responses = responses{1}; 0133 end