TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
BVP_feat_extr.m
Go to the documentation of this file.
1 %This file is part of TEAP.
2 %
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.
7 %
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.
12 %
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/>.
15 %
16 %> @file BVP_feat_extr.m
17 %> @brief Computes @b BVP features
18 %> @param BVPsignal: the BVP signal.
19 %> @param varargin: you can choose which features to extract (see featureSelector)
20 %>
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
32 %> 0.15-0.5Hz bands
33 %> - tachogram_LF: Tachogram's low freqneucy spectral content
34 %> <0.08Hz
35 %> - tachogram_MF: Tachogram's medium freqneucy spectral content
36 %> 0.08Hz> and <0.15Hz
37 %> - tachogram_HF: Tachogram's high freqneucy spectral content
38 %> 0.15Hz> and <0.5Hz
39 %> - tachogram_energy_ratio: energu ratio tachogram_MFSP/(tachogram_HSP+tachogram_LFSP)
40 %
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
47 %> returned
48 %
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)
52 
53 % Check inputs and define unknown values
54 narginchk(1, Inf);
55 
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
59  Bulk = [];
60 end
61 
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{:});
68 
69 %If some features are selected
70 if(~isempty(BVP_feats_names))
71 
72  %Compute the results
73 
74  %First compute IBI if needed by the requested features
75  if(any(ismember(featuresNamesIBI,BVP_feats_names)))
76  %Compute IBI
77  BVPSignal = BVP__compute_IBI( BVPSignal );
78  IBI = Signal__get_raw(BVPSignal.IBI);
79  IBI_sp = Signal__get_samprate(BVPSignal.IBI);
80 
81  %Update the Bulk with the new BVP signal
82  if(~isempty(Bulk))
83  Bulk = Bulk_update_signal(Bulk, Signal__get_signame(BVPSignal), BVPSignal);
84  end
85  end
86 
87  %Get the raw signals
88  rawSignal = Signal__get_raw(BVPSignal);
89  samprate = Signal__get_samprate(BVPSignal);
90  %averaged BVP - ralated to blood pressure
91 
92  if any(strcmp('mean_',BVP_feats_names))
93  mean_ = mean(rawSignal);
94  end
95 
96  if any(strcmp('meanIBI',BVP_feats_names)) || any(strcmp('HRV',BVP_feats_names))
97  HRV = std(IBI);
98  meanIBI = mean(IBI);
99  end
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
103  end
104 
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');
109  end
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;
113 
114  else
115 
116  if any(strncmp('sp',BVP_feats_names,2))
117  [P, f] = pwelch(rawSignal, welch_window_size_BVP, [], [], samprate);
118  P=P/sum(P);
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);
125  end
126  end
127  if length(IBI)< welch_window_size_BVP +IBI_sp
128  warning('Signal is too short for this welch window size');
129  end
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;
134  else
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);
151  end
152  end
153 
154  %Write the values to the final vector output
155  BVP_feats = [];
156  for (i = 1:length(BVP_feats_names))
157  eval(['BVP_feats = cat(2, BVP_feats, ' BVP_feats_names{i} ');']);
158  end
159 
160  else %no features selected
161  BVP_feats = [];
162  end
163 end
BVP_feat_extr
function BVP_feat_extr(in BVPSignal, in varargin)