 Timestamp:
 Oct 12, 2009, 1:20:14 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

misclib/fmri_analysis.m
r83 r94 1 function [maps_out,fdrth ] = fmri_analysis(DATA,varargin)1 function [maps_out,fdrth,H] = fmri_analysis(DATA,varargin) 2 2 % 3 3 % This function does a quick fMRI analysis in SPM style. No motion … … 17 17 show_wbar = true; % Show waitbar by default 18 18 mask = []; % Don't mask by default 19 regr = []; % Don't use additional regressors by default 19 20 20 21 if rem(length(varargin),2)~=0 … … 57 58 error('Invalid mask!') 58 59 end 59 end 60 end 61 62 if length(durat)==1 63 durat = repmat(durat,1,length(onset)); 64 durat = durat(:); 65 elseif length(durat)~=length(onset) 66 error('Mismatch in the number of elements in onset and duration!') 60 case 'regressor' 61 regr = value; 62 end 63 end 64 65 if ~isempty(onset) 66 if length(durat)==1 67 durat = repmat(durat,1,length(onset)); 68 durat = durat(:); 69 elseif length(durat)~=length(onset) 70 error('Mismatch in the number of elements in onset and duration!') 71 end 67 72 end 68 73 … … 83 88 if isempty(mask) 84 89 mask = true(size(data,1),size(data,2),size(data,3)); 90 end 91 92 % Check regressors 93 if ~isempty(regr) 94 if size(regr,1)~=size(data,4) && size(regr,2)~=size(data,4) 95 error('The lengths of the regressors do not match data length!') 96 end 97 if size(regr,1)~=size(data,4) 98 regr = regr.'; 99 end 100 regr = detrend(regr); 85 101 end 86 102 … … 164 180 if ~isempty(omitVols) 165 181 166 % Calculate new onsets and durations 167 ton = onset; 168 tof = onset+durat+1; 169 tmp=zeros(1,size(data,4)); 170 tmp(ton)=1;tmp(tof)=1; 171 sf=cumsum(tmp); 172 sf(omitVols)=[]; 173 tmp=diff([0 sf]); 174 onset = find(tmp==1); 175 durat = find(tmp==1); 176 if isempty(durat) 177 % Block stays up... 178 durat=length(tmp)onset(end)1; 179 elseif length(durat) < length(onset) 180 durat(1:end1) = durat(1:end1)onset(1:end1)1; 181 durat(end) = length(tmp)onset(end)1; 182 else 183 durat = duratonset1; 184 end 185 onset=onset(:); 186 durat=durat(:); 182 fprintf(1,'Skipping requested volumes...\n'); 183 if ~isempty(onset) 184 % Calculate new onsets and durations 185 ton = onset; 186 tof = onset+durat+1; 187 tmp=zeros(1,size(data,4)); 188 tmp(ton)=1;tmp(tof)=1; 189 sf=cumsum(tmp); 190 sf(omitVols)=[]; 191 tmp=diff([0 sf]); 192 onset = find(tmp==1); 193 durat = find(tmp==1); 194 if isempty(durat) 195 % Block stays up... 196 durat=length(tmp)onset(end)1; 197 elseif length(durat) < length(onset) 198 durat(1:end1) = durat(1:end1)onset(1:end1)1; 199 durat(end) = length(tmp)onset(end)1; 200 else 201 durat = duratonset1; 202 end 203 onset=onset(:); 204 durat=durat(:); 205 206 fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']); 207 fprintf(1,['New duration: ',num2str(durat(:)'),'\n']); 208 end 187 209 188 fprintf(1,'Skipping requested volumes...\n');189 fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']);190 fprintf(1,['New duration: ',num2str(durat(:)'),'\n']);210 if ~isempty(regr) 211 regr(omitVols,:)=[]; 212 end 191 213 data(:,:,:,omitVols) = []; 214 192 215 end 193 216 … … 196 219 T = 16; 197 220 dt = TR/T; 198 u = onset.^0; 199 if ~any(durat) 200 u = u/dt; 201 end 202 ton = round(onset*TR/dt) + 32; % onsets 203 tof = round(durat*TR/dt) + ton + 1; % offset 204 sf = zeros((k*T + 128),size(u,2)); 205 ton = max(ton,1); 206 tof = max(tof,1); 207 for j = 1:length(ton) 208 if numel(sf)>ton(j), 209 sf(ton(j),:) = sf(ton(j),:) + u(j,:); 210 end; 211 if numel(sf)>tof(j), 212 sf(tof(j),:) = sf(tof(j),:)  u(j,:); 213 end; 214 end 215 sf = cumsum(sf); % integrate 216 217 % Convolve stimulus with the HRF 218 conv_sf = conv(sf,bf); 219 220 % Resample the convolved stimulus 221 R = conv_sf([0:(k1)]*16+1+32); 222 R=R(:); 221 if ~isempty(onset) 222 u = onset.^0; 223 if ~any(durat) 224 u = u/dt; 225 end 226 ton = round(onset*TR/dt) + 32; % onsets 227 tof = round(durat*TR/dt) + ton + 1; % offset 228 sf = zeros((k*T + 128),size(u,2)); 229 ton = max(ton,1); 230 tof = max(tof,1); 231 for j = 1:length(ton) 232 if numel(sf)>ton(j), 233 sf(ton(j),:) = sf(ton(j),:) + u(j,:); 234 end; 235 if numel(sf)>tof(j), 236 sf(tof(j),:) = sf(tof(j),:)  u(j,:); 237 end; 238 end 239 sf = cumsum(sf); % integrate 240 241 % Convolve stimulus with the HRF 242 conv_sf = conv(sf,bf); 243 244 % Resample the convolved stimulus 245 R = conv_sf([0:(k1)]*16+1+32); 246 R=R(:); 247 else 248 R=[]; 249 end 250 223 251 224 252 % Smooth data … … 243 271 S = 0; 244 272 end 245 H = [R ones(size(R))];273 H = [R regr ones(k,1)]; 246 274 H = HS*(S'*H); % Filter design matrix 247 275 … … 333 361 end 334 362 335 if nargout ==2363 if nargout>=2 336 364 fdrth = thresFDR; 337 365 end
Note: See TracChangeset
for help on using the changeset viewer.