# Stroke Volume Estimation from Respiratory Inductive Plethysmography: Double Empirical Decomposition

Authors: Enas Abdulhay; Pierre-Yves Gumery; R. Ilango; S. Hariharasitaraman; M. Thilagaraj; N. Arunkumar; Gustavo Ramirez-Gonzalez
Journal: Scientific Programming
 (2022)

---

## Abstract

In this study, we have developed a “double-empirical mode decomposition algorithm” to estimate cardiac stroke volume from respiratory inductive plethysmography (RIP) signals. The algorithm consists of first an ensemble empirical mode decomposition (EEMD) to extract the cardiorespiratory components. Then, it is followed by an empirical mode decomposition (EMD) to extract only the cardiac components. This double approach permits (a) solving problems of mixing between cardiac and respiratory components (mode and scale mixing), (b) cardiogenic oscillations extraction in the respiratory inductive plethysmography signal, and (c) subsequent estimation of stroke volume. The algorithm is applied to simulated and real RIP signals. The simulated signals are generated by a cardiorespiratory model previously published by the authors. The real signals are measured via a developed inductive vest. In the real case, the values of estimated stroke volumes are compared to the values obtained by thoracocardiographic filter-based method. In the simulated case, the values are compared to the simulated cardiac activity. The results of comparison through Bland and Altman indicate an error lying in the range ±10%. In contrast to thoracocardiography, the proposed method consists of a promising tool for continuous noninvasive adaptive cardiac monitoring that does not need adjusting parameters or cut-off based on ECG. Also, in comparison to echocardiography and impedance-based methods, it does not necessitate the presence of an expert and is not too sensitive to current penetration.


---

## Body

## 1. Introduction

### 1.1. Background

The mechanical cardiac activity is a vital sign that is essentially monitored in several health cases. However, several recent methods for continuous quantification of ventricular volume or cardiac output include intravascular catheters, contact with radioactivity, availability of practiced expert, or necessity of continuous holding of sensing device. Those requirements make the measurement protocol less effective for continuous monitoring. Therefore, noninvasive procedures have been implemented, e.g., transthoracic electrical bioimpedance and echocardiography. Nonetheless, those techniques are not commonly valid for routine use for the reason that their efficacy in clinical situations and in continuous monitoring is still under examination [1]. On the other hand, a noninvasive technique, named thoracocardiography (TCG), has been introduced in [2]. It is centered on the concept of inductive plethysmography (IP). The outcome of inductive plethysmography is the sum of all variations in the volume surrounded by the device/sensor. According to TCG during spontaneous respiration, movements due to breathing induce about 95% of the signal recorded at the placement of the xiphoid process, while the movement due to left heart ventricle, i.e., cardiogenic oscillations accounts for 5% [2]. Hence, TCG targets the noninvasive recording of left ventricular stroke volume by ECG-triggered ensemble averaging and band-pass filtering [0.7 ∗ (heart rate) Hz–10 Hz] [3] of the resulted signal with the purpose of elimination of low frequency, associated with breathing and body artifacts, as well as of high-frequency noise. TCG leads to accurate quantification of stroke volume and cardiac output. Nevertheless, the cut-off values of the used filter are highly dependent on heart rate, thus TCG cannot be suitable in strong nonstationary conditions.Several candidate methods have been scanned in literature in order to search for a suitable alternative adaptive continuous noninvasive cardiac event detection technique. In [4], a residual deep neural network is customized and employed to detect passive seismic events. The method yields accurate detection. However, the denoising step is achieved by an IIR Wiener filter, i.e., not adaptively by the method itself. In [5], a proposed self-training technique is applied to sound event detection and takes on the task of identifying the presence of specific sound events in a complex audio recording. Nevertheless, it is particularly effective in extending the labeled database with concurrent sound events in case of a lack of fully described class labels. In [6], temporal event information was detected by utilizing features at different scales. Although the approach is promising, its achieved accuracy is not satisfactory for the clinical field of stroke volume estimation. A detection method in [7] is elaborated specifically for noninvasive stroke volume estimation via a novel hierarchical genetic fuzzy logic model. The outcomes are encouraging, but they rather indicate a dependence of stroke volume on the variations in skin warming, age, and pulse pressure of the patient; hence, it is not an automatic technique of stroke volume quantification. In [8], another technique has been suggested particularly for noninvasive stroke volume estimation via indirect myocardial work evaluation. Nonetheless, the basic concept is based on echocardiography that is not appropriate for continuous assessment. In [9], empirical mode decomposition has been applied to detect the variations due to heart beating reflected in breathing waveforms with the aim of tracking heart rate via a piezoelectric transducer. Empirical mode decomposition (EMD) is a nonlinear method that has been developed by Huang et al. [10–12] to decompose nonstationary waveforms into amplitude modulated frequency modulated functions. EMD does not necessitate stationarity/linearity to find out all the embedded oscillatory components [13, 14]. Empirical mode decomposition decomposes complex signals into intrinsic mode functions (IMF), which are high-frequency and low-frequency constituents. [10] illustrates that EMD outperforms wavelet decomposition as it is data driven and does not need a predefined mother model [12, 15]. Furthermore, EMD of biomedical signals yields IMFs based on the real underlying physiological processes. For example, [16] has illustrated the efficacy of EMD in analysis of R-R interval variation caused by breathing—which is the result of nonlinear processes—while other techniques had many drawbacks.The present work introduces a new data-driven technique for the detection of cardiogenic oscillations and stroke volume estimation in IP without the necessity of ECG and in a continuous manner. The suggested technique is formulated in order to solve the limitations mentioned above: it is adapted to nonlinearity and nonstationarity, does not use an ECG reference, fully adaptive, and does not necessitate user-dependent adjusting of parameters (as in TCG), noninvasive in contrast to many methods used in clinical field, appropriate for continuous monitoring with low cost and minimal need of expert presence (in contrast to transthoracic electrical bioimpedance and echocardiography). The suggested adaptive method is a double empirical mode decomposition applied to simulated as well as recorded cardiorespiratory IP waveforms. The signals are hence double-decomposed: (1) in a first step, signals are decomposed by ensemble empirical mode decomposition (explained in the following sections). Then, (2) in a second step, selected resulted components are redecomposed via empirical mode decomposition (explained in the following sections) in order to separate the extracted mechanical cardiac activity. The following subsections explain in details the information in literature related to empirical mode decomposition, confronted limitations (mode mixing and scale mixing), ensemble empirical mode decomposition, and double decomposition.


### 1.2. Empirical Mode Decomposition

A waveform can be named IMF if it satisfies the following requirements: the counts of extrema and zero crossings are equal or differ by one. Also, the local mean of the curves passing by the maxima and minima is about zero. These features of IMF help subsequent accurate calculation of local frequency and amplitude [17].First, maxima and minima are detected in order to interpolate two different curves between the points. Then, the local mean of those curves is calculated and subtracted from the original signal. The remainder is IMF1 if it satisfies the mentioned conditions, and IMF2 is extracted from (original signal–IMF1). Otherwise, the procedure is re-iterated on the remainder to find IMF1. All other IMFs are found out in analogy to the above steps (i.e., sifting) [18]. Stopping criteria are the maximum number of iterations and lowest acceptable local mean [11].The original signals(t) is decomposed as follows:(1)st=∑i=1MIMFit+rt.M: number of IMF, r(t): monotonic residue, and IMFt is the IMF [18].


### 1.3. Mode Mixing

Intermittent waveform is an oscillation—at a specific frequency—appearing or vanishing from a signal abruptly. Hence, in EMD, intermittency takes place when the frequency in an IMF jumps due to the appearance or disappearance of the intermittent section, which leads to an IMF with different frequencies, i.e., mode mixing. Additionally, it has been illustrated in literature that mode mixing is not caused only by mode intermittency but also by EMD dependence on frequency/amplitude ratio of oscillations in a signal. This matter of dependence has been discussed in [19]. The efficacy of separation by EMD decreases when the frequency ratio of two components is lower than a threshold that is a function of amplitude ratio. This is highly noticed when one of the two components has higher frequency and lower amplitude; it is important to emphasize that this is the condition of breathing and cardiac components in a respiratory IP indicated by Figures 1(a) and 1(b). The present work targets the processing of the cardiorespiratory signal in order to move it from the zone of mode mixing to the zone of no mode mixing.Figure 1
(a) Region of work in the case of cardiorespiratory signal obtained by RIP (modified image [19]). The dotted lines define the region of physiological variations of amplitude (a) and frequency (f) ratio. The bad separation of modes exists in zones 2 and 3. Our suggested double decomposition represents passing from A to B. (b) The target of the proposed algorithm is to eliminate mode mixing by changing the amplitude and frequency ratios of respiratory and cardiac activities from the mode mixing situation to the no-mode mixing situation.
(a)(b)


### 1.4. Ensemble Empirical Mode Decomposition

To avoid mode mixing, a new version of EMD has been implemented in [20]: Ensemble empirical mode decomposition (EEMD) method. EEMD inserts a random signal in order to prevent mode mixing. Hence, IMF is actually found out as the average of different counterpart IMFs extracted in different decompositions; every decomposition is the EMD of (original signal + white noise), so that the different white noises annul each other. Consequently, only the real constituents remain.


### 1.5. Double Decomposition

Scale mixing is a drawback of EEMD [20]. Scale mixing occurs when resulted components do not satisfy the strict conditions of IMF due to a ceiling of sifting iterations. In scale mixing, an achieved mode is the sum of two successive IMFs. Hence, a possible solution suggested in the present work is to apply a new EMD sifting to the mode achieved by EEMD.


## 1.1. Background

The mechanical cardiac activity is a vital sign that is essentially monitored in several health cases. However, several recent methods for continuous quantification of ventricular volume or cardiac output include intravascular catheters, contact with radioactivity, availability of practiced expert, or necessity of continuous holding of sensing device. Those requirements make the measurement protocol less effective for continuous monitoring. Therefore, noninvasive procedures have been implemented, e.g., transthoracic electrical bioimpedance and echocardiography. Nonetheless, those techniques are not commonly valid for routine use for the reason that their efficacy in clinical situations and in continuous monitoring is still under examination [1]. On the other hand, a noninvasive technique, named thoracocardiography (TCG), has been introduced in [2]. It is centered on the concept of inductive plethysmography (IP). The outcome of inductive plethysmography is the sum of all variations in the volume surrounded by the device/sensor. According to TCG during spontaneous respiration, movements due to breathing induce about 95% of the signal recorded at the placement of the xiphoid process, while the movement due to left heart ventricle, i.e., cardiogenic oscillations accounts for 5% [2]. Hence, TCG targets the noninvasive recording of left ventricular stroke volume by ECG-triggered ensemble averaging and band-pass filtering [0.7 ∗ (heart rate) Hz–10 Hz] [3] of the resulted signal with the purpose of elimination of low frequency, associated with breathing and body artifacts, as well as of high-frequency noise. TCG leads to accurate quantification of stroke volume and cardiac output. Nevertheless, the cut-off values of the used filter are highly dependent on heart rate, thus TCG cannot be suitable in strong nonstationary conditions.Several candidate methods have been scanned in literature in order to search for a suitable alternative adaptive continuous noninvasive cardiac event detection technique. In [4], a residual deep neural network is customized and employed to detect passive seismic events. The method yields accurate detection. However, the denoising step is achieved by an IIR Wiener filter, i.e., not adaptively by the method itself. In [5], a proposed self-training technique is applied to sound event detection and takes on the task of identifying the presence of specific sound events in a complex audio recording. Nevertheless, it is particularly effective in extending the labeled database with concurrent sound events in case of a lack of fully described class labels. In [6], temporal event information was detected by utilizing features at different scales. Although the approach is promising, its achieved accuracy is not satisfactory for the clinical field of stroke volume estimation. A detection method in [7] is elaborated specifically for noninvasive stroke volume estimation via a novel hierarchical genetic fuzzy logic model. The outcomes are encouraging, but they rather indicate a dependence of stroke volume on the variations in skin warming, age, and pulse pressure of the patient; hence, it is not an automatic technique of stroke volume quantification. In [8], another technique has been suggested particularly for noninvasive stroke volume estimation via indirect myocardial work evaluation. Nonetheless, the basic concept is based on echocardiography that is not appropriate for continuous assessment. In [9], empirical mode decomposition has been applied to detect the variations due to heart beating reflected in breathing waveforms with the aim of tracking heart rate via a piezoelectric transducer. Empirical mode decomposition (EMD) is a nonlinear method that has been developed by Huang et al. [10–12] to decompose nonstationary waveforms into amplitude modulated frequency modulated functions. EMD does not necessitate stationarity/linearity to find out all the embedded oscillatory components [13, 14]. Empirical mode decomposition decomposes complex signals into intrinsic mode functions (IMF), which are high-frequency and low-frequency constituents. [10] illustrates that EMD outperforms wavelet decomposition as it is data driven and does not need a predefined mother model [12, 15]. Furthermore, EMD of biomedical signals yields IMFs based on the real underlying physiological processes. For example, [16] has illustrated the efficacy of EMD in analysis of R-R interval variation caused by breathing—which is the result of nonlinear processes—while other techniques had many drawbacks.The present work introduces a new data-driven technique for the detection of cardiogenic oscillations and stroke volume estimation in IP without the necessity of ECG and in a continuous manner. The suggested technique is formulated in order to solve the limitations mentioned above: it is adapted to nonlinearity and nonstationarity, does not use an ECG reference, fully adaptive, and does not necessitate user-dependent adjusting of parameters (as in TCG), noninvasive in contrast to many methods used in clinical field, appropriate for continuous monitoring with low cost and minimal need of expert presence (in contrast to transthoracic electrical bioimpedance and echocardiography). The suggested adaptive method is a double empirical mode decomposition applied to simulated as well as recorded cardiorespiratory IP waveforms. The signals are hence double-decomposed: (1) in a first step, signals are decomposed by ensemble empirical mode decomposition (explained in the following sections). Then, (2) in a second step, selected resulted components are redecomposed via empirical mode decomposition (explained in the following sections) in order to separate the extracted mechanical cardiac activity. The following subsections explain in details the information in literature related to empirical mode decomposition, confronted limitations (mode mixing and scale mixing), ensemble empirical mode decomposition, and double decomposition.


## 1.2. Empirical Mode Decomposition

A waveform can be named IMF if it satisfies the following requirements: the counts of extrema and zero crossings are equal or differ by one. Also, the local mean of the curves passing by the maxima and minima is about zero. These features of IMF help subsequent accurate calculation of local frequency and amplitude [17].First, maxima and minima are detected in order to interpolate two different curves between the points. Then, the local mean of those curves is calculated and subtracted from the original signal. The remainder is IMF1 if it satisfies the mentioned conditions, and IMF2 is extracted from (original signal–IMF1). Otherwise, the procedure is re-iterated on the remainder to find IMF1. All other IMFs are found out in analogy to the above steps (i.e., sifting) [18]. Stopping criteria are the maximum number of iterations and lowest acceptable local mean [11].The original signals(t) is decomposed as follows:(1)st=∑i=1MIMFit+rt.M: number of IMF, r(t): monotonic residue, and IMFt is the IMF [18].


## 1.3. Mode Mixing

Intermittent waveform is an oscillation—at a specific frequency—appearing or vanishing from a signal abruptly. Hence, in EMD, intermittency takes place when the frequency in an IMF jumps due to the appearance or disappearance of the intermittent section, which leads to an IMF with different frequencies, i.e., mode mixing. Additionally, it has been illustrated in literature that mode mixing is not caused only by mode intermittency but also by EMD dependence on frequency/amplitude ratio of oscillations in a signal. This matter of dependence has been discussed in [19]. The efficacy of separation by EMD decreases when the frequency ratio of two components is lower than a threshold that is a function of amplitude ratio. This is highly noticed when one of the two components has higher frequency and lower amplitude; it is important to emphasize that this is the condition of breathing and cardiac components in a respiratory IP indicated by Figures 1(a) and 1(b). The present work targets the processing of the cardiorespiratory signal in order to move it from the zone of mode mixing to the zone of no mode mixing.Figure 1
(a) Region of work in the case of cardiorespiratory signal obtained by RIP (modified image [19]). The dotted lines define the region of physiological variations of amplitude (a) and frequency (f) ratio. The bad separation of modes exists in zones 2 and 3. Our suggested double decomposition represents passing from A to B. (b) The target of the proposed algorithm is to eliminate mode mixing by changing the amplitude and frequency ratios of respiratory and cardiac activities from the mode mixing situation to the no-mode mixing situation.
(a)(b)


## 1.4. Ensemble Empirical Mode Decomposition

To avoid mode mixing, a new version of EMD has been implemented in [20]: Ensemble empirical mode decomposition (EEMD) method. EEMD inserts a random signal in order to prevent mode mixing. Hence, IMF is actually found out as the average of different counterpart IMFs extracted in different decompositions; every decomposition is the EMD of (original signal + white noise), so that the different white noises annul each other. Consequently, only the real constituents remain.


## 1.5. Double Decomposition

Scale mixing is a drawback of EEMD [20]. Scale mixing occurs when resulted components do not satisfy the strict conditions of IMF due to a ceiling of sifting iterations. In scale mixing, an achieved mode is the sum of two successive IMFs. Hence, a possible solution suggested in the present work is to apply a new EMD sifting to the mode achieved by EEMD.


## 2. Materials and Methods

### 2.1. Subjects and Protocol

Totally, 25 healthy subjects—seated—contributed to the recording phase. Each subject signed the informed consent document. The recording is approved by the committee of the University Hospital Centre of Grenoble-Alpes. The protocol is to breathe spontaneously and quietly during 15 min.


### 2.2. Recording and Calculation of Volumes

The system Visuresp® developed by RBI, Meylan, France, is used to measure the variations of cross-section area of thoracic and abdominal sections through Respiratory IP vest. The vest is fabricated by an elastic material with a particular characteristic of only horizontal extensibility. The inductive spirals are incorporated in the vest at thoracic and abdominal levels with two additional straps to avoid vertical movement. The vest is connected to two primary and two secondary oscillators; the signals pass then to two frequency/voltage converters, two filters, and an Analog to digital converter (PowerLab, ADInstruments). The system is connected by an RS232 cable to a computer equipped with an acquisition software (Chart, ADInstruments). A system of isolation is also included. The zero calibration is automatically performed. ECG and pneumo-tachogram (Fleisch no.1) are also monitored.The calibration of the RIP volume measured by the system is conducted based on the method developed by laboratory team PRETA-France in [21]. At the beginning of recording, the volunteer is asked to perform a number of successive high-amplitude respiratory cycles in order for Visuresp® to determine the range of amplitude to be measured. At the end of recording, RIP signal is considered as ƛ ∗ Abdominal activity + µ ∗ Thoracic activity. The activity means the discrete signal representing the cross section variation measured by the vest. The optimization of ƛ and µ is achieved by minimizing the error between the calculated RIP and the integral value of the airflow measured by pneumotachography.


### 2.3. Cardiorespiratory Model

In the present work, the cardiorespiratory model developed previously by our team in [22–24] is used in order to simulate the variations of respiratory, cardiorespiratory, and cardiac volumes. Briefly, the model includes two modules: respiratory (breathing pattern, alveolar volume, and pleural pressure) and cardiac (cardiac activity, chest wall mechanics, and volume variations (i.e., simulated IP)). The cardiac module is a wave generator. The model takes also into consideration modulation of left ventricle (LV) stroke volume due to breathing.


### 2.4. Parameter Identification and Sensitivity Analysis

Parameter identification and qualitative sensitivity analysis have been conducted using Berkeley Madonna® software in order to determine the most sensitive parameters. Those parameters are varied with the purpose of generating a set of physiological simulations to which our algorithm will be applied.


### 2.5. Simulations

According to sensitivity analysis, the most sensitive parameters are cardiac rate (HR) and breathing frequency (BF), which lead to frequency ratio (f). Also, amplitude ratio (a) between respiratory and cardiac activities is important. In our work, the variations of these parameters from a simulation to another generated a set of seven stationary and nonstationary amplitude modulated simulated RIP signals (more than 1000 cycles). It is worthy to note that these parameters are related to RIP components frequency and amplitude ratio that affects RIP empirical decomposition behavior. The simulations ratios are presented in Table 1.Table 1
Simulated situations. HR = heart rate, Fr: breathing frequency,f: frequency ratio, a: amplitude ratio.
Situation1/f1/aSituation 1 (stationary)4 (HR = 1; fr = 0.25)[0.1–0.17]Situation 2 (nonstationary)f 1 = 7.1; f 2 = 7.8 (HR1 = 1, HR2 = 2.2; Fr1 = 0.14, Fr2 = 0.28)[0.05–0.15]Situation 3 (non–stationary)f 1 = 3.6; f 2 = 7.2 (HR1 = 1, HR2 = 2; Fr1 = 0.28, Fr2 = 0.28)[0.04–0.16]Situation 4 (stationary)3 (HR = 0.75; Fr = 0.25)[0.1–0.15]Situation 5 (stationary)7 (HR = 1; Fr = 0.14)[0.05–0.15]Situation 6 (nonstationary)f 1 = 7.1; f 2 = 7.8 (HR1 = 1, HR2 = 1.8; Fr1 = 0.14, Fr2 = 0.23)[0.05–0.15]Situation 7 (nonstationary)(Progressive rise of HR: [1.4–3.3]; progressive rise of fr: [0.2–0.5])[0.05–0.08]Those simulations are created based on physiological scenarios where basic amplitude ratio (tidal volume/stroke volume) lies approximately in the range [5–20], and basic frequency ratio (breathing frequency/heart rate) lies approximately in the range [1/3–1/8]. For example, superficial high frequency breathing leads to rise in f and fall of a; rapid and profound breathing leads to rise in f and a; drug intoxication leads to a fall in f and a; and emotional stress leads to a rise in a and a fall in f. Note that the values presented in Table 1 are 1/f and 1/a. Moreover, the simulations include the transit from one physiological case to another, i.e., nonstationary situations.


### 2.6. Ensemble Empirical Mode Decomposition

Simulated RIP and real recorded inductive plethysmography signal Vth are decomposed by EEMD. The number of sifting iterations is limited to 10 in order to avoid over-sifting [20]. A set of 10000 white noise signals with amplitude of 1.6 times the r.m.s of RIP signal has been used. The r.m.s VR of RIP signal has been approximated to the r.m.s of the respiratory component:(2)VR≅1M∑k=1MSk2.S(k) is thek th sample of RIP signal. The r.m.s (VR) has been calculated over the number of all samples (M). The signal/noise ratio has been determined as follows:(3)σbVR=1.6,where σb is the added white noise standard deviation.The average of the corresponding IMFs from several tests is then calculated. In the final result, a signal/noise (εr ) ratio was acceptable if near 10%.(4)σrVc=0,1=εr,where σr is the residual noise standard deviation and Vc is the r.m.s of the cardiac component estimated using the amplitude ratio a. We can consider, as indicated in [20] that(5)σr=σbN.Hence,(6)N=σb∗aVR∗εr=1.6∗aεr.The physiological range ofa is between 5 and 20 (e.g., if a preliminary value of a = 6 is chosen, N is then equal to 104).


### 2.7. Double Decomposition and Cardiac Component Reconstruction

As illustrated in the work flowchart in Figure2, a new round of sifting has been applied to modes containing CR mixing.Figure 2
Work flowchart (double decomposition strategy). Vh is the cardiac volume.As mentioned, the algorithm generates a mixture of EEMD scales on certain IMF. This phenomenon is related to the criterion of stopping the sifting operation that takes place on a static number of iterations. This creates a loss of resolution in the separation of different components. This mixing has been detected automatically by Teager–Kaiser energy operator in order to detect the variations in energy present in a signal containing a mixing of two components. The calculated value at an indexk is the difference between the squared amplitude of the component at that k and the product of the amplitude at k − 1 with the amplitude at k + 1. When the Teager–Kaiser value exceeds the threshold 0.00005, a mixing is identified.It is noteworthy that several simulations have been conducted in the present work and all of them gave the same type of results with a scale mixing in IMF4 and IMF5. The sum of those two IMF gives a cardiac signal containing certain breathing harmonics. The solution proposed here is the double decomposition. Figure1(a) illustrated the interest thereof. While the original signal was located in zone 3, the signal of the sum of IMF4 and IMF5 is in Zone 1. The EMD can then separate effectively cardiac and respiratory components [25, 26]. The main parameters used in EMD application are type of extrema detection, type of envelope interpolation, maximum number of sifting iterations, maximum accepted value for the difference between upper and lower envelopes in an IMF, maximum number of obtained IMF, and masking-assisted decomposition. The values of parameters used in EMD were chosen as follows: extrema detection algorithm used in [19], Spline interpolation, 2000, zero, unlimited, and not used, respectively.


### 2.8. Stroke Volume Comparison

For simulated data, left ventricle volume signal (VIv) is used as the reference cardiac activity. Stroke volumes of extracted cardiac activity are hence compared with stroke volumes (SV) of simulated cardiac activity. The used statistical test is Bland and Altman. On the other hand, for real stationary data, cardiac component extracted by a thoracocardiography filter (FIR with order 200) has been used as the reference. More than 170 cycles were used for comparison. The zone in ± R-R interval around the inspiration-expiration transition has been excluded because it suffers from CR harmonics overlapping.The method of evaluation of concordance by Bland and Altman consists of the calculation of the difference between the values obtained by two different techniques as well as the average of those values. When the concordance between two series of values is perfect, the summation of arithmetic differences between the values is null. Also, the average of differences is defined as the bias. The variability of those differences indicates whether the two methods have different natures. In brief, Bland and Altman helps efficiently in judging the interchangeability between the two methods. In the present work, two thresholds of concordance are fixed based on the 95% interval of confidence (±30%). The differences are expressed as percentages of values averages. Bland and Altman graphs have been plotted via MedCalc®.


## 2.1. Subjects and Protocol

Totally, 25 healthy subjects—seated—contributed to the recording phase. Each subject signed the informed consent document. The recording is approved by the committee of the University Hospital Centre of Grenoble-Alpes. The protocol is to breathe spontaneously and quietly during 15 min.


## 2.2. Recording and Calculation of Volumes

The system Visuresp® developed by RBI, Meylan, France, is used to measure the variations of cross-section area of thoracic and abdominal sections through Respiratory IP vest. The vest is fabricated by an elastic material with a particular characteristic of only horizontal extensibility. The inductive spirals are incorporated in the vest at thoracic and abdominal levels with two additional straps to avoid vertical movement. The vest is connected to two primary and two secondary oscillators; the signals pass then to two frequency/voltage converters, two filters, and an Analog to digital converter (PowerLab, ADInstruments). The system is connected by an RS232 cable to a computer equipped with an acquisition software (Chart, ADInstruments). A system of isolation is also included. The zero calibration is automatically performed. ECG and pneumo-tachogram (Fleisch no.1) are also monitored.The calibration of the RIP volume measured by the system is conducted based on the method developed by laboratory team PRETA-France in [21]. At the beginning of recording, the volunteer is asked to perform a number of successive high-amplitude respiratory cycles in order for Visuresp® to determine the range of amplitude to be measured. At the end of recording, RIP signal is considered as ƛ ∗ Abdominal activity + µ ∗ Thoracic activity. The activity means the discrete signal representing the cross section variation measured by the vest. The optimization of ƛ and µ is achieved by minimizing the error between the calculated RIP and the integral value of the airflow measured by pneumotachography.


## 2.3. Cardiorespiratory Model

In the present work, the cardiorespiratory model developed previously by our team in [22–24] is used in order to simulate the variations of respiratory, cardiorespiratory, and cardiac volumes. Briefly, the model includes two modules: respiratory (breathing pattern, alveolar volume, and pleural pressure) and cardiac (cardiac activity, chest wall mechanics, and volume variations (i.e., simulated IP)). The cardiac module is a wave generator. The model takes also into consideration modulation of left ventricle (LV) stroke volume due to breathing.


## 2.4. Parameter Identification and Sensitivity Analysis

Parameter identification and qualitative sensitivity analysis have been conducted using Berkeley Madonna® software in order to determine the most sensitive parameters. Those parameters are varied with the purpose of generating a set of physiological simulations to which our algorithm will be applied.


## 2.5. Simulations

According to sensitivity analysis, the most sensitive parameters are cardiac rate (HR) and breathing frequency (BF), which lead to frequency ratio (f). Also, amplitude ratio (a) between respiratory and cardiac activities is important. In our work, the variations of these parameters from a simulation to another generated a set of seven stationary and nonstationary amplitude modulated simulated RIP signals (more than 1000 cycles). It is worthy to note that these parameters are related to RIP components frequency and amplitude ratio that affects RIP empirical decomposition behavior. The simulations ratios are presented in Table 1.Table 1
Simulated situations. HR = heart rate, Fr: breathing frequency,f: frequency ratio, a: amplitude ratio.
Situation1/f1/aSituation 1 (stationary)4 (HR = 1; fr = 0.25)[0.1–0.17]Situation 2 (nonstationary)f 1 = 7.1; f 2 = 7.8 (HR1 = 1, HR2 = 2.2; Fr1 = 0.14, Fr2 = 0.28)[0.05–0.15]Situation 3 (non–stationary)f 1 = 3.6; f 2 = 7.2 (HR1 = 1, HR2 = 2; Fr1 = 0.28, Fr2 = 0.28)[0.04–0.16]Situation 4 (stationary)3 (HR = 0.75; Fr = 0.25)[0.1–0.15]Situation 5 (stationary)7 (HR = 1; Fr = 0.14)[0.05–0.15]Situation 6 (nonstationary)f 1 = 7.1; f 2 = 7.8 (HR1 = 1, HR2 = 1.8; Fr1 = 0.14, Fr2 = 0.23)[0.05–0.15]Situation 7 (nonstationary)(Progressive rise of HR: [1.4–3.3]; progressive rise of fr: [0.2–0.5])[0.05–0.08]Those simulations are created based on physiological scenarios where basic amplitude ratio (tidal volume/stroke volume) lies approximately in the range [5–20], and basic frequency ratio (breathing frequency/heart rate) lies approximately in the range [1/3–1/8]. For example, superficial high frequency breathing leads to rise in f and fall of a; rapid and profound breathing leads to rise in f and a; drug intoxication leads to a fall in f and a; and emotional stress leads to a rise in a and a fall in f. Note that the values presented in Table 1 are 1/f and 1/a. Moreover, the simulations include the transit from one physiological case to another, i.e., nonstationary situations.


## 2.6. Ensemble Empirical Mode Decomposition

Simulated RIP and real recorded inductive plethysmography signal Vth are decomposed by EEMD. The number of sifting iterations is limited to 10 in order to avoid over-sifting [20]. A set of 10000 white noise signals with amplitude of 1.6 times the r.m.s of RIP signal has been used. The r.m.s VR of RIP signal has been approximated to the r.m.s of the respiratory component:(2)VR≅1M∑k=1MSk2.S(k) is thek th sample of RIP signal. The r.m.s (VR) has been calculated over the number of all samples (M). The signal/noise ratio has been determined as follows:(3)σbVR=1.6,where σb is the added white noise standard deviation.The average of the corresponding IMFs from several tests is then calculated. In the final result, a signal/noise (εr ) ratio was acceptable if near 10%.(4)σrVc=0,1=εr,where σr is the residual noise standard deviation and Vc is the r.m.s of the cardiac component estimated using the amplitude ratio a. We can consider, as indicated in [20] that(5)σr=σbN.Hence,(6)N=σb∗aVR∗εr=1.6∗aεr.The physiological range ofa is between 5 and 20 (e.g., if a preliminary value of a = 6 is chosen, N is then equal to 104).


## 2.7. Double Decomposition and Cardiac Component Reconstruction

As illustrated in the work flowchart in Figure2, a new round of sifting has been applied to modes containing CR mixing.Figure 2
Work flowchart (double decomposition strategy). Vh is the cardiac volume.As mentioned, the algorithm generates a mixture of EEMD scales on certain IMF. This phenomenon is related to the criterion of stopping the sifting operation that takes place on a static number of iterations. This creates a loss of resolution in the separation of different components. This mixing has been detected automatically by Teager–Kaiser energy operator in order to detect the variations in energy present in a signal containing a mixing of two components. The calculated value at an indexk is the difference between the squared amplitude of the component at that k and the product of the amplitude at k − 1 with the amplitude at k + 1. When the Teager–Kaiser value exceeds the threshold 0.00005, a mixing is identified.It is noteworthy that several simulations have been conducted in the present work and all of them gave the same type of results with a scale mixing in IMF4 and IMF5. The sum of those two IMF gives a cardiac signal containing certain breathing harmonics. The solution proposed here is the double decomposition. Figure1(a) illustrated the interest thereof. While the original signal was located in zone 3, the signal of the sum of IMF4 and IMF5 is in Zone 1. The EMD can then separate effectively cardiac and respiratory components [25, 26]. The main parameters used in EMD application are type of extrema detection, type of envelope interpolation, maximum number of sifting iterations, maximum accepted value for the difference between upper and lower envelopes in an IMF, maximum number of obtained IMF, and masking-assisted decomposition. The values of parameters used in EMD were chosen as follows: extrema detection algorithm used in [19], Spline interpolation, 2000, zero, unlimited, and not used, respectively.


## 2.8. Stroke Volume Comparison

For simulated data, left ventricle volume signal (VIv) is used as the reference cardiac activity. Stroke volumes of extracted cardiac activity are hence compared with stroke volumes (SV) of simulated cardiac activity. The used statistical test is Bland and Altman. On the other hand, for real stationary data, cardiac component extracted by a thoracocardiography filter (FIR with order 200) has been used as the reference. More than 170 cycles were used for comparison. The zone in ± R-R interval around the inspiration-expiration transition has been excluded because it suffers from CR harmonics overlapping.The method of evaluation of concordance by Bland and Altman consists of the calculation of the difference between the values obtained by two different techniques as well as the average of those values. When the concordance between two series of values is perfect, the summation of arithmetic differences between the values is null. Also, the average of differences is defined as the bias. The variability of those differences indicates whether the two methods have different natures. In brief, Bland and Altman helps efficiently in judging the interchangeability between the two methods. In the present work, two thresholds of concordance are fixed based on the 95% interval of confidence (±30%). The differences are expressed as percentages of values averages. Bland and Altman graphs have been plotted via MedCalc®.


## 3. Results and Discussion

Figure3 illustrates the results obtained after the application of EEMD to a stationary simulated CR signal. IMF1 and IMF2 are likely composed of noise. The cardiac signal seems to spread over the IMF3, IMF4, and IMF5. Mode mixing effect does not appear. However, there is a respiratory component on IMF4 and IMF5 (scale mixing). The remaining IMFs concern respiration.Figure 3
IMF1-10 obtained after the application of EEMD to a stationary simulated RIP signal.The algorithm generates scale mixing on IMF4 and IMF5. The sum of these two IMF gives a cardiac signal containing some respiratory harmonics. The solution proposed here is the double decomposition. Figure4 illustrates the result of EMD applied to the sum of IMF4 and IMF5. The figure indicates an effective separation. Only IMF1 concerns the cardiac signal. In order to proceed to the reconstruction of the cardiac component, the considered IMF are IMF2 and IMF3 of the EEMD result and IMF1 of the EMD result.Figure 4
Double decomposition. EMD applied to the sum of IMF4 and IMF5 affected by scale mixing. From top to bottom: the sum of IMF4 and IMF5, the IMFs obtained by EMD and the residue (res).The test of Bland and Altman (Figure5) emphasizes the matching between the results and the criteria set out and searched in literature [2]. The relative difference between the simulated stroke volume (SV_model) and estimated stroke volume (SV_DD) has a value of less than ± 30% (95% confidence interval).Figure 5
Test of Bland and Altman used to compare stroke volumes estimated by double decomposition of stationary simulated RIP with simulated reference SV_model.The second objective of the present work is to analyze the performance of the algorithm in a case of variation in cardiorespiratory behaviour. Figure6 shows an absence of mode mixing in the application of EEMD to nonstationary RIP signal. It also shows that changes in heart rate and breathing are correctly detected on the IMFs.Figure 6
IMF1-10 obtained after the application of EEMD to a nonstationary simulated RIP signal.Figure7 illustrates the second decomposition applied to the sum of IMF4 and IMF5. Only IMF1 concerns the cardiac component. The test of Bland and Altman (Figure 8) certifies obtaining a limit of agreement consistent with the required criterion of interchangeability [2].Figure 7
Double decomposition. EMD applied to the sum of IMF4 and IMF5 (in Figure6) affected by scale mixing. From top to bottom: the sum of IMF4 and IMF5, the IMFs obtained by EMD and the residue (res).Figure 8
Test of Bland and Altman used to compare stroke volumes estimated by double decomposition of nonstationary simulated RIP with simulated reference SV_model.The third objective of the present work is to test the algorithm performance on a population of healthy subjects. Figure9 shows the result of EEMD application to a real cardiorespiratory signal. Figure 10 illustrates the double decomposition applied to the sum of IMF4 and IMF5. In order to reconstruct the cardiac wave, the following IMFs are considered: IMF 2 and 3 of the EEMD added to IMF1-3 of the EMD applied to the sum (IMF4 + IMF5).Figure 9
IMF1-10 obtained after the application of EEMD to a real cardiorespiratory signal.Figure 10
Double decomposition. EMD applied to the sum of IMF4 and IMF5 (in Figure 9) affected by scale mixing. From top to bottom: the sum of IMF4 and IMF5, the IMFs obtained by EMD and the residue (res).Figure11 illustrates the values SV_DD (estimated stroke volumes) cycle to cycle depending on their position relative to the beginning of the corresponding respiratory cycle (65 cardiac cycles for the subject n°7). We eliminate the cycles in the transition zone inspiration-expiration. The stroke volume estimated in the inspiratory phase and preceding the “forbidden zone” values are lower than those estimated at the end of expiration. This observation confirms that the findings are consistent with those described in [2, 3].Figure 11
Stroke volumes (cycle to cycle), subject n°7. (a) Several respiratory cycles of real cardiorespiratory signal. (b) SV calculated as a percentage of its average value. The red circles represent the SV estimated by double decomposition SV_DD. The blue stars represent the SV estimated by TCG [3]. The dotted lines surround the forbidden zone.
(a)(b)Figure12 illustrates the interchangeability between estimated stroke volumes (of real CR signals) and TCG (of real CR signals). This figure has been derived from the analysis of 170 cycles obtained on a set of seven healthy subjects and detected automatically. The automatic method has a good accuracy of cycle detection; 10 cycles were false positives and 2 cycles were false negatives.Figure 12
Bland and Altman comparison between stroke volume estimated by TCG and stroke volume estimated by double decomposition (real CR signals).For real stationnary data, we use thoracocardiography as our reference (170 cycles). In the case of nonstationary real data, we will have, first, to determine a suitable nonstationary protocol and second to compare our stroke volumes results to echocardiography. The algorithm has been efficient when applied to our simulated situations. We will have to raise this number in order to construct a population of simulated situations taking into account other sensitive parameters variability. This population will present the tests by which the algorithm application can be optimised.In addition, our algorithm should be optimised in order to account for stroke volume measurement in harmonics overlapping zone. Finally, we need to decrease the calculation time of the algorithm (2 hours for a recording of 2minutes). As a preliminary approach to find a solution to time consumption, the EEMD has been conducted in a “complementary” style; in other words, a complementary ensemble empirical mode decomposition (CEEMD) has been applied to signals instead of standard EEMD. This approach has been applied—as a trial—to nonstationary simulated scenario (situation 2) through achieving 50 pairs of complementary ensemble modes with positive and negative added white noises (noises that cancel each other). Hence, IMF is actually found out as the average of different counterpart IMFs extracted in different decompositions. The decomposition yielded good IMFs with eliminated residue noise similar to the level achieved by EEMD with a fast procedure that did not exceed (1/60) of the time consumed by EEMD.The novelty of the work lies mainly in two levels: (1) on the level of used instrument: the user-friendly, continuous, wearable, low-cost, and noninvasive monitoring of stroke volume via cardiorespiratory plethysmography. The main noninvasive methods currently used for stroke volume estimation are thoracic electrical bioimpedance [27] and transthoracic echocardiography [28]. However, in transthoracic echocardiography protocol, at least one intensivist in each unit who is an expert on transthoracic echocardiography should be available. Moreover, pitfalls and limits of transthoracic echocardiography need to be known by the users to avoid misinterpretations [29]. Also, in thoracic electrical bioimpedance protocol [30], body composition, skin status, electrodes adhesiveness, and obesity are critical factors that affect current penetration. (2) On the level of implemented algorithm (a) the suggested double empirical decomposition solves one of the challenging issues in wearable technology that is the absolute need of artefact detection and reduction [31]; the proposed method eliminates artefact adaptively based on the inherent physiological model. In addition, (b) the cardiac activity is extracted and separated from other activities in an intuitive accurate way, and (c) it is more adaptive than TCG that is triggered by heart rate and ECG.In the present work, several aspects have been taken into account for improvement of results (1) noise and artefacts are removed as undesired intrinsic mode functions. Hence, their effects are highly eliminated by discarding those IMFs as unwanted. (2) Also, the amplitude of the measured cardiorespiratory signals is calibrated via careful calculation of the correct coefficients of the weighted summation of the thoracic and the abdominal cross-sectional variations based on initial values of pneumotachograph airflow [32]. (3) In addition, the signal to noise ratio and the limit of iterations and noise additions in EEMD have been examined deeply with the intention of finding the best parameter values, since the method is very sensitive to noise level. (4) The comparison between our results and TCG results considers only the zones of signals that are studied by TCG in order to standardize the investigation.Several points represent sources of challenging complexity in the suggested method, starting from the phase of recording to the phases of calibration, simulation, decomposition, extraction, and comparison with reference: (1) the level of sensitivity of the implemented vest is a critical aspect as it determines the quality of recording and hence the correctness of measured cardiogenic oscillations. Without that required high level of correctness, the processed signals would not reflect a real cardiac activity; (2) the precise calibration of the weights of thoracic and abdominal variations is a highly important step. That step ensures a truthful representation of both cardiac and breathing activities in the processed signal, otherwise the whole subsequent work would be based on an erroneous cardiac reference; (3) the robust simulation of a cardiac activity embedded in a cardiorespiratory activity: the developed model represents a reference that must be realistic. For example, the simulated cardiac amplitude, cardiac frequency, and their variations through inspiration, transit, and expiration should absolutely be very well manipulated in analogy with physiology; (4) since the comparison—applied to real signals—considers the TCG reference, it was necessary to choose a test that judges interchangeability not superiority, this is due to the fact that a better method cannot be found out without a golden invasive method. The Bland and Altman test was therefore used; (5) the selected parameters in EEMD and in EMD in first and second decompositions, respectively, are too critical and should be cautiously tuned. The level of noise, the number of counterpart decompositions, the number of sifting, the noise tolerance, the IMF asymmetry tolerance, and the type of interpolation are all essential parameters that must be carefully adjusted in order to direct the technique accurately towards the desired cardiac components; (6) the automatic technique and the threshold that should be used to track the mode and scale mixing were the object of a successful exhaustive experimentation by the authors. EMD and EEMD are every so often “too empirical,” which yields to a difficulty of catching a general rule about the outcomes; (7) the skipping of real cardiorespiratory overlapping zones in inspiration-expiration transition periods by the suggested method is a disappointing choice that should be made since TCG cannot consist of a reference in those zones; and (8) finally, the processing cost and time are high. However, as indicated above, initial attempts of CEEMD indicate a big reduction of time and calculation cost.The future extended work will be the comparison with electrical bioimpedance and transthoracic echocardiography as they are the available tools for cardiorespiratory monitoring. Also, the investigation and optimization of the used device will be conducted in order to examine the used electrodes and their positions, as well as the possibility of better personalization of the vest [33]. Moreover, a pertinent automatic selection of cardiac IMF should be implemented.Table2 presents the main limitations of the available techniques that were targeted in the present work.Table 2
Comparison with other noninvasive approaches.
ReferenceUsed techniqueMain limitations[34]Finger plethysmographyHigh error. Bland and Altman (43%)[35]Transthoracic echocardiographyAt least one intensivist in each unit who is an expert on transthoracic echocardiography should be available. Moreover, pitfalls and limits of transthoracic echocardiography need to be known by the users to avoid misinterpretations[36]Electrical bioimpedanceBody composition, skin status, electrodes adhesiveness, and obesity are critical factors that affect current penetration.[37]ThoracocardiographyLack of adaptive nature as it depends on a band filter with cut-off values specified by heart rate and ECG.Our work:Cardiorespiratory inductive plethysmography + adaptive cardiac extraction. It does not need the presence of an expert. The device can be easily personalized and calibrated. The implemented algorithm is completely adaptive.


## 4. Conclusion

In this article, an implemented vest has been used according to the cardiorespiratory inductive plethysmography principle. Then, a double decomposition algorithm has been applied to simulated and real measured cardiorespiratory signal in order to extract cardiac activity. The simulated signals are generated by a previous model built by the authors. In the CR model, a first step has been achieved to the simulation of the physiological shapes of mechanical cardiac activity as well as respiratory waves. The double decomposition includes two steps: first step is an ensemble empirical mode decomposition. The second step is an empirical mode decomposition of selected components extracted from step 1. Both steps need the careful setting of specific parameters like the level of added noise and limit of iterations. The results show that this algorithm is a promising nonlinear method for efficient cardiogenic oscillations extraction in inductive plethysmography signal. It is more user friendly for continuous monitoring compared with available techniques such as transthoracic echocardiography. Also, it is more adaptive as compared with thoracocardiography.The suggested method consists an encouraging step towards noninvasive detection of cardiac activity in a continuous manner with low cost and need of expert presence. Also, it is a promising tool—for quantification of stroke volume—embedded in a fully mobile device.A further optimization step is necessary in order to measure the cardiac activity in CR harmonics overlapping zones. Moreover, the explored inductive vest should be studied for better personalization. Finally, a quantitative comparison should be made with approved techniques like electrical bioimpedance.



---

*Source: 9942367-2022-01-31.xml*