-
Notifications
You must be signed in to change notification settings - Fork 5
/
ts_corr_basic.m
89 lines (66 loc) · 2.01 KB
/
ts_corr_basic.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
function betas = ts_corr_basic(ts,ons,pmods,varargin)
%
% compute BOLD-regressor correlations (the "Behrens" plot)
%
% ARGS:
% ts - a time series extracted with spm_regions [vector]
% ons - event onsets (in secs) [vector]
% pmods - parametric values (can be a matrix nTrial*nPmods)
%
% options
% 'dur' - how long is the window to be deplayed (in secs) [scalar]
% 'orth' - orthognalize the PMODs (when nPM >= 2)
%
% Earlier version of this script was developed by Gerhard Jocham [[email protected]]
%
%%
% 'defaults' if not given
param.dur = 10;
param.orth = 0;
% record input arges
i = 1;
while i < nargin - 3
arg = varargin{i};
val = varargin{i+1};
param = fillpar(param, arg, val);
i = i+2;
end
if length(ons) ~= length(pmods)
error('onset vector and parametric value vactor have to be of same length and orientation')
end
if param.orth == 1 && size(pmods,2) < 2
error('only set orth to be true when there are at least 2 PMs!')
end
if param.orth == 1
pmods = spm_orth(pmods);
end
%%
upsrate = 10; % up sampling rate
TR = 2.51;
% prepare time series
ts = ts - (mean(ts)); % de-mean
% ts = ts(mean(ts)); % de-mean
% upsample time series
t = 0:length(ts)-1;
t_ups = 0:1/upsrate:length(ts)-1;
ts_ups = spline(t,ts,t_ups);
% create scanindex for data to extract from time series
onset = round(ons * upsrate / TR);
duration = round(param.dur * upsrate / TR);
%scanidx = repmat(onset,1,length(duration)) + repmat(0:duration-1,length(onset),1);
scanidx = repmat(onset,1, duration) + repmat(0:duration-1,length(onset),1);
% extract data from time sereis (sorted by trials (=row) * duration (=column)
ps = ts_ups(scanidx);
% normalize data and pmods
pmodn = normalise(pmods);
psn = normalise(ps);
betas = [pinv(pmodn) * psn]';
function param = fillpar(param, arg, val)
switch lower(arg)
case 'dur'
param.dur = val;
case 'orth'
param.orth = val;
end
end
end