1 %This file is part of TEAP.
3 %TEAP is free software: you can redistribute it and/or modify
4 %it under the terms of the GNU General Public License as published by
5 %the Free Software Foundation, either version 3 of the License, or
6 %(at your option) any later version.
8 %TEAP is distributed in the hope that it will be useful,
9 %but WITHOUT ANY WARRANTY; without even the implied warranty of
10 %MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 %GNU General Public License
for more details.
13 %You should have received a copy of the GNU General Public License
14 %along with TEAP. If not, see <http://www.gnu.org/licenses/>.
17 %> @brief Computes @b BVP features
18 %> @param BVPsignal: the BVP signal.
19 %> @param varargin: you can choose which features to extract (see featureSelector)
21 %> the list of available features is:
22 %> - mean_: averaged BVP - ralated to blood pressure
23 %> - HRV: heart rate variability calculated based on the standard
24 %> deviation IBI differences
25 %> - meanIBI: mean heart rate (beat per minute)
26 %> - MSE: Multi-Scale Entropy at 5 levels 1-5
27 %> - sp0001: spectral power in 0.0-0.1Hz band
28 %> - sp0102: spectral power in 0.1-2.1Hz band
29 %> - sp0203: spectral power in 0.2-3.1Hz band
30 %> - sp0304: spectral power in 0.3-4.1Hz band
31 %> - sp_energyRatio: spectral power ratio between 0.0-0.08Hz and
33 %> - tachogram_LF: Tachogram
's low freqneucy spectral content
35 %> - tachogram_MF: Tachogram's medium freqneucy spectral content
36 %> 0.08Hz> and <0.15Hz
37 %> - tachogram_HF: Tachogram
's high freqneucy spectral content
39 %> - tachogram_energy_ratio: energu ratio tachogram_MFSP/(tachogram_HSP+tachogram_LFSP)
41 %> @retval BVP_feats: list of features values
42 %> @retval BVP_feats_names: names of the computed features (it is good pratice to
43 %> check this vector since the order of requested features
44 %> can be different than the requested one)
45 %> @retval Bulk: if the input to the function is a Bulk than the Bulk is returned
46 %> with the updated BVP signal, including IBI. Otherwise NaN is
49 %> @author Copyright Guillaume Chanel 2013
50 %> @author Copyright Frank Villaro-Dixon, 2014
51 function [BVP_feats, BVP_feats_names, Bulk] = BVP_feat_extr(BVPSignal,varargin)
53 % Check inputs and define unknown values
56 %Make sure we have an BVP signal
57 [BVPSignal, Bulk] = BVP__assert_type(BVPSignal);
58 if(nargout < 3) %No bulk requested -> do not need to keep it
62 % Define full feature list and get features selected by user
63 featuresNames = {'mean_
', 'HRV
', 'meanIBI
', 'MSE
', ...
64 'sp0001
', 'sp0102
', 'sp0203
', 'sp0304
', 'sp_energyRatio
', ...
65 'tachogram_LF
', 'tachogram_MF
', 'tachogram_HF
', 'tachogram_energy_ratio
'};
66 featuresNamesIBI = {'mean_
', 'HRV
', 'meanIBI
', 'tachogram_LF
', 'tachogram_MF
', 'tachogram_HF
', 'tachogram_energy_ratio
'};
67 BVP_feats_names = featuresSelector(featuresNames,varargin{:});
69 %If some features are selected
70 if(~isempty(BVP_feats_names))
74 %First compute IBI if needed by the requested features
75 if(any(ismember(featuresNamesIBI,BVP_feats_names)))
77 BVPSignal = BVP__compute_IBI( BVPSignal );
78 IBI = Signal__get_raw(BVPSignal.IBI);
79 IBI_sp = Signal__get_samprate(BVPSignal.IBI);
81 %Update the Bulk with the new BVP signal
83 Bulk = Bulk_update_signal(Bulk, Signal__get_signame(BVPSignal), BVPSignal);
88 rawSignal = Signal__get_raw(BVPSignal);
89 samprate = Signal__get_samprate(BVPSignal);
90 %averaged BVP - ralated to blood pressure
92 if any(strcmp('mean_
',BVP_feats_names))
93 mean_ = mean(rawSignal);
96 if any(strcmp('meanIBI
',BVP_feats_names)) || any(strcmp('HRV
',BVP_feats_names))
100 if any(strncmp('MSE
',BVP_feats_names,3))
101 %multi-scale entropy for 5 scales on hrv
102 [MSE] = multiScaleEntropy(IBI,5); %FIXME: this is applied on a 4Hz signal -> but probably the features were originally computed on an IBI sequence
105 welch_window_size_BVP = samprate*20;
106 welch_window_size_IBI = IBI_sp*20;
107 if length(rawSignal)< welch_window_size_BVP +samprate
108 warning('Signal is too
short for this welch window size
');
110 if length(rawSignal)< welch_window_size_BVP +1
111 warning('Signal is too
short for this welch window size and the PSD features cannot be calculated
');
112 sp0001 = NaN;sp0102 = NaN;sp0203 = NaN; sp0304 = NaN;sp_energyRatio = NaN;
116 if any(strncmp('sp
',BVP_feats_names,2))
117 [P, f] = pwelch(rawSignal, welch_window_size_BVP, [], [], samprate);
119 %power spectral featyres
120 sp0001 = log(sum(P(f>0.0 & f<=0.1))+eps);
121 sp0102 = log(sum(P(f>0.1 & f<=0.2))+eps);
122 sp0203 = log(sum(P(f>0.2 & f<=0.3))+eps);
123 sp0304 = log(sum(P(f>0.3 & f<=0.4))+eps);
124 sp_energyRatio = log(sum(P(f<0.08))/sum(P(f>0.15 & f<0.5))+eps);
127 if length(IBI)< welch_window_size_BVP +IBI_sp
128 warning('Signal is too
short for this welch window size
');
130 if length(IBI)< welch_window_size_BVP +1
131 warning('Signal is too
short for this welch window size and the PSD features cannot be calculated
');
132 tachogram_LF = NaN;tachogram_MF = NaN;tachogram_HF = NaN;
133 tachogram_energy_ratio = NaN;
135 %tachogram features; psd features on inter beat intervals
136 %R. McCraty, M. Atkinson, W. Tiller, G. Rein, and A. Watkins, "The
137 %effects of emotions on short-term power spectrum analysis of
138 %heart rate variability," The American Journal of Cardiology, vol. 76,
139 %no. 14, pp. 1089 -1093, 1995
140 if any(strcmp('tachogram_LF
',BVP_feats_names)) ...
141 || any(strcmp('tachogram_MF
',BVP_feats_names)) || any(strcmp('tachogram_HF
',BVP_feats_names)) ...
142 || any(strcmp('tachogram_energy_ratio
',BVP_feats_names))
143 [Pt, ft] = pwelch(IBI, welch_window_size_BVP, [], [], IBI_sp);
144 clear tachogram %TODO: delete ? Why is it useful ?
145 %WARN: check that this is possible with the IBI sampling rate
146 %WARN: these values are sometimes negative because of the log, doesn't it appear as strange
for a user ?
147 tachogram_LF = log(sum(Pt(ft>0.01 & ft<=0.08))+eps);
148 tachogram_MF = log(sum(Pt(ft>0.08 & ft<=0.15))+eps);
149 tachogram_HF = log(sum(Pt(ft>0.15 & ft<=0.5))+eps);
150 tachogram_energy_ratio = tachogram_MF/(tachogram_LF+tachogram_HF);
154 %Write the values to the
final vector output
156 for (i = 1:length(BVP_feats_names))
157 eval([
'BVP_feats = cat(2, BVP_feats, ' BVP_feats_names{i}
');']);
160 else %no features selected