Skip to content

Commit 4fe85f5

Browse files
MSendenMSenden
authored
Merge pull request #4 from ccnmaastricht/matlab
Matlab
2 parents dc2a8d5 + d5f3f20 commit 4fe85f5

File tree

4 files changed

+153
-468
lines changed

4 files changed

+153
-468
lines changed

‎code/matlab/IRM.m

Lines changed: 26 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,8 @@
2525
%
2626
% Input-referred model (IRM) mapping tool.
2727
%
28-
% irm = IRM(params) creates an instance of the IRM class.
29-
% params is a structure with 5 required fields
28+
% irm = IRM(parameters) creates an instance of the IRM class.
29+
% parameters is a structure with 5 required fields
3030
% - f_sampling: sampling frequency (1/TR)
3131
% - n_samples : number of samples (volumes)
3232
% - n_rows : number of rows (in-plane resolution)
@@ -55,7 +55,7 @@
5555
% (e.g. help IRM.mapping)
5656
%
5757
% typical workflow:
58-
% 1. irm = IRM(params);
58+
% 1. irm = IRM(parameters);
5959
% 2. irm.set_stimulus();
6060
% 3. irm.create_timecourse(FUN,xdata);
6161
% 4. results = irm.mapping(data);
@@ -87,24 +87,24 @@
8787

8888
methods (Access = public)
8989

90-
function self = IRM(params,varargin)
90+
function self = IRM(parameters,varargin)
9191
% constructor
9292
p = inputParser;
93-
addRequired(p,'params',@isstruct);
93+
addRequired(p,'parameters',@isstruct);
9494
addOptional(p,'hrf',[]);
95-
p.parse(params,varargin{:});
95+
p.parse(parameters,varargin{:});
9696

9797
self.is = 'input-referred modeling tool';
9898

9999
self.two_gamma = @(t) (6*t.^5.*exp(-t))./gamma(6)...
100100
-1/6*(16*t.^15.*exp(-t))/gamma(16);
101101

102-
self.f_sampling = p.Results.params.f_sampling;
102+
self.f_sampling = p.Results.parameters.f_sampling;
103103
self.p_sampling = 1/self.f_sampling;
104-
self.n_samples = p.Results.params.n_samples;
105-
self.n_rows = p.Results.params.n_rows;
106-
self.n_cols = p.Results.params.n_cols;
107-
self.n_slices = p.Results.params.n_slices;
104+
self.n_samples = p.Results.parameters.n_samples;
105+
self.n_rows = p.Results.parameters.n_rows;
106+
self.n_cols = p.Results.parameters.n_cols;
107+
self.n_slices = p.Results.parameters.n_slices;
108108
self.n_total = self.n_rows*self.n_cols*self.n_slices;
109109

110110
if ~isempty(p.Results.hrf)
@@ -228,7 +228,7 @@ function create_timecourse(self,FUN,xdata)
228228
%
229229
% optional inputs are
230230
% - threshold: minimum voxel intensity (default = 100.0)
231-
% - mask : binary mask for selecting voxels
231+
% - mask : binary mask for selecting voxels
232232

233233
text = 'mapping input-referred model...';
234234
fprintf('%s\n',text)
@@ -247,7 +247,7 @@ function create_timecourse(self,FUN,xdata)
247247
data = reshape(data(1:self.n_samples,:,:,:),...
248248
self.n_samples,self.n_total);
249249
mean_signal = mean(data);
250-
data = zscore(data);
250+
251251

252252
if isempty(mask)
253253
mask = mean_signal>=threshold;
@@ -256,53 +256,51 @@ function create_timecourse(self,FUN,xdata)
256256
voxel_index = find(mask);
257257
n_voxels = numel(voxel_index);
258258

259-
mag_d = sqrt(sum(data(:,mask).^2));
259+
data = zscore(data(:,mask));
260+
mag_d = sqrt(sum(data.^2));
260261

261262
results.R = zeros(self.n_total,1);
262263
results.P = zeros(self.n_total,self.n_predictors);
263264

264265
if size(self.hrf,2)==1
265-
hrf_fft = fft(repmat([self.hrf;...
266-
zeros(self.n_samples-self.l_hrf,1)],...
267-
[1,self.n_points]));
266+
hrf_fft = fft([self.hrf;...
267+
zeros(self.n_samples-self.l_hrf,1)]);
268268
tc = zscore(ifft(self.tc_fft.*hrf_fft))';
269269

270270
mag_tc = sqrt(sum(tc.^2,2));
271271
for m=1:n_voxels
272272
v = voxel_index(m);
273273

274-
CS = (tc*data(:,v))./...
275-
(mag_tc*mag_d(v));
274+
CS = (tc*data(:,m))./...
275+
(mag_tc*mag_d(m));
276276
id = isinf(CS) | isnan(CS);
277277
CS(id) = 0;
278278
[results.R(v),j] = max(CS);
279279
for p=1:self.n_predictors
280280
results.P(v,p) = self.xdata{p}(self.idx(j,p));
281281
end
282282

283-
waitbar(v/self.n_total,wb)
283+
waitbar(m/n_voxels,wb)
284284
end
285285
else
286-
hrf_fft_all = fft([self.hrf;...
287-
zeros(self.n_samples-self.l_hrf,self.n_total)]);
286+
hrf_fft_all = fft([self.hrf(:,mask);...
287+
zeros(self.n_samples-self.l_hrf,n_voxels)]);
288288
for m=1:n_voxels
289289
v = voxel_index(m);
290290

291-
hrf_fft = repmat(hrf_fft_all(:,v),...
292-
[1,self.n_points]);
293-
tc = zscore(ifft(self.tc_fft.*hrf_fft))';
291+
tc = zscore(ifft(self.tc_fft.*hrf_fft_all(:,m)))';
294292
mag_tc = sqrt(sum(tc.^2,2));
295293

296-
CS = (tc*data(:,v))./...
297-
(mag_tc*mag_d(v));
294+
CS = (tc*data(:,m))./...
295+
(mag_tc*mag_d(m));
298296
id = isinf(CS) | isnan(CS);
299297
CS(id) = 0;
300298
[results.R(v),j] = max(CS);
301299
for p=1:self.n_predictors
302300
results.P(v,p) = self.xdata(self.idx(j,p),p);
303301
end
304302

305-
waitbar(v/self.n_total,wb)
303+
waitbar(m/n_voxels,wb)
306304
end
307305
end
308306
results.R = reshape(results.R,self.n_rows,self.n_cols,self.n_slices);

‎code/matlab/PEA.m

Lines changed: 68 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,8 @@
2525
%
2626
% Phase-encoding analysis tool.
2727
%
28-
% pea = PEA(params) creates an instance of the PEA class.
29-
% params is a structure with 7 required fields
28+
% pea = PEA(parameters) creates an instance of the PEA class.
29+
% parameters is a structure with 7 required fields
3030
% - f_sampling: sampling frequency (1/TR)
3131
% - f_stim : stimulation frequency
3232
% - n_samples : number of samples (volumes)
@@ -46,7 +46,7 @@
4646
% function (e.g. help PEA.fitting)
4747
%
4848
% typical workflow:
49-
% 1. pea = PEA(params);
49+
% 1. pea = PEA(parameters);
5050
% 2. pea.set_delay(delay);
5151
% 3. pea.set_direction(direction);
5252
% 4. results = pea.fitting(data);
@@ -74,17 +74,17 @@
7474

7575
methods (Access = public)
7676

77-
function self = PEA(params)
77+
function self = PEA(parameters)
7878
% constructor
7979
self.is = 'PEA tool';
8080

81-
self.f_sampling = params.f_sampling;
82-
self.f_stim = params.f_stim;
81+
self.f_sampling = parameters.f_sampling;
82+
self.f_stim = parameters.f_stim;
8383
self.p_sampling = 1/self.f_sampling;
84-
self.n_samples = params.n_samples;
85-
self.n_rows = params.n_rows;
86-
self.n_cols = params.n_cols;
87-
self.n_slices = params.n_slices;
84+
self.n_samples = parameters.n_samples;
85+
self.n_rows = parameters.n_rows;
86+
self.n_cols = parameters.n_cols;
87+
self.n_slices = parameters.n_slices;
8888
self.n_total = self.n_rows*self.n_cols*self.n_slices;
8989
self.t = (0:self.p_sampling:self.p_sampling*self.n_samples-1);
9090
self.delay = 0;
@@ -132,7 +132,7 @@ function set_direction(self,direction)
132132
end
133133
end
134134

135-
function results = fitting(self,data)
135+
function results = fitting(self,data,varargin)
136136
% identifies phase and amplitude at stimulation frequency for
137137
% each voxel and returns a structure with the following fields
138138
% - Phase
@@ -146,50 +146,77 @@ function set_direction(self,direction)
146146
% required inputs are
147147
% - data : a matrix of empirically observed BOLD timecourses
148148
% whose columns correspond to time (volumes).
149+
%
150+
% optional inputs are
151+
% - threshold: minimum voxel intensity (default = 100.0)
152+
% - mask : binary mask for selecting voxels
149153

150154
text = 'performing phase encoding analysis...';
151155
fprintf('%s\n',text)
152156
wb = waitbar(0,text,'Name',self.is);
153157

158+
p = inputParser;
159+
addRequired(p,'data',@isnumeric);
160+
addOptional(p,'threshold',100);
161+
addOptional(p,'mask',[]);
162+
p.parse(data,varargin{:});
163+
164+
data = single(p.Results.data);
165+
threshold = p.Results.threshold;
166+
mask = p.Results.mask;
167+
154168
F = exp(self.direction*2i*pi*self.f_stim*(self.t-self.delay));
155169
X = zscore([real(F)',imag(F)']);
156-
XX = (X'*X)\X';
157-
data = zscore(reshape(single(data(1:self.n_samples,:,:,:)),...
158-
self.n_samples,self.n_total));
159-
std_signal = std(data);
170+
171+
data = reshape(single(data(1:self.n_samples,:,:,:)),...
172+
self.n_samples,self.n_total);
173+
174+
mean_signal = mean(data);
175+
176+
if isempty(mask)
177+
mask = mean_signal>=threshold;
178+
end
179+
mask = mask(:);
180+
voxel_index = find(mask);
181+
n_voxels = numel(voxel_index);
182+
183+
data = zscore(data(:, mask));
184+
185+
beta = (X' * X) \ X' * data;
186+
187+
Y_ = X * beta;
188+
residuals = data - Y_;
189+
160190
results.Phase = zeros(self.n_total,1);
161191
results.Amplitude = zeros(self.n_total,1);
162192
results.F_stat = zeros(self.n_total,1);
163193
results.P_value = ones(self.n_total,1);
164194

165195
df1 = 1;
166196
df2 = self.n_samples-2;
167-
for v=1:self.n_total
168-
if std_signal>0
169-
b = XX * data(:,v);
170-
y = X*b;
171-
172-
% estimate and correct for autocorrelation
173-
r = data(:,v) - y;
174-
T = [[0;r(1:end-1)],[zeros(2,1);r(1:end-2)]];
175-
W = [1; -((T'* T) \ T' * r)];
176-
Xc(:,1) = [X(:,1),[0;X(1:end-1,1)],[zeros(2,1);X(1:end-2,1)]] * W;
177-
Xc(:,2) = [X(:,2),[0;X(1:end-1,2)],[zeros(2,1);X(1:end-2,2)]] * W;
178-
Dc = [data(:,v),[0;data(1:end-1,v)],[zeros(2,1);data(1:end-2,v)]] * W;
179-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180-
181-
b = (Xc' * Xc) \ Xc' * Dc;
182-
y = Xc * b;
183-
mu = mean(DC);
184-
MSM = (y-mu)'*(y-mu)/df1;
185-
MSE = (y-Dc)'*(y-Dc)/df2;
186-
187-
results.Phase(v) = angle(b(1)+b(2)*1i);
188-
results.Amplitude(v) = abs(b(1)+b(2)*1i);
189-
results.F_stat(v) = MSM/MSE;
190-
results.P_value(v) = max(1-fcdf(MSM/MSE,df1,df2),1e-20);
191-
end
192-
waitbar(v/self.n_total,wb)
197+
for m=1:n_voxels
198+
v = voxel_index(m);
199+
200+
% estimate and correct for autocorrelation
201+
T = [[0;residuals(1:end-1,m)],[zeros(2,1); residuals(1:end-2,m)]];
202+
W = [1; -((T'* T) \ T' * residuals(:,m))];
203+
Xc(:,1) = [X(:,1), [0;X(1:end-1,1)], [zeros(2,1);X(1:end-2,1)]] * W;
204+
Xc(:,2) = [X(:,2), [0;X(1:end-1,2)], [zeros(2,1);X(1:end-2,2)]] * W;
205+
Dc = [data(:,m), [0;data(1:end-1,m)], [zeros(2,1);data(1:end-2,m)]] * W;
206+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
207+
208+
b = (Xc' * Xc) \ Xc' * Dc;
209+
y = Xc * b;
210+
mu = mean(DC);
211+
MSM = (y-mu)'*(y-mu)/df1;
212+
MSE = (y-Dc)'*(y-Dc)/df2;
213+
214+
results.Phase(v) = angle(b(1)+b(2)*1i);
215+
results.Amplitude(v) = abs(b(1)+b(2)*1i);
216+
results.F_stat(v) = MSM/MSE;
217+
results.P_value(v) = max(1-fcdf(MSM/MSE,df1,df2),1e-20);
218+
219+
waitbar(m/n_voxels,wb)
193220
end
194221

195222
results.Phase = reshape(results.Phase,...

0 commit comments

Comments
 (0)