function noisereg = analyzePRFcomputeGLMdenoiseregressors(stimulus,data,tr) <stimulus>,<data>,<tr> are from analyzePRF this is an internal function called by analyzePRF.m. this function constructs a GLM design matrix based on <stimulus> and then uses GLMdenoisedata.m to determine the optimal set of noise regressors. these regressors are returned in <noisereg>. in case you want to see the results of GLMdenoisedata.m, the figure output and the results of GLMdenoisedata.m are saved to temporary files (as reported to the command window).
0001 function noisereg = analyzePRFcomputeGLMdenoiseregressors(stimulus,data,tr) 0002 0003 % function noisereg = analyzePRFcomputeGLMdenoiseregressors(stimulus,data,tr) 0004 % 0005 % <stimulus>,<data>,<tr> are from analyzePRF 0006 % 0007 % this is an internal function called by analyzePRF.m. this function constructs a 0008 % GLM design matrix based on <stimulus> and then uses GLMdenoisedata.m to determine 0009 % the optimal set of noise regressors. these regressors are returned in <noisereg>. 0010 % 0011 % in case you want to see the results of GLMdenoisedata.m, the figure output and 0012 % the results of GLMdenoisedata.m are saved to temporary files (as reported to the 0013 % command window). 0014 0015 % internal constants 0016 corrthresh = .9; % used in determining which apertures are the same 0017 0018 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure out GLM design matrix 0019 0020 % concatenate everything and drop the dummy column (X is pixels x frames) 0021 X = catcell(1,stimulus)'; 0022 X(end,:) = []; 0023 0024 % figure out where blanks are (logical vector, 1 x frames) 0025 blanks = all(X==0,1); 0026 0027 % normalize each frame (in preparation for computing correlations) 0028 X = unitlength(X,1,[],0); 0029 X(:,blanks) = 0; 0030 0031 % initialize the result (this will grow in size on-the-fly) 0032 glmdesign = zeros(size(X,2),0); 0033 0034 % do the loop 0035 wh = find(~blanks); 0036 cnt = 1; 0037 while ~isempty(wh) 0038 ix = wh(1); % pick the first one to process 0039 corrs = X(:,ix)' * X; % compute correlation with all frames 0040 spots = find(corrs > corrthresh); % any frame with r > corrthresh counts as the same 0041 %%% fprintf('cnt=%d: numspots=%d\n',cnt,length(spots)); 0042 glmdesign(spots,cnt) = 1; % add to design matrix 0043 X(:,spots) = 0; % blank out (since we're done with those columns) 0044 blanks(spots) = 1; % update the list of blanks 0045 wh = find(~blanks); 0046 cnt = cnt + 1; 0047 end 0048 0049 % finally, un-concatenate the results 0050 glmdesign = splitmatrix(glmdesign,1,cellfun(@(x) size(x,1),stimulus)); 0051 0052 % clean up 0053 clear X; 0054 0055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% run GLMdenoise to get the noise regressors 0056 0057 % what directory to save results to? 0058 glmdenoisedir = [tempname]; 0059 assert(mkdir(glmdenoisedir)); 0060 0061 % call GLMdenoise 0062 fprintf('using GLMdenoise figure directory %s\n',[glmdenoisedir '/GLMdenoisefigures']); 0063 0064 % prep 0065 hrfmodel = 'assume'; 0066 hrfknobs = []; 0067 0068 % run it 0069 results = GLMdenoisedata(glmdesign,data,tr,tr,hrfmodel,hrfknobs, ... 0070 struct('numboots',0), ... 0071 [glmdenoisedir '/GLMdenoisefigures']); 0072 0073 % get the noise regressors 0074 noisereg = cellfun(@(x) x(:,1:results.pcnum),results.pcregressors,'UniformOutput',0); 0075 0076 % save 'results' to a file 0077 file0 = [glmdenoisedir '/GLMdenoise.mat']; 0078 fprintf('saving GLMdenoise results to %s (in case you want them).\n',file0); 0079 results = rmfield(results,{'pcweights' 'models' 'modelse'}); % remove some boring fields 0080 save(file0,'results','noisereg','glmdesign'); 0081 0082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% quick visualization 0083 0084 for p=1:length(glmdesign) 0085 figureprep([100 100 700 700]); hold on; 0086 imagesc(glmdesign{p}); colormap(gray); 0087 axis image tight; 0088 set(gca,'YDir','reverse'); 0089 title(sprintf('run %02d',p)); 0090 figurewrite(sprintf('glmdesign%02d',p),[],[],glmdenoisedir); 0091 end