% Example 2: Inspect the model fits obtained by analyzePRF
0001 %% Example 2: Inspect the model fits obtained by analyzePRF 0002 0003 %% Download dataset (if necessary) and add analyzePRF to the MATLAB path 0004 0005 setup; 0006 0007 %% Load and analyze the data (this follows the example1.m script) 0008 0009 % Start parallel MATLAB to speed up execution 0010 if matlabpool('size')==0 0011 matlabpool open; 0012 end 0013 0014 % Load in the data 0015 load('exampledataset.mat'); 0016 0017 % Upsample the data to match the stimulus rate 0018 data = tseriesinterp(data,2,1,2); 0019 0020 % Analyze the data 0021 results = analyzePRF(stimulus,data,1,struct('seedmode',[0 1],'display','off')); 0022 %% 0023 0024 %% Perform some setup 0025 0026 % Define some variables 0027 res = [100 100]; % row x column resolution of the stimuli 0028 resmx = 100; % maximum resolution (along any dimension) 0029 hrf = results.options.hrf; % HRF that was used in the model 0030 degs = results.options.maxpolydeg; % vector of maximum polynomial degrees used in the model 0031 0032 % Pre-compute cache for faster execution 0033 [d,xx,yy] = makegaussian2d(resmx,2,2,2,2); 0034 0035 % Prepare the stimuli for use in the model 0036 stimulusPP = {}; 0037 for p=1:length(stimulus) 0038 stimulusPP{p} = squish(stimulus{p},2)'; % this flattens the image so that the dimensionality is now frames x pixels 0039 stimulusPP{p} = [stimulusPP{p} p*ones(size(stimulusPP{p},1),1)]; % this adds a dummy column to indicate run breaks 0040 end 0041 0042 % Define the model function. This function takes parameters and stimuli as input and 0043 % returns a predicted time-series as output. Specifically, the variable <pp> is a vector 0044 % of parameter values (1 x 5) and the variable <dd> is a matrix with the stimuli (frames x pixels). 0045 % Although it looks complex, what the function does is pretty straightforward: construct a 0046 % 2D Gaussian, crop it to <res>, compute the dot-product between the stimuli and the 0047 % Gaussian, raise the result to an exponent, and then convolve the result with the HRF, 0048 % taking care to not bleed over run boundaries. 0049 modelfun = @(pp,dd) conv2run(posrect(pp(4)) * (dd*[vflatten(placematrix(zeros(res),makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0) / (2*pi*abs(pp(3))^2))); 0]) .^ posrect(pp(5)),hrf,dd(:,prod(res)+1)); 0050 0051 % Construct projection matrices that fit and remove the polynomials. 0052 % Note that a separate projection matrix is constructed for each run. 0053 polymatrix = {}; 0054 for p=1:length(degs) 0055 polymatrix{p} = projectionmatrix(constructpolynomialmatrix(size(data{p},2),0:degs(p))); 0056 end 0057 0058 %% Inspect the data and the model fit 0059 0060 % Which voxel should we inspect? Let's inspect the second voxel. 0061 vx = 2; 0062 0063 % For each run, collect the data and the model fit. We project out polynomials 0064 % from both the data and the model fit. This deals with the problem of 0065 % slow trends in the data. 0066 datats = {}; 0067 modelts = {}; 0068 for p=1:length(data) 0069 datats{p} = polymatrix{p}*data{p}(vx,:)'; 0070 modelts{p} = polymatrix{p}*modelfun(results.params(1,:,vx),stimulusPP{p}); 0071 end 0072 0073 % Visualize the results 0074 figure; hold on; 0075 set(gcf,'Units','points','Position',[100 100 1000 100]); 0076 plot(cat(1,datats{:}),'r-'); 0077 plot(cat(1,modelts{:}),'b-'); 0078 straightline(300*(1:4)+.5,'v','g-'); 0079 xlabel('Time (s)'); 0080 ylabel('BOLD signal'); 0081 ax = axis; 0082 axis([.5 1200+.5 ax(3:4)]); 0083 title('Time-series data');