source: misclib/fmri_analysis.m @ 64

Last change on this file since 64 was 64, checked in by tjniskan, 11 years ago
  • Fixed trend estimation in "show voxel time series" window with EPI

data

  • the fmri_analysis function no longer removes the reference image of

the VNMR EPIs. The basic_fmri_analysis -plugin however removes the
reference image.

M misclib/fmri_analysis.m
M an2_revision.m
M aedes.m
M plugins/basic_fmri_analysis.m

File size: 4.9 KB
Line 
1function maps_out = 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];
13 
14
15if rem(length(varargin),2)~=0
16  error('parameters/values missing!')
17end
18
19% Parse varargin
20for ii=1:2:length(varargin)
21  param = lower(varargin{ii});
22  value = varargin{ii+1};
23  switch param
24    case 'TR'
25      TR = value;
26    case {'onset','onsets'}
27      onset = value;
28      onset = onset(:);
29    case {'durat','durations'}
30      durat = value;
31      durat = durat(:);
32    case 'smooth'
33      smooth_kernel = value;
34    case 'contrast'
35      contr = value;
36  end
37end
38
39if length(durat)==1
40  durat = repmat(durat,1,length(onset));
41elseif length(durat)~=length(onset)
42  error('Mismatch in the number of elements in onset and duration!')
43end
44
45
46% Check data
47if isstruct(DATA) && isfield(DATA,'FTDATA')
48%   if isfield(DATA,'PROCPAR') && isfield(DATA.PROCPAR,'readres')
49%     data = DATA.FTDATA(:,:,:,2:end);
50%   else
51    data = DATA.FTDATA;
52%   end
53elseif isnumeric(DATA) && ndims(DATA)>2
54  data=DATA;
55else
56  error('Input data has to be 3D or 4D numeric matrix!')
57end
58
59if ndims(data)==3
60  data = permute(data,[1 2 4 3]);
61end
62
63% Load Rat HRF (8 seconds long)
64bf = [
650
66   0.000008876945146
67   0.000387621003097
68   0.003012441669394
69   0.011548214895067
70   0.030056522080168
71   0.061233628208021
72   0.105349995559045
73   0.160158500633323
74   0.221528149631163
75   0.284405279593501
76   0.343761753554314
77   0.395323573307322
78   0.436002661868442
79   0.464045998242252
80   0.478965423658889
81   0.481327079129854
82   0.472473786978796
83   0.454237528221769
84   0.428680153860941
85   0.397883140401152
86   0.363793656845973
87   0.328124931280142
88   0.292303457117353
89   0.257453130446147
90   0.224406054242129
91   0.193730677047637
92   0.165769521647024
93   0.140680556701329
94   0.118477987483032
95   0.099069734765398
96   0.082290068881472
97   0.067926763712502
98   0.055742762013807
99   0.045492744488986
100   0.036935220196777
101   0.029840852217657
102   0.023997740489123
103   0.019214335904674
104   0.015320580887648
105   0.012167779438633
106   0.009627606069476
107   0.007590575511228
108   0.005964217696074
109   0.004671136982604
110   0.003647081075751
111   0.002839102788446
112   0.002203865389816
113   0.001706118264261
114   0.001317352436400
115   0.001014633776781
116   0.000779604140824
117   0.000597636251764
118   0.000457125954290
119   0.000348904854607
120   0.000265756797153
121   0.000202022712087
122   0.000153279812313
123   0.000116082720651
124   0.000087755728843
125   0.000066226941806
126   0.000049896490408
127   0.000037532277385
128];
129
130% Create stimulus function (32 bin offset)
131k = size(data,4); % Number of scans
132T     = 16;
133dt    = TR/T;
134u     = onset.^0;
135if ~any(durat)
136  u  = u/dt;
137end
138ton       = round(onset*TR/dt) + 32;                    % onsets
139tof       = round(durat*TR/dt) + ton + 1;                       % offset
140sf        = zeros((k*T + 128),size(u,2));
141ton       = max(ton,1);
142tof       = max(tof,1);
143for j = 1:length(ton)
144  if numel(sf)>ton(j),
145    sf(ton(j),:) = sf(ton(j),:) + u(j,:);
146  end;
147  if numel(sf)>tof(j),
148    sf(tof(j),:) = sf(tof(j),:) - u(j,:);
149  end;
150end
151sf        = cumsum(sf);                                 % integrate
152
153% Convolve stimulus with the HRF
154conv_sf = conv(sf,bf);
155
156% Resample the convolved stimulus
157R = conv_sf([0:(k-1)]*16+1+32);
158R=R(:);
159H = [R ones(size(R))];
160
161% Smooth data
162if all(smooth_kernel)
163  %fprintf(1,'Smoothing data...');
164  smooth_data = fmri_smooth(data,smooth_kernel);
165  %fprintf(1,'...done\n');
166else
167  fprintf(1,'fMRI analysis warning: Could not smooth data!\n');
168  smooth_data = data;
169end
170
171% Calculate parametric map(s)
172nPlanes = size(smooth_data,3);
173nCols = size(smooth_data,2);
174pmap1 = zeros(size(smooth_data,1),size(smooth_data,2),size(smooth_data,3));
175pmap2 = zeros(size(pmap1));
176tmap = zeros(size(pmap1));
177c = repmat(contr,nCols,1);
178HTH = pinv(H'*H);
179R = eye(size(H,1))-H*HTH*H';
180for ii=1:nPlanes
181  fprintf(1,'Processing plane %d/%d\n',ii,nPlanes);
182  for kk=1:nCols
183    if kk==1
184      fprintf(1,'.....Processing column %03d/%03d',kk,nCols);
185    else
186      fprintf(1,'\b\b\b\b\b\b\b');
187      fprintf(1,'%03d/%03d',kk,nCols);
188    end
189    drawnow
190    col_data = squeeze(smooth_data(:,kk,ii,:)).';
191    col_means = mean(col_data,1);
192    col_data = detrend(col_data,'linear');
193    col_data = col_data+repmat(col_means,k,1);
194    th=H\col_data;
195    pmap1(:,kk,ii)=th(1,:).';
196    pmap2(:,kk,ii)=th(2,:).';
197   
198    r = col_data-H*th;
199    rr=diag(r'*r);
200    rr=rr(:);
201    sig2 = rr./trace(R);
202    sig2 = repmat(sig2,1,2);
203    T = diag(c*th)./sqrt(c.*sig2*HTH*H'*H*HTH*contr');
204    tmap(:,kk,ii)=T.';
205  end
206  fprintf(1,'\n');
207end
208
209maps_out.pmap1 = pmap1;
210maps_out.pmap2 = pmap2;
211maps_out.tmap = tmap;
212
213% Calculate effective degrees of freedom
214%v = (trace(R).^2)/trace(R*R);
215
216
217
218% % Calculate T-maps and contrasts
219% c = contr;
220%
221% % Projection and residual
222% HTH = pinv(H'*H);
223% R = eye(size(H,1))-H*HTH*H';
224% r = Z-H*th;
225% sig2 = (r'*r)./trace(R);
226%
227% % Calculate T-value
228% T = (c*th)/sqrt(c*sig2*HTH*H'*H*HTH*c');
229%
230
231
232
233
Note: See TracBrowser for help on using the repository browser.

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