-
Notifications
You must be signed in to change notification settings - Fork 4
/
heart_buildHRVregressors.m
126 lines (95 loc) · 4.61 KB
/
heart_buildHRVregressors.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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
function heart_buildHRVregressors(subj_idx)
% 0.15 to 0.4
% 0.04 to 0.15
% filter
% hilbert
% store envelope
% make fisrt level GLM
cfgMain = global_getcfgmain
% subj_idx = 13
load(['/media/ignacio/IREBOLLOENS/HRVts/IBIts_S' sprintf('%.2d',subj_idx) '.mat'])
%%
%
% HRV_timeseries_filename = global_filename(subj_idx,cfgMain,'HRV_timeseries')
% load(HRV_timeseries_filename)
%
%
% HRV_notdownsampled = ibi_int';
% load(strcat(global_path2root,'\scripts4paper\files\sampleFieldtripStruc.mat'))
load('/home/ignacio/Matlab/scripts4paper/files/sampleFieldtripStruc.mat')
% labelsChannelsMAIN = {'HRV'};
% labelsChannels = labelsChannelsMAIN;
% clusterRegionsComparisons = HRV_notdownsampled;
% dataStructure.hdr = EGG_downsampled.hdr;
% dataStructure.fsample = 1;
% dataStructure.time{1,1} = [0:1/dataStructure.fsample:(size(clusterRegionsComparisons,2))-1];
% dataStructure.label = labelsChannels;%channelStr;
% dataStructure.cfg = EGG_downsampled.cfg;
% dataStructure.trial{1,1} = clusterRegionsComparisons;
% dataStructure.sampleinfo = [1 length(HRV_notdownsampled)]
% disp('Resampling...')
% cfg = []; %initialize configuration structure
% cfg.detrend = 'no'; % remove linear trend from the data (done per trial)
% cfg.demean = 'yes';
% cfg.resamplefs= 0.5; % 4 x top-freq (15 cpm = 0.25 Hz) - Nyquist = 30 cpm frequency at which the data will be resampled
% HRV_downsampled = ft_resampledata(cfg,dataStructure); % This procedure also lowpass filter the data at half the new sr
%
% HRV_polyremoval= ft_preproc_polyremoval (ibi_int, 2);
% HRV_downsampled.trial{1,1} = ft_preproc_standardize (HRV_downsampled. trial{1,1});
%
%% LF HRV
filter_frequency_spread=0.06 ; % In hz Half widfth or fullwidth????
centerFrequency = 0.1
sr = 1; % 1 TR = 2s
filterOrder=(cfgMain.fOrder*fix(sr/(centerFrequency-filter_frequency_spread))-1);%in nsamples
transition_width= cfgMain.transitionWidth/100; % in normalised units
LFfilteredHeart=tools_bpFilter(ibi_int,sr,filterOrder,centerFrequency,filter_frequency_spread,transition_width,cfgMain.filterType);
HilbertLFHRV = hilbert(LFfilteredHeart);
phaseLFHRV = angle(HilbertLFHRV); % Cut data to have the same length as EGG (cut this way to get rid of fmri edge artifact on EGG)
AmplitudeEnvelopeLFHRV = abs(HilbertLFHRV); % Cut data to have the same length as EGG (cut this way to get rid of fmri edge artifact on EGG)
figure
plot(LFfilteredHeart)
hold on
plot(AmplitudeEnvelopeLFHRV,'r')
%% HF HRV
filter_frequency_spread=0.125 ; % In hz Half widfth or fullwidth????
centerFrequency = 0.275
sr = 1; % 1 TR = 2s
filterOrder=(cfgMain.fOrder*fix(sr/(centerFrequency-filter_frequency_spread))-1);%in nsamples
transition_width= cfgMain.transitionWidth/100; % in normalised units
HFfilteredHeart=tools_bpFilter(ibi_int,sr,filterOrder,centerFrequency,filter_frequency_spread,transition_width,cfgMain.filterType);
HilbertHFHRV = hilbert(HFfilteredHeart);
phaseHFHRV = angle(HilbertHFHRV); % Cut data to have the same length as EGG (cut this way to get rid of fmri edge artifact on EGG)
AmplitudeEnvelopeHFHRV = abs(HilbertHFHRV); % Cut data to have the same length as EGG (cut this way to get rid of fmri edge artifact on EGG)
figure
plot(ibi_int)
hold on
plot(AmplitudeEnvelopeLFHRV,'r')
plot(AmplitudeEnvelopeHFHRV,'G')
% downsample lFHRV envelope
% load(strcat(global_path2root,'\scripts4paper\files\sampleFieldtripStruc.mat'))
load('/home/ignacio/Matlab/scripts4paper/files/sampleFieldtripStruc.mat')
labelsChannelsMAIN = {'LFHRV','HFHRV'};
labelsChannels = labelsChannelsMAIN;
clusterRegionsComparisons = [AmplitudeEnvelopeLFHRV';AmplitudeEnvelopeHFHRV'];
dataStructure.hdr = EGG_downsampled.hdr;
dataStructure.fsample = 1;
dataStructure.time{1,1} = [0:1/dataStructure.fsample:(size(clusterRegionsComparisons,2))-1];
dataStructure.label = labelsChannels;%channelStr;
dataStructure.cfg = EGG_downsampled.cfg;
dataStructure.trial{1,1} = clusterRegionsComparisons;
dataStructure.sampleinfo = [1 length(AmplitudeEnvelopeLFHRV)]
disp('Resampling...')
cfg = []; %initialize configuration structure
cfg.detrend = 'no'; % remove linear trend from the data (done per trial)
cfg.demean = 'yes';
cfg.resamplefs= 0.5; % 4 x top-freq (15 cpm = 0.25 Hz) - Nyquist = 30 cpm frequency at which the data will be resampled
HRV_downsampled = ft_resampledata(cfg,dataStructure); % This procedure also lowpass filter the data at half the new sr
figure
% plot(ibi_int)
hold on
plot(HRV_downsampled.trial{1,1}(1,:),'-or')
plot(HRV_downsampled.trial{1,1}(2,:),'-og')
HFHRVregressor = zscore(HRV_downsampled.trial{1,1}(2,:)')
LFHRVregressor = zscore(HRV_downsampled.trial{1,1}(1,:)')
save(['/media/ignacio/IREBOLLOENS/HRVts/HRVregressors_S' sprintf('%.2d',subj_idx) '.mat'])