-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBP_HRV_Frequency_Domain.m
165 lines (136 loc) · 6.75 KB
/
BP_HRV_Frequency_Domain.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
function [BRSPSD, BRSFFT] = ...
BP_HRV_Frequency_Domain(RpLoc, rrInt, SpLoc, Sp, position, debug, dataNum)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates HRV, SPV power spectra and BRS values in the frequency domain
% using both pwelch and fft
%
% INPUT:
% rrInt - vector of RR in ms
% RpLoc - vector of R peaks locations in ms
% SpLoc - vector of systolic pressure peaks locations in ms
% Sp - vector of systolic pressure values in mmHg
% position - supine or tilt, for plotting
% debug - logical value, 0 no plot / 1 plot
% dataNum - subject ID, for plotting
%
% RETURN VALUES:
% density of RR Tachogram in ms²
% BRSPSD - values using pwelch. Struct with fields:
% BRSPSD.psd_spectrum - vector of HRV power spectral density (psd) in
% ms²/Hz
% BRSPSD.psd_f - vector of frequencies at which psd was calculated
% in Hz
% BRSPSD.HrvLf - spectral power in the LF range in ms²
% BRSPSD.HrvHf - spectral power in the HF range in ms²
% BRSPSD.total_power - total spectral power (VLF + LF + HF) in ms²
% BRSPSD.HrvLf_nu - normalized spectral power in the LF range
% BRSPSD.HrvHf_nu - normalized spectral power in the HF range
% BRSPSD.HrvSv - LF/HF ratio, sympathovagal balance
% BRSFFT - values using pwelch. Struct with same fields as BRSPSD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% frequency range
HRV_veryLF_low = 0.01;
HRV_veryLF_high = 0.04;
HRV_LF_low = 0.04;
HRV_LF_high = 0.15;
HRV_HF_low = 0.15;
HRV_HF_high = 0.4;
debug = 0; % change here for 0 plot/ 1 no plot
f_int = 4 ; % frequency of interpolation in Hz
beatpos = RpLoc - RpLoc(1); % in order to start at time 0
beatpos = beatpos./1000 ; % from ms to s
t_rrInt =beatpos(1):1/f_int:beatpos(end-1);
%resampling at 4 Hz
rrIntRes = interp1(beatpos(1:end-1), rrInt, t_rrInt);
if debug == 1 % Plot of equidistant RR time series (after interpolation)
figure;
stem(beatpos(1:end-1), rrInt, 'r');
hold on;
stem(t_rrInt, rrIntRes, 'g*');
hold off
xlabel("Time (s)");
ylabel("RR-interval (ms)");
title(['HRV ' position ' position No:' num2str(dataNum)])
legend('Original HRV','Interpolated HRV')
end
% calculate HRV power spectral density using psd
% use of a hamming window of length 512 with 50% overlap (207)
[RR_FFT, f_FFT] = pwelch((rrIntRes(1:end-10) - nanmean(rrIntRes)), hamming(512), 207, 1024, f_int);
BRSPSD.psd_spectrum = RR_FFT;
BRSPSD.psd_f = f_FFT;
BRSPSD.HrvLf = sum(RR_FFT((f_FFT <= HRV_LF_high) & (f_FFT >= HRV_LF_low))) * f_int/1024 ;
BRSPSD.HrvHf = sum(RR_FFT((f_FFT <= HRV_HF_high) & (f_FFT >= HRV_HF_low))) * f_int/1024 ;
nu = (sum(RR_FFT((f_FFT <= HRV_HF_high) & (f_FFT >= HRV_veryLF_low))) - sum(RR_FFT((f_FFT <= HRV_veryLF_high) & (f_FFT >= HRV_veryLF_low))) ) * f_int/1024 ;
BRSPSD.total_power = nu;
BRSPSD.HrvLf_nu = BRSPSD.HrvLf / nu;
BRSPSD.HrvHf_nu = BRSPSD.HrvHf / nu;
BRSPSD.HrvSv = BRSPSD.HrvLf_nu / BRSPSD.HrvHf_nu;
% calculate HRV power spectral density using fft
L = length(rrIntRes)-10;
Y = fft(rrIntRes(1:end-10) - nanmean(rrIntRes) );
Y_FFT = abs( Y(1:L/2 + 1)./ L );
Y_FFT(2:end-1) = 2 .* Y_FFT(2:end-1);
Y_FFT = Y_FFT .* Y_FFT ;
F_FFT = (0:L/2) * f_int/L;
BRSFFT.fft_spectrum = Y_FFT;
BRSFFT.fft_f = F_FFT;
BRSFFT.HrvLf = sum(Y_FFT((F_FFT <= HRV_LF_high) & (F_FFT >= HRV_LF_low))) ;
BRSFFT.HrvHf = sum(Y_FFT((F_FFT <= HRV_HF_high) & (F_FFT >= HRV_HF_low))) ;
nu = (sum(Y_FFT((F_FFT <= HRV_HF_high) & (F_FFT >= HRV_veryLF_low))) - sum(Y_FFT((F_FFT <= HRV_veryLF_high) & (F_FFT >= HRV_veryLF_low))) ) ;
BRSFFT.total_power = nu;
BRSFFT.HrvLf_nu = BRSFFT.HrvLf / nu;
BRSFFT.HrvHf_nu = BRSFFT.HrvHf / nu;
BRSFFT.HrvSv = BRSFFT.HrvLf_nu / BRSFFT.HrvHf_nu;
% Ps
pspos = SpLoc - SpLoc(1); % to start at time 0
pspos = pspos./1000 ; % from ms to s
t_ps_int =pspos(1):1/f_int:pspos(end);
%resampling with 4 Hz
ps_int = interp1(pspos, Sp, t_ps_int);
if debug ==1 % Plot of equidistant Sp time series (after interpolation)
figure;
stem(pspos, Sp, 'r');
hold on;
stem(t_ps_int, ps_int, 'g*');
hold off
xlabel("Time (s)");
ylabel("Ps (mmHg)");
title(['Systolic blood pressure ' position ' position' ' No:' num2str(dataNum)])
legend('Original Ps','Interpolated Ps')
end
% calculate SPV power spectral density using psd
% use of a hamming window of length 512 with 50% overlap (207) -> this
% should match HRV
[PS_FFT, fp_FFT] = pwelch(ps_int(1:end-10) - nanmean(ps_int), hamming(512), 207, 1024, f_int);
BRSPSD.PsvLf = sum(PS_FFT((fp_FFT <= HRV_LF_high) & (fp_FFT >= HRV_LF_low))) * f_int/401 ;
BRSPSD.PsvHf = sum(PS_FFT((fp_FFT <= HRV_HF_high) & (fp_FFT >= HRV_HF_low))) * f_int/401 ;
BRSPSD.alphaLf = sqrt(BRSPSD.HrvLf/BRSPSD.PsvLf) ;
BRSPSD.alphaHf = sqrt(BRSPSD.HrvHf/BRSPSD.PsvHf) ;
L = length(ps_int)-10;
X = fft(ps_int(1:end-10) - nanmean(ps_int) );
X_FFT = abs( X(1:L/2 + 1)./ L );
X_FFT(2:end-1) = 2 .* X_FFT(2:end-1);
X_FFT = X_FFT .* X_FFT ;
Fp_FFT = (0:L/2) * f_int/L;
BRSFFT.PsvLf = sum(X_FFT((Fp_FFT <= HRV_LF_high) & (Fp_FFT >= HRV_LF_low))) ;
BRSFFT.PsvHf = sum(X_FFT((Fp_FFT <= HRV_HF_high) & (Fp_FFT >= HRV_HF_low))) ;
BRSFFT.alphaLf = sqrt(BRSFFT.HrvLf/BRSFFT.PsvLf) ;
BRSFFT.alphaHf = sqrt(BRSFFT.HrvHf/BRSFFT.PsvHf) ;
[CS_FFT, fcs_FFT] = cpsd( rrIntRes( 1:min( [length(ps_int), length(rrIntRes)])-10 ) - nanmean(rrIntRes), ...
ps_int( 1:min( [length(ps_int), length(rrIntRes)] )-10 ) - nanmean(ps_int), hamming(512) , 207, 1024, f_int);
CS_FFT = sqrt(abs(CS_FFT));
CsLf = mean( CS_FFT( (fcs_FFT <= HRV_LF_high) & (fcs_FFT >= HRV_LF_low) ) ) ;
CsHf = mean( CS_FFT( (fcs_FFT <= HRV_HF_high) & (fcs_FFT >= HRV_HF_low) ) ) ;
BRSPSD.CsLf = CsLf ;
BRSPSD.CsHf = CsHf ;
BRSPSD.cohereLf = 0;
BRSPSD.cohereHf = 0;
[cohere_FFT,fcohere_FFT] = mscohere( rrIntRes( 1:min( [length(ps_int), length(rrIntRes)])-10 ) - nanmean(rrIntRes), ...
ps_int( 1:min( [length(ps_int), length(rrIntRes)] )-10 ) - nanmean(ps_int), hamming(512), 207, 1024, f_int);
BRSPSD.cohereLf = sum( cohere_FFT( (fcohere_FFT <= HRV_LF_high) & (fcohere_FFT >= HRV_LF_low) ) ...
.* fcohere_FFT( (fcohere_FFT <= HRV_LF_high) & (fcohere_FFT >= HRV_LF_low) ) ) ...
/ sum( fcohere_FFT( (fcohere_FFT <= HRV_LF_high) & (fcohere_FFT >= HRV_LF_low) ) );
BRSPSD.cohereHf = sum( cohere_FFT( (fcohere_FFT <= HRV_HF_high) & (fcohere_FFT >= HRV_HF_low) ) ...
.* fcohere_FFT( (fcohere_FFT <= HRV_HF_high) & (fcohere_FFT >= HRV_HF_low) ) ) ...
/ sum( fcohere_FFT( (fcohere_FFT <= HRV_HF_high) & (fcohere_FFT >= HRV_HF_low) ) );
end