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:
$$x\left( t \right) = \mathop \sum \limits_{i = 1}^{A} IMF_{i} \left( t \right) + r(t)$$
(1)
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:
$$\left[ {\begin{array}{*{20}c} {S_{ + } (t)} \\ {S_{  } (t)} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 1 & 1 \\ 1 & {  1} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {x(t)} \\ {w(t)} \\ \end{array} } \right]$$
(2)
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:
$$\bar{E}_{i} = \frac{1}{K}\int\limits_{t = 1}^{K} {IMF_{i}^{2} (t)_{{}} dt}$$
(3)
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:
$$\bar{T}_{i} = \frac{{\mathop \sum \nolimits_{j = 1}^{a  1} \left( {P_{i} \left[ {j + 1} \right]  P_{i} \left[ j \right]} \right)}}{a  1}$$
(4)
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:
$$\bar{f}_{i} = \frac{1}{{\overline{{T_{i} }} }}$$
(5)
Instantaneous phase difference calculation
The instantaneous phases of a main component signal are calculated using normalized direct quadrature in the following equation [33]:
$$\theta_{i}^{{}} (t) = \tan^{  1} \frac{{\sqrt {1  {\text{N}}IMF_{i}^{2} (t)} }}{{{\text{N}}IMF_{i} (t)}}$$
(6)
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:
$$\theta_{IPD} = \left {\theta_{AWM}  \theta_{TWM} } \right$$
(7)
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:
$$s_{1} (t) = A\sin \left( {\frac{2\pi }{T}  t} \right)$$
(8)
$$s_{2} (t) = A\sin \left( {\frac{2\pi }{T}t + \theta_{0} } \right)$$
(9)
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:
$$s_{3} (t) = \left\{ {\begin{array}{*{20}r} \hfill {\frac{4A}{T}t,} & \hfill \quad {0 \le (t  mT) \le \frac{T}{4}} \\ \hfill {  \frac{4A}{T}t + 2A,} & \hfill \quad {\frac{T}{4} \le (t  mT) \le \frac{3T}{4}} \\ \hfill {\frac{4A}{T}t  4A,} & \hfill \quad {\frac{3T}{4} \le (t  mT) \le T} \\ \end{array} } \right.$$
(10)
where m is the number of cycles (integers 0–59). The formula for s
_{4}(t) is presented as follows:
$$s_{4} (t) = \left\{ {\begin{array}{*{20}r} \hfill {\frac{4A}{T}t + \theta_{0} ,} & \hfill \quad {0 \le (t  mT) \le \frac{T}{4}} \\ \hfill {  \frac{4A}{T}t + 2A + \theta_{0} ,} & \hfill \quad {\frac{T}{4} \le (t  mT) \le \frac{3T}{4}} \\ \hfill {\frac{4A}{T}t  4A + \theta_{0} ,} & \hfill \quad {\frac{3T}{4} \le (t  mT) \le T} \\ \end{array} } \right.$$
(11)
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.
$$n_{1} \left( t \right) = \sigma_{1} \left( t \right) + w_{1} (t)$$
(12)
$$n_{2} \left( t \right) = \sigma_{2} \left( t \right) + w_{2} (t)$$
(13)
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.
$$\sigma_{1} \left( t \right) = {\text{Asin}}\left( {\frac{2\pi }{T}t} \right)$$
(14)
$$\sigma_{2} \left( t \right) = {\text{Asin}}\left( {\frac{2\pi }{T}t + \theta_{0} } \right)$$
(15)
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.

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).