 Research
 Open access
 Published:
Instantaneous phase difference analysis between thoracic and abdominal movement signals based on complementary ensemble empirical mode decomposition
BioMedical Engineering OnLine volume 15, Article number: 112 (2016)
Abstract
Background
Thoracoabdominal asynchrony is often adopted to discriminate respiratory diseases in clinics. Conventionally, Lissajous figure analysis is the most frequently used estimation of the phase difference in thoracoabdominal asynchrony. However, the temporal resolution of the produced results is low and the estimation error increases when the signals are not sinusoidal. Other previous studies have reported timedomain procedures with the use of bandpass filters for phaseangle estimation. Nevertheless, the bandpass filters need calibration for phase delay elimination.
Methods
To improve the estimation, we propose a novel method (named as instantaneous phase difference) that is based on complementary ensemble empirical mode decomposition for estimating the instantaneous phase relation between measured thoracic wall movement and abdominal wall movement. To validate the proposed method, experiments on simulated time series and humansubject respiratory data with two breathing types (i.e., thoracic breathing and abdominal breathing) were conducted. Latest version of Lissajous figure analysis and automatic phase estimation procedure were compared.
Results
The simulation results show that the standard deviations of the proposed method were lower than those of two other conventional methods. The proposed method performed more accurately than the two conventional methods. For the humansubject respiratory data, the results of the proposed method are in line with those in the literature, and the correlation analysis result reveals that they were positively correlated with the results generated by the two conventional methods. Furthermore, the standard deviation of the proposed method was also the smallest.
Conclusions
To summarize, this study proposes a novel method for estimating instantaneous phase differences. According to the findings from both the simulation and humansubject data, our approach was demonstrated to be effective. The method offers the following advantages: (1) improves the temporal resolution, (2) does not introduce a phase delay, (3) works with nonsinusoidal signals, (4) provides quantitative phase estimation without estimating the embedded frequency of breathing signals, and (5) works without calibrated measurements. The results demonstrate a higher temporal resolution of the phase difference estimation for the evaluation of thoracoabdominal asynchrony.
Background
The clinical background of phase difference estimation in breathing
Breathing is one of the most essential processes for sustaining human life [1, 2]. In the quiet tidal breathing of a healthy human, the thoracic wall movement (TWM) and abdominal wall movement (AWM) occur nearly simultaneously and move synchronously (i.e., phase difference between 0° and 20°) [3–9]. The asynchrony (i.e., a phase difference) between the TWM and AWM has been found under different breathing conditions. For instance, various asynchronies have been found in people engaging in breathing exercise [7, 10]. In clinical practice, thoracoabdominal asynchrony has been measured to determine the underlying mechanism of patients’ respiratory system and to diagnose respiratory diseases [3, 4, 11]. Estimating thoracoabdominal asynchrony presents information about the time shift of AWM or other respiratory movements to the TWM. This shift demonstrates the disability of the respiratory mechanism. Specially, the estimation of asynchrony can be used to assess the function of the respiratory system. Thoracoabdominal asynchrony has been found in patients with chronic obstructive pulmonary disease [4], diaphragmatic weakness [12, 13], and upper airway obstruction [9], as well as after cardiothoracic surgery [3]. In addition, the “respiratory rate variability” (in analogy to heart rate variability or pulse rate variability) has not been established yet since it is difficult to determine the instantaneous respiratory rate unambiguously [14]. It is difficult because most of the measurements of thoracoabdominal asynchrony need to segment respiratory signal into breathbybreath [15]. As the recent advances in the field of the pulse rate variability research have demonstrated that the “instantaneous pulse rate variability” improved the temporal resolution and also helped in frequency exploration of cardiovascular autoregulation (e.g., autonomic nervous system function) [16], it was thought that the development of the instantaneous respiratory pattern could also be helpful for exploring how the breathing modulates autonomic nervous system.
The signal processing required for phase difference estimation
To estimate the phase relation of thoracoabdominal motions (i.e., TWM and AWM signals), various approaches that entail using respiratory inductance plethysmography (RIP) have been proposed [3, 17]. On the research path of thoracoabdominal motions, the Lissajous figure analysis (also called loop analysis) technique is a conventional method for estimating the phase difference between TWM and AWM signals. The technique, which has been widely used in clinical applications [3, 18, 19] and research [20, 21], first assumes that both the TWM and AWM signals are sinusoidal [20]. For calculating the phase difference, first, a filter is applied to these signals for noise reduction. Second, when the relation of these signals is drawn into an XY plot, in which the xaxis represents AWM and the yaxis represents TWM, a loop can be found. Finally, analyzing the loop then produces one value of a phase difference for each cycle of breathing. This temporal resolution (i.e., only a single value for each breath) of the phase analysis leaves room for improvement. Another problem is that the TWM and AWM signals are nonsinusoidal patterns in most cases, and estimation errors may increase when the signals deviate from being sinusoidal [22, 23]. To enhance the performance of phaseangle calculation in various respiratory patterns and to improve the temporal resolution of the estimation, a crosscorrelation method was used for phase analysis [23]. In this method, both the TWM and AWM signals are first partitioned into segments. The optimal timeshift of each segment is then searched for estimating the phase difference between the TWM and AWM signals. The breathbybreath segmentation and search of the optimal timeshift enable the intensive computation of phase difference calculation [15]. Other previous studies have reported timedomain procedures with the use of bandpass filters for analyzing asynchrony [15, 24]. The difficulties among all conventional approaches are that the thoracoabdominal motions measured using RIP have timevarying amplitudes and frequencies, and the measurements are easily impeded by noises such as movement artifacts (lowfrequency, highamplitude noise) and cardiogenic oscillations (highfrequency, lowamplitude noise). Thus, filtering noise from the original signal or reducing the sampling rate are necessary to increase the accuracy of phaseangle estimation [23]. Nevertheless, the bandpass filter, which is used to extract the original respiratory signals, induces a phase delay during phaseangle estimation, as mentioned in [15] (i.e., “the order of the linearphase bandpass FIR filter was set to 200, corresponding to a 2s time delay”). To prevent such a delay, a zerophase filter which processes the input data in both the forward and reverse directions can be used. The other studies [25, 26] demonstrated that an empirical mode decomposition (EMD)based filter can adaptively filter noises without any phase delay. Recently, Chen et al. [21] demonstrated a combination approach of complementary ensemble empirical mode decomposition (CEEMD) and Lissajous figure analysis for analyzing the phase difference of thoracoabdominal motions measured using RIP. The results of their 29humansubject experiment indicated a significant difference in the phase difference between thoracic breathing (TB) and abdominal breathing (AB). However, the temporal resolution of the phase analysis result was still limited to being per breathing cycle.
The aim and the organization of the current study
In this paper, we present a novel procedure for estimating the instantaneous phase difference (IPD) of RIP signals between TWM and AWM. This procedure offers five advantages: (1) improves the temporal resolution, (2) does not introduce a phase delay, (3) works with nonsinusoidal signals, (4) provides quantitative phase estimation without estimating the embedded frequency of breathing signals, and (5) works without calibrated measurements. Experiments on two simulated time series (i.e., sinusoidal and triangular signals [23]) and humansubject respiratory data were conducted for validation. In addition to the proposed IPD method, the newer version of loop analysis reported in [21] and the automatic phase estimation procedure (APEP) invented in [15] were applied to the data for comparison. Throughout this paper, the TWM signal measured using an RIP sensor is denoted as RIP_{TWM}, and the AWM signal measured using the RIP sensor is denoted as RIP_{AWM}. The mathematical symbol “±” is used to denote mean ± standard deviation.
Methods
A novel method for instantaneous phase difference estimation
To estimate the IPD from RIP_{TWM} and RIP_{AWM} signals, a novel method for IPD estimation was used as follows:

1.
CEEMD was first used to decompose original RIP signals into intrinsic mode functions (IMFs) [27].

2.
The main components of the RIP signals were then selected from the IMFs [28]; they are the main component of the RIP_{TWM} signal and the main component of the RIP_{AWM} signal.

3.
Finally, the IPD was computed according to the instantaneous phase of the main component extracted from the RIP_{TWM} and instantaneous phase of the main component extracted from the RIP_{AWM} signal.
Signal decomposition using CEEMD
CEEMD is a method derived from EMD [29], which is used to decompose biosignals into oscillatory modes (called IMFs) with different frequencies without introducing any phase delay [30]. Through EMD, a signal x(t) can be decomposed into IMFs and a final residue as follows:
where IMF _{ i }(t) represents the IMF, the i goes from 1 to A, r(t) is the residue of the raw data x(t), and t represents the time in seconds. The original signal is decomposed into sequential IMFs by using an EMDbased method. The EMDbased method is an iterative process and can be stopped whenever the goal of a user is met. For the proposed IPD method, one can have more or less IMFs extracted from an original as long as the main components are extracted in order to calculate the phase difference of the TWM and AWM signals. In this study, the K was chosen to be 10. Notably, the maximum number of the IMFs that can be extracted from a signal is limited by the sampling rate. The EMD process will be stopped when no more IMFs can be extracted from a signal.
The IMFs are the signals extracted from the original signal by using EMD; they satisfy two conditions [29] so that they can be used to calculate the instantaneous phase without introducing unwanted fluctuations induced by asymmetric wave forms. The final residue is the remainder of the original signal after the extraction of all the IMFs. The first condition requires that the number of extrema and the number of zero crossings in the data set either equal or differ at most by one. The second condition is that the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero at any point. The first condition is similar to the traditional narrow band requirements for a stationary Gaussian process. The second condition modifies the classical global requirement to a local one. In the EMD procedure, the frequency corresponding to the activities of the IMF obtained (i.e., 1st IMF, 2nd IMF, …, 10th IMF) becomes increasingly lower. This is because during the procedure, the IMFs are sequentially extracted from the original signal. The EMDbased method iteratively extracts the mean values of the envelope defined by the local maxima and local minima to produce an IMF. This process stops when the original becomes a monotonic function from which no more IMF can be extracted. The frequencies decrease with the increasing number of the IMF. For the final IMF, the frequency descends approximately to zero. The number of IMFs that can be extracted from a signal is limited by the sampling rate of the signal.
EMD consists of an inner loop that extracts an IMF from an input signal and an outer loop that subtracts the obtained IMF from the original signal and then submits the remaining signal into the inner loop for obtaining the next IMF (see Fig. 1). The procedure of the outer loop is as follows:

1.
The original signal x(t) is submitted into the inner loop to obtain an IMF.

2.
After an IMF is generated through the inner loop, the IMF is subtracted from the original signal (i.e., \(x_{i} \left( t \right) = x_{i  1} \left( t \right)  IMF_{i} \left( t \right)\)).

3.
The inner loop is restarted to obtain the next IMF according to the resulting signal x _{ i } (t).
EMD stops when no more IMFs can be extracted from a signal. The procedure of the inner loop is shown as follows:

1.
The local maxima and minima are detected in the input signal.

2.
The local maxima (all the peaks of the signal) and minima (all the valleys of the signal) are connected using cubic spline interpolation as the upper (i.e., u _{ i, g−1 } (t)) and lower (i.e., l _{ i, g–1 } (t)) envelopes, respectively.

3.
The mean envelope (i.e., m _{ i, g−1 } (t)) is calculated as the mean value of the upper and lower envelopes.

4.
The mean envelope is subtracted from the input signal (i.e., \(\it \it {\text{h}}_{i,\,g} \left( t \right) \, = \, h_{i,\,g  1} \left( t \right) \,  {\text{m}}_{i,\,g  1} (t)\)) and the procedure is repeated from the first step.
The inner loop iterates this procedure until the mean envelope approaches zero. Specifically, the inner loop is stopped on the basis of evaluation function \(\delta \left( t \right) = {{\left( {u_{i,\,g  1} \left( t \right) + l_{i,\,g  1} \left( t \right)} \right)} \mathord{\left/ {\vphantom {{\left( {u_{i,\,g  1} \left( t \right) + l_{i,\,g  1} \left( t \right)} \right)} {\left( {u_{i,\,g  1} \left( t \right)  l_{i,\,g  1} \left( t \right)} \right)}}} \right. \kern0pt} {\left( {u_{i,\,g  1} \left( t \right)  l_{i,\,g  1} \left( t \right)} \right)}}\). The inner loop is iterated until δ(t) < 0.2 for 95 % of the total length of δ(t), while δ(t) < 2 for the remaining 5 %. Additional details related to the stopping criteria of EMDbased methods are in [31]. The final subtraction result is denoted as an IMF.
CEEMD is an extension of EMD, which was proposed for solving the mode mixing problem and boundary effect problem of EMD [27]. The mode mixing problem is defined as a pattern of signals whose activities reside within the same frequency of different IMFs. When that happens, the IMFs extracted may not be monocomponent. This occurs when EMD fails to extract the activities within the same frequency to become a single IMF. CEEMD processes numerous instances of positive white noise and negative white noise to the signal during the extraction of IMFs. According to the literature [27, 32], adding white noise provides a uniform reference scale distribution and enhances EMD to prevent mode mixing. The boundary effect problem indicates that when generating the upper and lower envelopes during the EMD procedure, EMD uses zero as the amplitude of the envelope at the start and end of the signal. This manner of defining the values of the boundaries (i.e., the start and end of a signal) engenders an undesired pattern of the generated IMFs. CEEMD processes uniformly random values of boundaries. This was proved to prevent undesirable patterns in the boundaries of the extracted IMFs. The concept of CEEMD in decreasing the effect of mode mixing in EMD was applied by adding white noise [27] to the signal as follows:
where x(t) represents the raw data of the signal, w(t) denotes white noise, S _{+}(t) is the mixture of the raw data adding positive white noise, and S _{−}(t) is the mixture of the raw data adding negative white noise. The final IMF is the ensemble of IMF with both positive and negative noises. Figure 1 illustrates the CEEMD procedure. We added 50 instances of positive Gaussian white noise and 50 instances of negative Gaussian white noise. The white noise level was set to 25 % of the standard deviation of the respiratory signals.
Main component extraction
After the original signal was decomposed into IMFs, the main component of the RIP signal was determined according to the rule proposed in [28], in that the IMF with the highest energy density is the main component of the RIP signal. The energy density of an IMF is calculated as follows:
where \(\bar{E}_{i}\) is the energy density of the ith IMF and K is the total duration of IMF _{ i }(t). The frequency of an IMF is calculated as follows:
This equation was used to calculate the average period of the breathing cycles of the ith IMF by averaging all the periods determined according to the peaks of the ith IMF. P _{ i }[j] is the sample time of the jth peak of IMF _{ i }(t). The term \(\bar{T}_{i}\) represents the average period of the ith IMF, where a is the total peak number. The average frequency of an IMF was calculated as the inverse of its average period given by:
Instantaneous phase difference calculation
The instantaneous phases of a main component signal are calculated using normalized direct quadrature in the following equation [33]:
where NIMF _{i}(t) represents the normalized ith IMF. The IPD (θ _{ IPD }) between the instantaneous phase of AWM (θ _{ AWM }) and the instantaneous phase of TWM (θ _{ TWM }) is then calculated as follows:
The ith IMF can be normalized to NIMF _{i}(t) through the following procedure [33]:

1.
\(\left {IMF_{{{\text{i}},{\text{k}}}} \left( t \right)} \right\) is obtained for k = 0 by assigning \(\left {IMF_{\text{i}} \left( t \right)} \right\) to \(\left {IMF_{{{\text{i}},{\text{k}}}} \left( t \right)} \right\).

2.
All the local maxima of \(\left {IMF_{{{\text{i}},{\text{k}}}} \left( t \right)} \right\) are identified.

3.
All the maxima points are connected using cubic spline interpolation to generate the empirical envelope of \(\left {IMF_{{{\text{i}},{\text{k}}}} \left( t \right)} \right\) as e _{k}(t).

4.
IMF _{i, k}(t)/e _{k}(t) is assigned to \(\left {IMF_{{{\text{i}},{\text{k + 1}}}} \left( t \right)} \right\), and e _{k}(t) is used to pointwise normalize IMF _{i, k}(t).

5.
The procedure is repeated from step 1.
When all the \(\left {IMF_{{{\text{i}},{\text{k + 1}}}} \left( t \right)} \right\) values in step 5 are less than or equal to [−1, 1], the normalization is complete. \(\left {IMF_{{{\text{i}},{\text{k + 1}}}} \left( t \right)} \right\) is then assigned to NIMF _{i}(t) as the final result.
Experiment on simulated time series
The current study tested two types of signals (sinusoidal vs. triangular) × 6 degree levels (20°, 30°, 50°, 90°, 140°, and 170°) × 2 setups of noise (correlated noise vs. uncorrelated noise) as the simulated time series, resulting in 24 simulation setups.
Sinusoidal signals
For the first simulation, sinusoidal signals s _{1}(t) and s _{ 2 }(t) were respectively used to simulate pure TWM and AWM signals as follows:
where A is the amplitude (A = 1 in this study), T is the cycle period (T = 5 s in this study), and θ _{0} is the initial phase. In this simulation, θ _{0} = 20°, 30°, 50°, 90°, 140°, and 170° [3, 9, 13] were tested. The value T = 5 s was used because the normal range of the respiratory rate of adults is 12 (0.2 Hz) to 18 (0.3 Hz) breaths per minute [34–36].
Triangular signals
Triangular waves, s _{3}(t) and s _{4}(t), were also respectively used to simulate TWM and AWM (see Fig. 2). The formula for s _{3}(t) is shown as follows:
where m is the number of cycles (integers 0–59). The formula for s _{4}(t) is presented as follows:
where θ _{0} is the initial phase. In the current study, θ _{0} = 20°, 30°, 50°, 90°, 140°, and 170° [3, 9, 13] were tested.
Additive noise processes
Noises present in respiratory signals (e.g., electronic and sensor noises) measured using RIP (i.e., RIP_{TWM} and RIP_{AWM}) are nondeterministic; therefore, noise should be modeled as stochastic processes. Hence, in line with [15], a Gaussian white process with standard deviation 0.25 was used to model electronic and sensor noises, whereas the high amplitude with lowfrequency noise was used to model body movement artifacts. These generated noises were added to the signals to simulate the noise corruptions on these signals. Simulated TWM and AWM signals with correlated noises (i.e., n _{1}(t) = n _{2}(t)) and uncorrelated noises (i.e., n _{1}(t) ≠ n _{2}(t)) were tested.
Let w _{1}(t) and w _{2}(t) be the independent Gaussian white noises, and let σ _{ 1 }(t) and σ _{2}(t) be the highamplitude with lowfrequency noise that corrupts the TWM and AWM signals, respectively.
where A is the amplitude (A = 2), T is the cycle period (T = 75 s), and θ _{0} is the initial phase (θ _{0} = 180°) of the lowfrequency noise.
Experiment on human subjects
Ethics statement
The data used in this section was from the research project “Pacedrespiratory induced heart rate variability and cardiac output evaluation (Protocol No: 100015E),” approved by the Institution Review Board of the National Taiwan University Hospital Hsinchu Branch. The committee was organized under and operated in accordance with the Good Clinical Practice guidelines and governmental laws and regulations. Written Informed consents were obtained from all subjects before the experiment. The experiment and the use of the data obtained from human subject were performed in accordance with the approved protocol.
Subjects
There were 50 subjects that were college students selected from a university in Taiwan, ranging in age between 20 and 23 (21 ± 1 years; 42 men, 8 women), body height between 157.5 and 188 cm (172.05 ± 6.68 cm), body weight between 43 and 87 kg (64.61 ± 9.77 kg), body mass index (BMI) between 15.21 and 29.05 kg/m^{2} (21.77 ± 2.69 kg/m^{2}), thoracic circumference between 76 and 105 cm (88.08 ± 6.27 cm), and abdominal circumference between 63.50 and 89 cm (77.15 ± 7.52 cm), were asked to perform TB and AB during the experiment. The subjects were all instructed not to take alcoholic, caffeinecontaining drinks or big meals 4 h before the experiment.
Experimental procedure
The experiment was carried out in a small, quiet room (7.60 × 3.20 m) without people. The subject wore two RIP sensors (RIPmate Adult Thorax Alice 5 Inductance Kit, Ambu Inc., Denmark) and remained in a seated position in an office chair (0.50 × 0.51 m, height 0.43 m) during the experiment, and was instructed to perform TB and AB. The RIP_{TWM} sensor was worn below the axilla to record TWM, and the RIP_{AWM} sensor was placed on the navel to record AWM. Both RIP signals were acquired by a data acquisition hub (NI SCB68, National Instruments, USA) and a data acquisition card (NI USB 6255, National Instruments, USA) at a sampling rate of 50 Hz and were subsequently transferred to a computer (Acer Veriton M2610; processor: Intel Core i32120 3.3G/3M/65W; memory: 4 GB DDR31066; operating system: Microsoft Windows 7 Professional 64 bit). The procedure of the experiment was as follows:

1.
The experiment started with an instruction (“Please follow the metronome in the monitor to breath with thoracic breathing/abdominal breathing”) presented on the computer screen (ViewSonic VE700, 17in, 1280 × 1024in resolution) for 5 s. The subject then performed TB at a breathing rate (controlled by a metronome on the screen) of 12 breaths per minute (i.e., 0.2 Hz [34–36]) for 5 min. Because TB is the normal breathing type for humans, the subject was asked to breathe normally and focus on the movement of the chest.

2.
The subject was instructed to learn AB by using a respiratory signal monitoring system (same as that used in [21]). The system taught the subject AB through three steps.

a.
The system instructed the subject to place the left hand on the chest and place the right hand on the abdomen. In addition, the system instructed the subject to focus on the elevation and degradation of the AWM in the second step.

b.
The system instructed the subject to inhale deeply through the nose and hold the breath for 2 s, and then exhale through the mouth until the end of the expiration.

c.
The system then instructed the subject to repeat the second step continuously for 10 min.

a.

3.
After the training session, the subject performed AB at a breathing rate of 12 breaths per minute for an additional 5 min.
The RIP data were recorded during both the TB and AB periods. The program controlling the instructions and data acquisition was developed using LabVIEW platform (LabVIEW 2012, National Instruments, USA).
Data analysis and comparison
The result obtained using the proposed IPD method was compared with that derived using the improved version of loop analysis [21] and the APEP described in [15]. The loop analysis procedure is described as follows:

1.
The method first extracts the TWM and AWM signals from the RIP_{TWM} and RIP_{AWM} as their main components (i.e., the IMF with the highest energy among all the IMFs decomposed from the original signal [28]), respectively.

2.
The breathing cycles of the TWM and AWM signals are identified automatically by detecting the valleys of the signal by using a computer program.

3.
Loop analysis [3] is applied to the extracted TWM and AWM signals. The relation of two timedependent functions x(t) and y(t) is displayed in an x–y plot. For each breathing cycle, a Lissajous figure of the rib cage versus abdomen motion signals is drawn. The phase difference (θ) is then calculated using the following equation: sin(θ) = m/s, where m is the length of abdomen motion excursion, which is parallel to the abscissa in the Lissajous figure in the midrib cage motion, and s is the length of the overall abdomen motion excursion.
The APEP procedure is described as follows [15]:

1.
The TWM and AWM signals are submitted into a zerophase finiteimpulse response bandpass filter with parameters set to a highpass frequency of 0.4 Hz, a highstop frequency of 0.45 Hz, a lowpass frequency of 0.15 Hz, and a lowstop frequency of 0.1 Hz.

2.
The resulting signals are converted using a binary converter that operates pointwise on its input signal, \(\tilde{s}\left[ n \right]\), as follows:
$$\tilde{s}\left[ n \right] = \left\{ {\begin{array}{*{20}c} {0, \quad if \tilde{s}\left[ n \right] < 0} \\ {1, \quad if \tilde{s}\left[ n \right] > 0} \\ \end{array} } \right.$$(16) 
3.
The two resulting signals are submitted into an exclusiveOR gate that operates pointwise on its input signals.

4.
The phase difference (θ) of the resulting signal is estimated using a lowpass filter that operates pointwise on its input signal, u[n], as follows:
$$\theta [n] = \frac{1}{L}\mathop \sum \limits_{k = 0}^{L  1} u[n  k]$$(17)
where L denotes the number of sample points in the 5slong sliding window (L = 251).
Bland–Altman analysis was applied to the humansubject data to examine the relation between the phase difference estimation methods [37, 38]. The twostep procedure was carried out according to the recommendations of [39]. First, the correlation between the estimated values (i.e., the equation of the linear relationship and correlation coefficient) was investigated. Second, the relative differences between each pair of measures were plotted against the mean of the pair to determine whether most of the pairs fall within the 95 % limit of the agreement bound.
All numerical analyses in this study were performed using LabVIEW 2012 (National Instruments, USA). Both simulated and humansubject data were sampled at 50 Hz. The significance level of the entire statistical hypothesis tests used in this study, if not stated otherwise herein, was set to 0.05. All statistical analyses were performed using commercial statistics software (Statistical Package for Social Science, Version 22, SPSS Inc., USA).
Results
Experiment on simulated time series
Figure 2 shows the 60s sample results of the phase difference estimations on two triangular signals s_{3}(t) and s_{4}(t) with correlated (left part) noise (denoted by n_{1}(t)) and uncorrelated (right part) noise (denoted by n_{1}(t) and n_{2}(t)) obtained using the IPD method, loop analysis, and APEP. The frequency of s_{3}(t) and s_{4}(t) was set to 0.2 Hz, and the phase difference θ _{0} was set to 20°. The dashed line represents the gold standard (i.e., θ _{0} = 20^{∘}). Figure 2 indicates that the average value of the IPD result is closer to θ _{0} compared with the average value of loop analysis and APEP results (i.e., 19°, 15° and 19° for IPD, loop analysis, and APEP with correlated noise, respectively; and 20°, 17° and 30° for IPD, loop analysis, and APEP with uncorrelated noise, respectively). This figure also demonstrates that loop analysis generates only one output value per cycle of breathing at the end of each cycle (i.e., 11 output values in the 60s sample result), whereas both the IPD method and APEP provide pointwise output values (i.e., 60 s × 50 samples/s = 300 samples). The ratio of the number of output values between loop analysis and each of the other two methods is 11:300.
This difference indicates a higher temporal resolution of the phase estimation achieved using the IPD method and APEP, compared with the temporal resolution reached using loop analysis. Table 1 indicates that the IPD method was accurate in all the simulation setups. The results of each of the simulation setups were subjected to onesample t tests to determine whether significant differences existed between the gold standard \(\theta_{0}\) and the phase difference estimation results obtained using all the methods. For all the degree levels tested and types of signals used, the absolute differences between the gold standard \(\theta_{0}\) and mean value of the IPD results were 0°. The standard deviations of the IPD results were between 0.13° and 0.49°. No significant difference between θ _{0} and the IPD estimated phase difference was found in any of the 24 simulation setups. However, the absolute differences between the gold standard \(\theta_{0}\) and mean value of the loop analysis results were between 0° and 8°. The standard deviations of the loop analysis results were between 0.1° and 1.04°. Significant differences between θ _{0} and the phase difference estimated in loop analysis were found in 21 of the 24 simulation setups. Finally, the absolute differences between the gold standard \(\theta_{0}\) and mean value of the APEP results were between 0° and 20°. The standard deviations of the APEP results were between 0.26° and 0.87°. Significant differences between θ _{0} and the phase difference estimated in APEP were determined in 10 of the 24 simulation setups. For half of the ten simulation setups, the θ _{0} value was larger than 140°. The results indicated an inaccuracy of the APEP in phase degree θ _{0} larger than 140°, regardless of the signal or noise types. The absolute differences between the gold standard \(\theta_{0}\) and the mean values of the loop analysis results for the triangular signals with correlated noise and uncorrelated noise were higher than the absolute differences for the sinusoidal signals with correlated noise in 10 of the 12 simulation setup. However, no similar tendency was observed in the IPD results or APED results.
Experiment on humansubject
In total, two types of breathing (i.e., TB and AB) × 2 (RIP_{TWM} and RIP_{AWM} signal) × 50 (subjects) = 200 rows of the raw data were collected during the experiment. The RIP_{TWM} and RIP_{AWM} signals were decomposed into IMFs by using the CEEMD method described in the “Methods” section. According to the procedure described in the “Methods” section, the sixth IMFs of both the RIP_{TWM} and RIP_{AWM} signals were determined to be the main component (i.e., TWM) of the RIP_{TWM} signal and the main component (i.e., AWM) of the RIP_{AWM} signal, respectively. The frequencies of the extracted TWM signals in TB and the AWM signals in AB were nearly 0.2 Hz. Sample results of the phase difference analysis obtained using our method for IPD estimation, loop analysis, and APEP are presented in Fig. 3. These results indicate that for the IPD method and APEP, the temporal resolution of the result is higher. In other words, the IPD and APEP provided the phase difference estimation at a sampling rate of 50 Hz (i.e., there were 60 s × 50 samples/s = 300 samples for both the IPD and APEP, as shown in Fig. 3), whereas the result of loop analysis was only cycle based (i.e., less than 0.5 Hz; there were 10 output values for loop analysis, as shown in Fig. 3). In addition, it appears that all three methods could still estimate the phase difference between the RIP_{TWM} signal and the RIP_{AWM} signal even when the RIP_{TWM} signal is markedly lower in amplitude (i.e., 0.001 mV; first panel of Fig. 3) than the RIP_{AWM} signal (i.e., 0.004 mV; second panel of Fig. 3).
Figure 4 presents the results of all three phase analysis methods. The pairedsample t test results indicate that the estimated phase difference between RIP_{TWM} and RIP_{AWM} signals in AB was significantly greater than that in TB, for the estimation conducted through loop analysis (35 ± 18.77° vs. 22 ± 16.16°, p < .001), the proposed IPD method (24 ± 6.29° vs. 19 ± 5.68°, p < .001), and the APEP (55 ± 38.40° vs. 38 ± 30.08°, p < .05). The correlation analysis results between the IPD method and the other two methods are provided in the upper panels of Figs. 5 and 6 whereas the Bland–Altman analysis results between the IPD method and the other two methods are provided in the lower panels of Figs. 5 and 6.
The data used for Bland–Altman analysis were logarithmically transformed from the original data to ensure a normal distribution, according to the suggestions provided in [38]. Both the result of loop analysis (R^{2} = 0.80 under TB and R^{2} = 0.80 under AB) and the result of APEP (R^{2} = 0.39 under TB and R^{2} = 0.68 under AB) were positively correlated with the result of the proposed IPD method. The Bland–Altman analysis, which examines whether data points that represent the difference between IPD and another method fall within the means ± 2 × standard deviation, revealed that 94 % of the data points in TB and 98 % of the data points in AB for IPD versus loop analysis fell within the means ± 2 × standard deviation. By contrast, 94 % of the data points in TB and all the points in AB for IPD versus APEP fell within the means ± 2 × standard deviation.
Discussion
Experiment on simulated time series
Sinusoidal signals are usually used to guide subjects’ inhaled/exhaled volume in the literature [40, 41] because the sinusoidal pattern appears to be similar to respiratory movements of normal subjects [40]. However, respiratory movements are sometimes more triangular in shape than sinusoidal, for instance, in lightly anesthetized rhesus monkeys [23]. In reality, true breathing is neither always triangular nor always sinusoidal. The pattern may be affected by the breathing condition or the disease of a subject. Hence, both sinusoidal wave and triangular wave were used in our simulation. Moreover, according to our simulation result, it is suggested to use the IPD method or the APEP than loop analysis when the target respiratory signal is more triangular than sinusoidal.
The loop analysis was a conventional approach that is proved to be useful and thus widely used in the literature for estimating phase difference. However, the results of the simulated time series show that loop analysis, because of its limitation, is low in temporal resolution. Results of the proposed IPD method and APEP indicate a higher temporal resolution than that of loop analysis. The inaccuracy of loop analysis for the triangular signals is in line with the findings of previous studies [15, 23]. The IPD method and APEP were determined to be more accurate than loop analysis regarding the absolute differences between the gold standard θ _{0} and the mean of the estimated values in most of the simulation setups. However, the finding of the inaccuracy of the APEP in estimating θ _{0} larger than 140° suggests that the IPD may be a more favorable choice regarding methods for detecting a high degree of phase difference. The estimation of high degrees of phase difference (e.g., from 150° to 179°) is vital because they were found to be highly correlated with certain types of respiratory diseases [3, 9, 13]. Our finding suggests that the competitive performance of the proposed IPD method in estimating θ _{0} of sinusoidal and triangular signals was corrupted with correlated/uncorrelated noises, because for all simulation setups, no significant difference was found between the gold standard θ _{0} and IPD result. Furthermore, the absolute differences between the gold standard \(\theta_{0}\) and mean value of the IPD results were small, and the standard deviations of the IPD results were small. These results indicate the accuracy and the consistency of the proposed method.
The magnitude of the absolute difference found between the gold standard θ _{0} and the estimation results observed in the simulation is essential and required to be taken into consideration when developing clinical applications. It is because the magnitude of this difference could be proportional to the prediction error in medical diagnosis. This magnitude is not negligible even when it is only 1°, considering this 1° error may affect the analysis result largely in some applications. For instance, the difference in phase angles found in normal subjects and moderate chronic obstructive pulmonary disease patients was just 3.5° [4]. In this case, it is hard to determine the patients with moderate chronic obstructive pulmonary disease using phase angle. Figure 7 demonstrates the effects of CEEMD and a zerophase bandpass filter on two sinusoidal signals with correlated noises and phase difference θ _{0} = 50°. Both the result of the CEEMD filter and the zerophase bandpass filter appear to be similar to the original signals and have no phase delay to the original signals. Moreover, we also found that CEEMD filter appears to be lesser noisy than the zerophase bandpass filter in Fig. 7. In addition, despite the fact that the zerophase filter can be implemented for the APEP method to prevent phase delay, the filter must have time samples available before to start processing. This could be a limitation when developing realtime applications [42].
Experiment on humansubject
Because the subjects were requested to perform paced breathing at a breathing rate of 12 breaths per minute (i.e., 0.2 Hz), the 0.2Hz frequency of the CEEMDextracted TWM signal in TB and the AWM signals in AB seems to be correct. In line with the finding reported by previous research [5–7, 10], all three methods estimated a greater phase difference in AB than in TB (Fig. 4). This finding suggests that the proposed IPD in estimating phase difference in real respiratory data shows a pattern similar to that from using the other two methods. The Bland–Altman and correlation analysis results also indicate that the estimation result of the proposed IPD method shares the same trend with those of the other two methods regarding real respiratory data.
However, the phase difference results estimated through loop analysis (i.e., 22 ± 16.16° for TB and 35 ± 18.77° for AB) and the IPD method (i.e., 19 ± 5.68° for TB and 24 ± 6.29° for AB) are in line with previous findings, in that the phase difference of TWM and AWM under the TB and AB conditions fell within 0° to 40° for healthy subjects [4–9]. This renders the result of the APEP (i.e., 38 ± 30.08° for TB and 55 ± 38.40° for AB) questionable. Furthermore, the IPD had the smallest standard deviation (i.e., 5.68 for TB and 6.29 for AB) among the three methods, whereas the APEP had the largest standard deviation (i.e., 30.08 for TB and 38.40 for AB). This implies that the IPD is a more stable approach to estimating IPD compared with the other two methods.
The correlation analysis results shown in the upper panels of Figs. 5 and 6 present positive correlations between IPD and the other two methods; the R^{2} between the IPD method and APEP (R^{2} = 0.39) was lower than the R^{2} between the IPD method and loop analysis (R^{2} = 0.80) in TB. The large standard deviation (30.08° for TB and 38.40° for AB) and the unexpectedly large phase difference (38° for TB and 55° for AB) of the APEP result (compare with those in the literature [4–9]) may be the reason that the IPD result is correlated with loop analysis result but less correlated with the APEP result. Furthermore, although the IPD result is less correlated with the APEP result compared with that of loop analysis, the IPD results from Bland–Altman analysis were still favorable (lower panels of Figs. 5, 6). Most of the data points of the IPD result fell within the means ± 2 × standard deviation of the loop analysis result and the APEP result. This means that for the real data, the pattern of the phase estimated using the IPD method is similar to those of the other two methods. The negative correlations are depicted in the lower panels of Figs. 5 and 6. In this data set, the exact value of the phase difference estimated using IPD tends to be higher than those estimated using the other two methods when all three methods yield a relatively low phase difference. By contrast, when all three methods yield a relatively high phase difference, the IPD tends to yield lower values than those estimated through loop analysis and the APEP. This leaves room for further research because there is currently no ground truth of the “phase difference”, a quantity that cannot be measured directly by using any type of device. The only current certainty is that the result generated by the IPD method shares a similar tendency with those of the other two methods for real data.
Research limitations
Despite the experimental results appeared to be promising, the proposed IPD method has certain limitations. First, the use of EMDbased data processing in our method introduces a certain amount of complexity of using the method (same as the recent version of Lissajous figure analysis method proposed in [21]). The EMDbased data processing may sometimes require a researcher to manually determine multiple parameters (e.g., parameters related to the stopping criteria of the iterative decomposition processes) to avoid mode mixing problem and the boundary effect problem [27]. However, in ordinary these parameters can be intuitively determined based on the decomposition result. Another limitation is that the iterative process of the EMDbased data processing may be time consuming [43]. Notably, the order of the computational complexity of the EMD has been shown to be equivalent to FFT, and can be optimized in order to operate in realtime applications [43, 44]. For applications and data analysis tasks operate offline, parallelization [45] and the use of graphics processing unit [46] can further improve the speed of data processing.
Conclusions
Because thoracoabdominal asynchrony is often used to discriminate respiratory diseases in clinics, the current study presents a novel method (called IPD estimation) based on CEEMD for estimating the phase relation between TWM and AWM signals measured using RIP sensors. Both simulation and humansubject experimental respiratory data show that the proposed IPD estimation method demonstrated competitive performance compared with the recent version of Lissajous figure analysis and the APEP. We believe that the findings of the current study can provide future studies with a new phase difference estimation method that (1) has high temporal resolution, (2) does not introduce phase delay, (3) works with nonsinusoidal signals, (4) provides quantitative phase estimates without estimating the embedded frequency of the breathing signals, and (5) works without calibrated measurements. Thus, the accuracy of the research result may be enhanced. Furthermore, according to the higher temporal resolution of the estimation of phase difference, more information related to the interaction between TWM and AWM can be revealed. We believe that IPD estimation with higher temporal resolution could be a useful aid to the exploration of how the breathing pattern modulates cardiovascular system and autonomic nervous system in the near future.
Abbreviations
 AB:

abdominal breathing
 AWM:

abdominal wall movement
 APEP:

automatic phase estimation procedure
 CEEMD:

complementary ensemble empirical mode decomposition
 EMD:

empirical mode decomposition
 IMFs:

intrinsic mode functions
 IPD:

instantaneous phases difference
 RIP:

respiratory inductance plethysmography
 TB:

thoracic breathing
 TWM:

thoracic wall movement
References
John BW. Respiratory physiology: the essentials. 9th ed. Alphen aan den Rijn: Wolters Kluwer; 2012.
Ganong WF. Review of medical physiology. 22nd ed. New York: Mcgrawhill; 2005.
Hammer J, Newth CJL. Assessment of thoracoabdominal asynchrony. Paediatr Respir Rev. 2009;10:75–80.
Chien JY, Ruan SY, Huang YCT, Yu CJ, Yang PC. Asynchronous thoracoabdominal motion contributes to decreased 6min walk test in patients with COPD. Respir Care. 2013;58:320–6.
Gallego J, Benammou S, Vardon G, Chambille B, Denjean A, Lorino H. Influence of thoracoabdominal pattern of breathing on respiratory resistance. Respir Physiol. 1997;108:143–52.
CancellieroGaiad KM, Ike D, Pantoni CB, BorghiSilva A, Costa D. Respiratory pattern of diaphragmatic breathing and pilates breathing in COPD subjects. Braz J Phys Ther. 2014;18:291–9.
Tomich G, França D, Diório A, Britto R, Sampaio R, Parreira V. Breathing pattern, thoracoabdominal motion and muscular activity during three breathing exercises. Braz J Med Biol Res. 2007;40:1409–17.
Allen JL, Greenspan JS, Deoras KS, Keklikian E, Wolfson MR, Shaffer TH. Interaction between chest wall motion and lung mechanics in normal infants and infants with bronchopulmonary dysplasia. Pediatr Pulmonol. 1991;11:37–43.
Sivan Y, Deakers TW, Newth CJ. Thoracoabdominal asynchrony in acute upper airway obstruction in small children. Am Rev Respir Dis. 1990;142:540–4.
Vieira D, Mendes L, Elmiro N, Velloso M, Britto R, Parreira V. Breathing pattern and thoracoabdominal motion during breathing exercises in healthy subjects. Eur Respir J. 2014;42:1329.
Reber A, Bobbia S, Hammer J, Frei F. Effect of airway opening manoeuvres on thoracoabdominal asynchrony in anaesthetized children. Eur Respir J. 2001;17:1239–43.
Bedenice D, Mazan MR, Kuehn H, Hoffman AM. Diaphragmatic paralysis due to phrenic nerve degeneration in a llama. J Vet Intern Med. 2002;16:603–6.
Goldman MD. Interpretation of thoracoabdominal movements during breathing. Clin Sci. 1982;62:7–11.
Cysarz D, Zerm R, Bettermann H, Frühwirth M, Moser M, Kröz M. Comparison of respiratory rates derived from heart rate variability, ECG amplitude, and nasal/oral airflow. Ann Biomed Eng. 2008;36:2085–94.
Motto AL, Galiana HL, Brown KA, Kearney RE. Automated estimation of the phase between thoracic and abdominal movement signals. IEEE Trans Biomed Eng. 2005;52:614–21.
Chang CC, Hsu HY, Hsiao TC. The interpretation of very high frequency band of instantaneous pulse rate variability during paced respiration. Biomed Eng Online. 2014;13:1–11.
Musante G, Schulze A, Gerhardt T, Everett R, Claure N, Schaller P, et al. Proportional assist ventilation decreases thoracoabdominal asynchrony and chest wall distortion in preterm infants. Pediatr Res. 2001;49:175–80.
Newth C, Hammer J. Measurements of thoracoabdominal asynchrony and work of breathing in children. Prog Respir Res. 2005;33:148–56.
Allen JL, Wolfson MR, McDowell K, Shaffer TH. Thoracoabdominal synchrony in infants with airflow obstruction. Am Rev Respir Dis. 1990;141:337–42.
Selbie R, Fletcher M, Arestis N, White R, Duncan A, Helms P, Duffty P. Respiratory function parameters in infants using inductive plethysmography. Med Eng Phys. 1997;19:501–11.
Chen YC, Hsiao TC, Chen JL. Better thoracoabdominal synchrony in abdominal breathing: evidence from complementary ensemble empirical mode decompositionbased Lissajous figure analysis. J Med Imaging Health Inf. 2015;5:400–5.
Laouani A, Rouatbi S, Saguem S, Calabrese P. Thorax and abdomen motion analysis in patients with obstructive diseases. J Pulm Respir Med. 2016;6:313.
Prisk GK, Hammer J, Newth CJL. Techniques for measurement of thoracoabdominal asynchrony. Pediatr Pulmonol. 2002;34:462–72.
Aoude AA, Kearney RE, Brown KA, Galiana HL, RoblesRubio CA. Automated offline respiratory event detection for the study of postoperative apnea in infants. IEEE Trans Biomed Eng. 2011;58:1724–33.
Chang CC, Hsiao TC, Hsu HY. Frequency range extension of spectral analysis of pulse rate variability based on Hilbert–Huang transform. Med Biol Eng Comput. 2014;52:343–51.
Chang KM. Arrhythmia ECG noise reduction by ensemble empirical mode decomposition. Sensors. 2010;10:6063–80.
Yeh JR, Shieh JS, Huang NE. Complementary ensemble empirical mode decomposition: a novel noise enhanced data analysis method. Adv Adapt Data Anal. 2010;2:135–56.
Chen YC, Hsiao TC, Hsu JH, Chen JL. Breathing pattern recognition of abdominal wall movement by using ensemble empirical mode decomposition. Adv Adapt Data Anal. 2014;6:1450002–18.
Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, Yen NC, Tung CC, Liu HH. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proc R Soc Lond A. 1998;454:903–95.
Chang CC, Kao SC, Hsiao TC, Hsu HY. Assessment of autonomic nervous system by using empirical mode decompositionbased reflection wave analysis during nonstationary conditions. Physiol Meas. 2014;35:1873–83.
Rilling G, Flandrin P, Goncalves P. On empirical mode decomposition and its algorithms. In: IEEEEURASIP workshop on nonlinear signal and image processing. 2003. p. 3:8–11.
Wu Z, Huang NE. Ensemble empirical mode decomposition: a noiseassisted data analysis method. Adv Adapt Data Anal. 2009;1:1–41.
Huang NE, Wu Z, Long SR, Arnold KC, Chen X, Blank K. On instantaneous frequency. Adv Adapt Data Anal. 2009;1:177–229.
Ogawa Y, Iwasaki K, Shibata S, Kato J, Ogawa S, Oi Y. Different effects on circulatory control during volatile induction and maintenance of anesthesia and total intravenous anesthesia: autonomic nervous activity and arterial cardiac baroreflex function evaluated by blood pressure and heart rate variability analysis. J Clin Anesth. 2006;18:87–95.
Lourens MS, van den Berg B, Aerts JG, Verbraak AF, Hoogsteden HC, Bogaard JM. Expiratory time constants in mechanically ventilated patients with and without COPD. Intensive Care Med. 2000;26:1612–8.
Reyes B, Reljin N, Kong Y, Nam Y, Chon K. Tidal volume and instantaneous respiration rate estimation using a smartphone camera. IEEE J Biomed Health Inform. 2016:1–15.
Bland JM, Altman D. Statistical methods for assessing agreement between two methods of clinical measurement. Lancet. 1986;327:307–10.
Giavarina D. Understanding bland altman analysis. Biochem Med. 2015;25:141–51.
Naidu MUR, Reddy BM, Yashmaina S, Patnaik AN, Rani PU. Validity and reproducibility of arterial pulse wave velocity measurement using new device with oscillometric technique: a pilot study. Biomed Eng Online. 2005;4:1–14.
da Silva Junior EP, Esteves GP, Dames KK, de Melo PL. A telemedicine instrument for Internetbased home monitoring of thoracoabdominal motion in patients with respiratory diseases. Rev Sci Instrum. 2011;82:014301.
Simons E, Wood RA. Effect of body size on breathing pattern and fineparticle deposition in children. Pediatrics. 2005;116:555.
Lyons RG. Understanding digital signal processing. 2nd ed. Upper Saddle River: Pearson Education; 2010.
Wang YH, Yeh CH, Young HWV, Hu K, Lo MT. On the computational complexity of the empirical mode decomposition algorithm. Phys A. 2014;400:159–67.
Lu TC, Chen PY, Yeh SW, Van LD. Multiple stopping criteria and highprecision EMD architecture implementation for Hilbert–Huang transform. In: IEEE biomedical circuits and systems conference proceedings. 2014. p. 200–203.
Chang LW, Lo MT, Anssari N, Hsu KH, Huang NE, Hwu WmW. Parallel implementation of multidimensional ensemble empirical mode decomposition. In: IEEE international conference on acoustics, speech and signal processing. 2011. p. 1621–24.
Chen D, Li D, Xiong M, Bao H, Li X. GPGPUaided ensemble empiricalmode decomposition for EEG analysis during anesthesia. IEEE Trans Inf Technol Biomed. 2010;14:1417–27.
Authors’ contributions
YCC carried out data acquisition, data analysis and prepared the manuscript. TCH obtained ethics approval and arranged resources, materials, and analysis tools. Both authors read and approved the final manuscript.
Acknowledgements
N. E. Huang (academician of Academia Sinica and member of National Academy of Engineering) is gratefully acknowledged for empirical mode decomposition studies.
Competing interests
The authors declare that they have no competing interests.
Availability of data and material
The informed consent received from subjects does not include an explicit permission of uploading the collected data to public database for anonymous download. Hence, we believe that to make the data available to all interested researchers upon request should be more appropriate. The reviewers and readers may request the data by contacting: TzuChien Hsiao (email: labview@cs.nctu.edu.tw).
Ethics approval and consent to participate
The data used in this paper was from the research project “Pacedrespiratory induced heart rate variability and cardiac output evaluation (Protocol No: 100015E),” approved by the Institution Review Board of the National Taiwan University Hospital Hsinchu Branch. The committee was organized under and operates in accordance with the Good Clinical Practice guidelines and governmental laws and regulations. Written Informed consents were obtained from all subjects before the experiment. The experiment and the use of the data obtained from human subject were performed in accordance with the approved protocol.
Funding
This work was fully supported by the Taiwan Ministry of Science and Technology under Grant numbers MOST 1032221E009139 and MOST 1052221E009159. This work was also supported in part by the “Aim for the Top University Plan” of Biomedical Electronics Translational Research Center in National Chiao Tung University and Ministry of Education, Taiwan.
Author information
Authors and Affiliations
Corresponding author
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Chen, YC., Hsiao, TC. Instantaneous phase difference analysis between thoracic and abdominal movement signals based on complementary ensemble empirical mode decomposition. BioMed Eng OnLine 15, 112 (2016). https://doi.org/10.1186/s1293801602337
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1293801602337