source: misclib/fmri_analysis.m @ 80

Last change on this file since 80 was 80, checked in by tjniskan, 10 years ago
  • Changed the historical "an2_" prefix to "aedes_" in all files. NOTE:

Any script or function relying to Aedes functions will be broken
because of this. Just do a search/replace from "an2_" to "aedes_" in
your files and all should be well...

  • Changed the name of an2_readtab.m to a more informative

aedes_readphasetable.m

File size: 7.8 KB
Line 
1function [maps_out,fdrth] = fmri_analysis(DATA,varargin)
2%
3% This function does a quick fMRI analysis in SPM style. No motion
4% correction or slice timing at this time...
5%
6
7% Defaults
8TR = 2.039;
9contr = [1 0];
10onset = [30 75 120].';
11durat = [15 15 15].';
12smooth_kernel = [2 2 1];
13d_cutoff = [];      % Don't temporally filter data by default
14qFDR = 0.05;        % Default threshold for FDR correction
15globalNorm = false; % Don't do global normalization by default
16omitVols = [];      % Don't omit any volumes from the analysis by default
17show_wbar = true;   % Show waitbar by default
18mask = [];          % Don't mask by default
19
20if rem(length(varargin),2)~=0
21  error('parameters/values missing!')
22end
23
24% Parse varargin
25for ii=1:2:length(varargin)
26  param = lower(varargin{ii});
27  value = varargin{ii+1};
28  switch param
29    case 'tr'
30      TR = value;
31    case {'onset','onsets'}
32      onset = value;
33      onset = onset(:);
34    case {'durat','durations'}
35      durat = value;
36      durat = durat(:);
37    case 'smooth'
38      smooth_kernel = value;
39    case 'contrast'
40      contr = value;
41    case 'fdrth'
42      qFDR = value;
43    case {'hipass','hicutoff'}
44      d_cutoff = value;
45    case 'omitvolumes'
46      omitVols = value;
47    case 'wbar'
48      if value==1
49        show_wbar = true;
50      else
51        show_wbar = false;
52      end
53    case 'mask'
54      if isnumeric(value) || islogical(value)
55        mask = value;
56      else
57        error('Invalid mask!')
58      end
59  end
60end
61
62if length(durat)==1
63  durat = repmat(durat,1,length(onset));
64  durat = durat(:);
65elseif length(durat)~=length(onset)
66  error('Mismatch in the number of elements in onset and duration!')
67end
68
69
70% Check data
71if isstruct(DATA) && isfield(DATA,'FTDATA')
72  data = DATA.FTDATA;
73elseif isnumeric(DATA) && ndims(DATA)>2
74  data=DATA;
75else
76  error('Input data has to be 3D or 4D numeric matrix!')
77end
78
79if ndims(data)==3
80  data = permute(data,[1 2 4 3]);
81end
82
83if isempty(mask)
84  mask = true(size(data,1),size(data,2),size(data,3));
85end
86
87% Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!!
88bf = [
890
90   0.000008876945146
91   0.000387621003097
92   0.003012441669394
93   0.011548214895067
94   0.030056522080168
95   0.061233628208021
96   0.105349995559045
97   0.160158500633323
98   0.221528149631163
99   0.284405279593501
100   0.343761753554314
101   0.395323573307322
102   0.436002661868442
103   0.464045998242252
104   0.478965423658889
105   0.481327079129854
106   0.472473786978796
107   0.454237528221769
108   0.428680153860941
109   0.397883140401152
110   0.363793656845973
111   0.328124931280142
112   0.292303457117353
113   0.257453130446147
114   0.224406054242129
115   0.193730677047637
116   0.165769521647024
117   0.140680556701329
118   0.118477987483032
119   0.099069734765398
120   0.082290068881472
121   0.067926763712502
122   0.055742762013807
123   0.045492744488986
124   0.036935220196777
125   0.029840852217657
126   0.023997740489123
127   0.019214335904674
128   0.015320580887648
129   0.012167779438633
130   0.009627606069476
131   0.007590575511228
132   0.005964217696074
133   0.004671136982604
134   0.003647081075751
135   0.002839102788446
136   0.002203865389816
137   0.001706118264261
138   0.001317352436400
139   0.001014633776781
140   0.000779604140824
141   0.000597636251764
142   0.000457125954290
143   0.000348904854607
144   0.000265756797153
145   0.000202022712087
146   0.000153279812313
147   0.000116082720651
148   0.000087755728843
149   0.000066226941806
150   0.000049896490408
151   0.000037532277385
152];
153
154fprintf(1,'\n******************************************\n');
155fprintf(1,'Starting fMRI analysis.\n');
156if isstruct(DATA) && isfield(DATA,'HDR') && isfield(DATA.HDR,'fpath')
157  filename = [DATA.HDR.fpath,DATA.HDR.fname];
158  fprintf(1,'File: %s\n\n',filename)
159else
160  fprintf(1,'\n');
161end
162
163% Omit volumes from data and stimulus function if requested
164if ~isempty(omitVols)
165 
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:end-1) = durat(1:end-1)-onset(1:end-1)-1;
181    durat(end) = length(tmp)-onset(end)-1;
182  else
183    durat = durat-onset-1;
184  end
185  onset=onset(:);
186  durat=durat(:);
187 
188  fprintf(1,'Skipping requested volumes...\n');
189  fprintf(1,['New onsets: ',num2str(onset(:)'),'\n']);
190  fprintf(1,['New duration: ',num2str(durat(:)'),'\n']);
191  data(:,:,:,omitVols) = [];
192end
193
194% Create stimulus function (32 bin offset)
195k = size(data,4); % Number of scans
196T     = 16;
197dt    = TR/T;
198u     = onset.^0;
199if ~any(durat)
200  u  = u/dt;
201end
202ton       = round(onset*TR/dt) + 32;                    % onsets
203tof       = round(durat*TR/dt) + ton + 1;                       % offset
204sf        = zeros((k*T + 128),size(u,2));
205ton       = max(ton,1);
206tof       = max(tof,1);
207for 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;
214end
215sf        = cumsum(sf);                                 % integrate
216
217% Convolve stimulus with the HRF
218conv_sf = conv(sf,bf);
219
220% Resample the convolved stimulus
221R = conv_sf([0:(k-1)]*16+1+32);
222R=R(:);
223
224% Smooth data
225if all(smooth_kernel)
226  if show_wbar
227    wbh = aedes_calc_wait('Smoothing data...');
228    drawnow
229  end
230  smooth_data = fmri_smooth(data,smooth_kernel);
231else
232  fprintf(1,'fMRI analysis warning: Could not smooth data!\n');
233  smooth_data = data;
234end
235
236if ~isempty(d_cutoff)
237  % Create filtering matrix
238  nCosine = ceil((2*k*TR)/(d_cutoff + 1));
239  S = sqrt(2/k)*cos([1:nCosine]'*pi*([1:k]/k)).';
240  KKT = eye(size(S,1))-2*S*S'+S*S'*S*S';
241else
242  KKT = eye(k);
243  S = 0;
244end
245H = [R ones(size(R))];
246H = H-S*(S'*H); % Filter design matrix
247
248% Calculate parametric map(s)
249nPlanes = size(smooth_data,3);
250nCols = size(smooth_data,2);
251nRows = size(smooth_data,1);
252pmap1 = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3));
253pmap2 = zeros(size(pmap1));
254tmap = zeros(size(pmap1));
255c = repmat(contr,nRows,1);
256HTH = pinv(H'*H);
257R = eye(size(H,1))-H*HTH*H';
258if show_wbar
259  close(wbh);
260end
261if show_wbar
262  wbh = aedes_wbar(0,sprintf('Estimating parameters. Processing plane 0/%d',nPlanes));
263  drawnow
264end
265for ii=1:nPlanes
266 
267  %fprintf(1,'Processing plane %d/%d\n',ii,nPlanes);
268  for kk=1:nCols
269    if show_wbar
270      aedes_wbar(ii/nPlanes,wbh,sprintf('Estimating parameters. Processing plane %d/%d, column %d/%d',ii,nPlanes,kk,nCols));
271    end
272    col_data = squeeze(smooth_data(:,kk,ii,:)).';
273    col_data = col_data-S*(S'*col_data);
274    th=H\col_data;
275    pmap1(:,kk,ii)=th(1,:).';
276    pmap2(:,kk,ii)=th(2,:).';
277   
278    r = col_data-H*th;
279    rr=diag(r'*r);
280    rr=rr(:);
281    sig2 = rr./trace(R*KKT);
282    sig2 = repmat(sig2,1,2);
283    T = diag(c*th)./sqrt(c.*sig2*HTH*H'*KKT*H*HTH*contr');
284    T(find(mask(:,kk,ii)==0))=0;
285    tmap(:,kk,ii)=T.';
286  end
287end
288if show_wbar
289  close(wbh);
290end
291
292if show_wbar
293  wbh = aedes_calc_wait('Calculating threshold...');
294  drawnow
295end
296
297% Calculate effective degrees of freedom
298dof = (trace(R*KKT).^2)/trace(R*KKT*R*KKT);
299
300% p-values
301pval_map = 1-tdist(tmap,dof);
302
303% Perform FDR (False Discovery Rate) correction
304cV = 1;
305
306pValues = pval_map(:);
307tValues = tmap(:);
308
309[pValuesSorted,sortInd] = sort(pValues);
310tValuesSorted = tValues(sortInd);
311nP = length(pValues);
312
313pFDR = [1:nP]'/nP*qFDR/cV; % FDR-correction
314thresFDRind = find(pValuesSorted<=pFDR,1,'last');
315%if isempty(thresFDRind)
316%       fprintf(1,'Warning: No significant values with the FDR-corrected threshold %f\n',...
317%               qFDR);
318%       return
319%end
320if ~isempty(thresFDRind)
321  thresFDR = tValuesSorted(thresFDRind);
322else
323  thresFDR = [];
324end
325
326maps_out.pmap1 = pmap1;
327maps_out.pmap2 = pmap2;
328%maps_out.pmap3 = pmap3;
329maps_out.tmap = tmap;
330
331if show_wbar
332  close(wbh);
333end
334
335if nargout==2
336  fdrth = thresFDR;
337end
338
339if ~isempty(thresFDR)
340  fprintf(1,['FDR threshold at p<',num2str(qFDR),': %.3f\n'],thresFDR)
341else
342  fprintf(1,['No significant voxels over FDR threshold at p<',num2str(qFDR),'!\n'])
343end
344
345fprintf(1,'******************************************\n');
346
347
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008