1. Introduction
The popularity and quality of physiological measurement devices have grown significantly in recent years [20]. These devices see increasing applicability because of fields like mobile health [14, 22]. In critical clinical scenarios like the intensive care unit or the operating room, we are able to afford multiple sensors, each one optimized to monitor a specific physiological system. A typical example is the patient monitor commonly seen at the bedside. However, this is not always the case. For example, in an ambulance, only a few sensors are available, and for most mobile health applications, only one or two sensors are used. Even in environments where patient monitors are present, due to technical problems, we may not be able to trust the quality of the information recorded. For example, when the respiratory flow channel is dead, we can count on the respiratory information hidden in the photoplethysmogram (PPG). Therefore, every bit of information is precious in any scenario.
We want to maximize the quantity of clinical information extracted from physiological signals. The main challenge is separating information which has been mixed together in one channel. For example, the respiratory information exists as the amplitude modulation of the electrocardiogram (ECG), while it exists as the lowfrequency component of the PPG signal. Furthermore, the cardiac information exists as the highfrequency component of the respiratory flow signal, etc. In general, when this mixing of signals appears as the sum of multiple oscillatory components, the task of separating them could be understood as the single channel blind source separation (scBSS) problem. Under some conditions, the problem can be easily solved by applying a bandpass filter. When extra channels are available, an adaptive filtering scheme can be helpful. A more sophisticated approach is required when only one channel is available.
The scBSS problem is more complicated when we enter the field of mobile health because of the presence of artifacts. The term “artifact” refers to an element of the recording that is irrelevant to the information in which we have interest. Usually, artifacts are not modeled as random noise but rather as a signal with some structure. When the artifact’s structure follows a set of rules, these rules can help us remove the artifact. If the rules are physiological, we may even turn the artifact into useful physiological information, and this is our target in this paper. Typical examples include the cardiogenic artifact in the respiratory signal or the respiratory artifact in the PPG signal. We focus on the scBSS problem when there is only a single recorded channel available. We use a nonlineartype timefrequency filtering strategy to extract as much information as possible from a singlechannel recording.
1.1. Cardiogenic artifacts
Take the commonly seen cardiogenic artifact in biomedical measurements as an example of maximizing information quantity. Cardiogenic artifacts come from the blood volume movement in the thorax [1] and may mask the physiological information in which we have interest. Examples of measurements in which cardiogenic artifacts may be observed include the impedance pneumography (IP) [24, 12, 5], the respiratory inductive plethysmography (RIP) [29], thoracic or abdominal movements recorded via piezoelectric band [10], the capnogram [19], the electroencephalogram [23] and the esophageal pressure signal [17]. Published algorithms intended for clinical deployment have focused on removing cardiogenic artifacts using techniques such as digital filtering [25, 15], adaptive filtering [12, 28], template subtraction via the electrocardiogram [17, 18], the synchrosqueezing transform [11], or blind source separation [23]
. Removing the cardiogenic artifact enhances and clarifies the information relayed to medical professionals; it accurately estimates the main information present in the recording and allows for its prudent utilization.
Filtering out the cardiogenic artifact to enhance the main signal is a strategy suitable for clinical deployment because of the variety of resources available in such settings; hemodynamic information may be measured using standard methods such as the electrocardiogram. However, when only a few channels are available, the cardiogenic artifact is valuable because it provides additional physiological information. Utilization of the cardiogenic artifact has been suggested and applied several times in the past. For example, in a thoracocardiography [16, 2], cardiac output is estimated by first separating the cardiac waveforms from the dominant respiratory signal. In [10], the authors suggest using the cardiogenic artifact in thoracic and abdominal movements via piezoelectric band to detect central sleep apnea.
In this paper, we apply two recently developed nonlineartype timefrequency (TF) analysis techniques, namely the synchrosqueezing transform (SST) and the deshape SST (dsSST), to achieve the scBSS task and hence maximize the quantity of physiological information which may be extracted from an individual signal. We demonstrate how to estimate instantaneous heart rate (IHR) from the cardiogenic artifact in a singlechannel IP, and we use this estimate to separate the IP into its respiratory and hemodynamic components. The algorithm is applied to nineteen patients undergoing broncoendoscopies for the sake of diagnosing pulmonary disease. The MATLAB code is made publicly available so that our methods may be reproduced and applied to other signals for other purposes.
2. Modeling the cardiogenic artifact
We model the cardiogenic artifact by the waveshape function and the adaptive nonharmonic model [26, 9]. The adaptive nonharmonic model is motivated by a need to describe oscillatory physiological phenomena such as respiration and cardiac activity. There are several facts about oscillatory physiological signals that interest us. The amplitude of the oscillation, the cycle period, and the oscillation pattern may vary from one cycle to another. Moreover, one physiological signal may not contain only one oscillatory component. Take the IP into account, whose primary purpose is to measure respiration. The amplitude becomes larger and the cycle period becomes longer when taking a deep breath, and the oscillatory pattern might change when the inhalation and exhalation pattern changes. In addition to the respiratory oscillation in the IP, the cardiac activity is also recorded as another oscillatory component. Since this recorded cardiac activity is not the main purpose of the IP, it is commonly viewed as nuisance and called the cardiogenic artifact. Due to heart rate variability, the cardiac activity also has a timedependent period and might have a timevarying amplitude and oscillatory pattern.
Clearly, these signals cannot be modeled using simple harmonic functions or without regard to the timevarying nature of their morphologies. In [26], the following adaptive nonharmonic model is proposed to capture signals of this kind:
(1) 
where is the recorded signal, describes the number of oscillatory components, describes the amplitude of the ^{th} oscillatory component (which is assumed to be positive), describes the oscillatory pattern of the ^{th} component (which is assumed to be periodic), describes the phase of the ^{th} component (which is assumed to be monotonically increasing so that describes the period of the cycle at time ), and is assumed to be the noise that contaminates the signal. We call the amplitude modulation (AM) of the ^{th} component. Note that the inverse of the period of a cycle is in general understood as the frequency, so is called the instantaneous frequency (IF) of the ^{th} component. Finally, is called the waveshape function, which captures the oscillatory pattern of the ^{th} signal. The model is called “nonharmonic” since the oscillation is nonsinusoidal. For the IP, since it contains not only the respiratory signal, but also the cardiogenic artifact. See Fig. 1 for an example.
This model is also considered in several other works, for example, [8, 7, 27]. Mild assumptions are needed for the adaptive nonharmonic model to be wellbehaved. The model in (1) is further expanded in [9] to capture the timevarying oscillatory pattern. The generalized model described in [9] is more complicated, but its essence is the same as that in (1), and (1) is enough for the purpose of this work. Therefore, to keep the discussion simple, we satisfy ourselves with (1) to avoid distracting the focus. We refer readers with interest to [9] for technical details.
3. The deshape synchrosqueezing transform
Motivated by biomedical signals and the adaptive nonharmonic model (1), the dsSST is designed in [9] as a signal processing tool to extract information from defined in (1). We summarize the dsSST here. Again, to avoid convoluting the main idea of this paper with technical details, the following description will be kept light, and the interested reader can find rigorous mathematical details in [4, 3, 9].
We start by discussing how the wellknown shorttime Fourier transform (STFT) behaves on a given signal
satisfying (1). In general, the STFT aims to capture how the signal oscillates at different times by truncating the signal into pieces. Specifically, with a chosen window function , such as a Gaussian function centered at the origin, the STFT is defined as(2) 
where indicates time and indicates frequency. We call the spectrogram of the signal at time , since it represents the power spectrum of the truncated signal around .
The motivation of the dsSST is the following fact. Due to the nonsinusoidal oscillation, at each time , we will see dominant peaks around the fundamental frequency and its multiples in . When there are multiple oscillatory components in , multiples of the fundamental frequency of one component may mask the fundamental frequency of another component. Moreover, it is often difficult to determine whether a ridge in the spectrogram corresponds to a multiple or a fundamental tone. We thus want to find a way to filter out these multiples while preserving the fundamental frequency of each component. The dsSST contains as key ingredients two nonlinear operators for this purpose. These two operators are based on exploring the information hidden in the STFT itself, namely the symmetry structure between phase and frequency, and the phase information.
The first nonlinear operator takes into account the symmetry structure between phase and frequency in order to decouple the dynamical information in which we have interest (such as the IF of each component) from the artifacts (multiples) which arise when the waveshape functions are nonsinusoidal. To achieve this goal, note that when the waveshape function of a component is nonsinusoidal, there is an oscillatory pattern in the power spectrum whose cycles show at the fundamental frequency and its multiples. A naive idea which follows is that the frequency of this oscillation provides information about the period of the signal . We take into account the old cepstrum idea in signal processing [13] and derive the shorttime cepstral transform (STCT) for use in our dynamical setup. The STCT is defined by
(3) 
where is sufficiently small and is called the quefrency (its unit is seconds or any feasible unit in the time domain). The reason for taking the ^{th} power of is delicate. While does oscillate, the amplitude of this oscillation changes from one cycle to another. To remove the influence of this amplitude modulation, we could take the natural logarithm of so that the amplitude modulation is decoupled as a “lowfrequency component.” (The remaining signal consists of cycles of constant amplitude.) However, taking the natural logarithm might be unstable numerically, so we use the approximation , called the “soft logarithm.” Systematic exploration in this regard can be found in [9]. Ultimately, we obtain the fundamental period and its multiples in .
The nonlinear mask for the spectrogram is given by
(4) 
where is given in . This nonlinear mask is designed by taking the fact that the two main quantities describing oscillation, namely period and frequency, are inverse to one another. Since captures the fundamental period and its multiples at time , captures the fundamental frequency and its divisions. Since the common information in and is now the fundamental frequency, we can remove multiples from the STFT by
(5) 
where is interpreted as frequency. The final TF representation is
, which is a nonlinearly filtered spectrogram. This step could be viewed as applying a nonlinear filter to the signal to remove the influence of the waveshape function.
The second nonlinear operator takes the phase information in the STFT into account to further sharpen the nonlinearly filtered spectrogram. This nonlinear operator is produced by applying the synchrosqueezing transform [4, 3], namely
(6) 
where and means Dirac measure, and the reassignment rule is determined by
(7) 
Here, is the derivative of the chosen window function , means the imaginary part, and gives a threshold so as to avoid instability in the computation when is small. is what we call the dsSST, and provides a sharpened spectrogram of the oscillatory signal at time that is free of the influence of the waveshape function.
3.1. Discrete Case
The continuous signal is uniformly sampled over a discrete set of time points with sampling interval . The sampling rate is hence . Suppose the recording starts at time
. Write the uniformly sampled signal as a column vector
, where is the number of samples and the th entry of is . We index our vectors and matrices beginning with . Choose a discrete window function , a discretization of a chosen window , which satisfies . We say that is the window length. Writefor the discretization of the derivative of the window function. For example, a discrete Gaussian window (and its derivative) with standard deviation
sampled over the interval at a sampling interval of could be(8)  
(9) 
where . Introduce a parameter so that is the chosen number of pixels in the frequency axis of our timefrequency representation. Evaluate the STFT of , a matrix whose entries are
(10) 
where when or , is the time index. and is the frequency index. Next, evaluate the STCT of . The STCT is a matrix whose entries are
(11) 
where is the chosen power parameter and is the quefrency index. We crop and consider only the first columns associated with the positive quefrency axis. The inverted STCT of is a matrix . For each time index , consider the function whose known values are
(12) 
The entries of the inverted STCT are calculated by interpolation:
(13) 
The deshape STFT of , a matrix , is given by the pointwise product
(14) 
where and . To sharpen , we first calculate
(15) 
We choose a threshold and calculate the reassignment operator
(16) 
The dsSST of , a matrix , is finally given by the formula
(17) 
We could also use the reassignment operator to obtain the SST of , denoted as , using the formula
(18) 
3.2. Numerical implementation
We provide numerical details for both the discrete dsSST and the discrete SST. First, to decrease the computational load, when suitable, we suggest to evaluate the STFT of at only a subset of its sampling points. For example, choosing a “hop” of samples means we calculate only those rows for which . Here, means the largest integer smaller than . Second, we take the real part of the STCT and lift negative entries to zero. Third, we select a least upper bound on the fundamental frequency of all components expected in . This bound allows us to eliminate envelope information from the STCT below the quefrency . When performing reassignment in both the dsSST and the SST, the threshold is selected by taking the ^{th}quantile of the values
(19) 
Selecting to be a higher quantile results in the removal of noise. Since the STFT may exhibit large changes in norm over time, the threshold is determined as a function of each time point . When calculating the inverted STCT, we use shapepreserving piecewise cubic interpolation. Finally, we suggest to remove the lowfrequency information below from all of our timefrequency representations because it tends to dominate the visualization. For example, in this paper, based on physiological knowledge, we remove information below the frequency . MATLAB code for performing all of the abovementioned timefrequency analysis techniques is available at https://github.com/jrvmalik.
4. Recycling algorithm
The recycling algorithm is demonstrated by extracting hemodynamic information from the cardiogenic artifact in the IP. Our recycling algorithm uses information solely from the IP. There are two steps to our algorithm. First, the IHR information is extracted from the dsSST of the IP. The parameters for the dsSST are given in Table 1. Parameters were chosen in an ad hoc fashion. The IHR is deemed to be the dominant curve in the dsSST exceeding beats per minute (). This dominant curve is extracted by performing the following optimization.
(20) 
The dominant curve curve is a function of a subset of the sampling points of and assumes positive integer values between and . The second term in the objective functional penalizes the curve for performing large jumps in the frequency axis between two consecutive time points; that is, it controls the regularity of the extracted heart rate. Our implementation relies on a choice of
. Applying a linear transformation using the sampling rate
of yields an estimate for the IHR in beatsperminute ():(21) 
where .
We evaluate the effectiveness of the IHR estimation in the following way. Using the accompanying electrocardiogram, the groundtruth IHR, which is viewed as a continuous function on , is sampled at with the value , where is the location (in seconds) of the peak in the electrocardiogram. Suppose there are detected R peaks. Using shapepreserving piecewise cubic interpolation over , we recover the groundtruth IHR series, , where the th entry of is the groundtruth IHR at time in . The following root mean square error (RMSE) metric is used to compare the estimated heart rate, , to the groundtruth heart rate, .
(22) 
Note that these values appear in units of beatsperminute (). In clinical practice, physicians make decisions based on the ongoing heart rate, which is an average over a maximum delay of seconds of the IHR, according to the ANSI/AAMI standard [6]. For this clinical need, we consider , which is the signal obtained after applying a 10second bidirectional moving average filter to , and we evaluate the 10second RMSE, defined as
(23) 
Parameter  Value ( extraction)  Value (source separation) 
Input Sampling Rate  
Window Length  samples  samples 
Window Type  Gaussian  Gaussian 
Frequency Resolution  
Reassignment Quantile  
Hop  samples  sample 
Power  N/A  
Lower Bound  
Upper Bound  
Heart Rate Search Range    N/A 
Curve Penalty  1  N/A 
The second step in our algorithm is the separation of the cardiac and respiratory signals. The respiratory signal is isolated using an approach which may be viewed as applying an adaptive timevarying bandpass filter, which is nonlinear in nature. The theoretical foundation has been established in [4, Theorem 3.3], and the technique was used previously in several places, for example, [10]. The respiratory signal is reconstructed from the SST of as
(24) 
where and denotes taking the real part. Removing the respiratory signal from the IP yields an approximation of the cardiac signal. Since the cardiac signal may contain unwanted lowfrequency information coming from the spectral linkage of the respiratory signal and noise, we apply a bidirectional highpass thirdorder Butterworth filter with cutoff frequency as a postprocessing step, and we add the lowfrequency component back to the separated respiratory signal.
4.1. Comparison with other approaches
To achieve a fair comparison, we compare our algorithm with other existing methods. While we are the first to focus on reconstructing the cardiac signal from the IP signal, other algorithms have indirectly achieved a reconstruction of the cardiac signal in their attempts to clarify the respiratory signal. In [15], the cardiogenic artifact is extracted from the thoracic impedance signal using a smoothing cubic spline filter. We implement the same smoothing cubic spline filter with a cutoff of Hz. Since the groundtruth cardiac activity present in the IP recording is difficult to assess, only a qualitative comparison of the two algorithms is provided.
We may, however, augment the filtering in [15] with one of two spectral methods to achieve an estimate for the IHR which rivals the dsSST approach. To show that the advantages of the dsSST include a cancellation of respiratory harmonics exceeding the threshold Hz, we do the following. First, remove the lowfrequency component of the IP and obtain the signal . Second, estimate the IHR by either: extracting the dominant curve from the SST in the range  bpm, or determining the highest peak in the power spectrum (DFT) of in the range 
bpm. To evaluate whether the proposed recycling algorithm surpasses these two additional estimates for the IHR, we apply a onesided Wilcoxon signedrank test under the null hypothesis that the difference between the pairs follows a symmetric distribution around zero. We consider the significance level of
. To handle the multiple comparison issue, we consider the Bonferroni correction, and we adjust the significance level to be .5. Result
We retrospectively analyze a data set containing IP signals from a prospective, randomized study conducted at the tertiary medical center, Chang Gung Memorial Hospital, in Linkou, New Taipei, Taiwan. The study protocol was approved by the Chang Gung Medical Foundation Institutional Review Board (No.1040872C). All of the enrolled patients provided written informed consent. IP signals were recorded using the Philips Patient Monitor MP60 from subjects undergoing procedural sedation during bronchoscopy examinations for the sake of pulmonary disease diagnosis. The IP signals were recorded at a sampling rate of 64 . Thirtyfive subjects were enrolled in the study.
The 35 recordings were examined for oversaturation (i.e. outofrange). We attempted to select a contiguous, inrange, minute segment from each recording. Eleven subjects were excluded from our analysis because no such segment existed. Five subjects were excluded because the corresponding electrocardiogram was of low quality. Nineteen subjects were left for our analysis. In Fig. 2, we show two impedance pneumographies; the first signal is a typical oversaturated signal and was rejected, and the second signal was judged to be suitable and was included in our analysis.
A portion of the dsSST of a minute IP is displayed in Fig. 3. In red, we show the groundtruth IHR obtained from the corresponding electrocardiogram. The estimate for the IHR is obtained by extracting the dominant curve in the dsSST. In blue, we superimpose the original IP. In gray, we superimpose the corresponding electrocardiogram. In Table 2, we show the RMSE and RMSE10 values for the 19 segments selected for analysis. The RMSE and RMSE10 values for all subjects were (mean standard deviation) and , respectively.
deShape SST  SST  HPF  
RMSE  RMSE10  RMSE  RMSE10  RMSE  RMSE10  
Subject 1  1.83  1.01  1.86  1.05  12.62  12.49 
Subject 2  1.69  0.70  1.69  0.69  1.91  1.09 
Subject 3  3.12  1.10  3.14  1.17  5.56  4.70 
Subject 4  2.41  1.42  2.55  1.65  3.11  2.39 
Subject 5  2.63  2.25  42.72  42.67  45.77  45.72 
Subject 6  1.28  1.18  1.27  1.06  4.18  4.01 
Subject 7  2.00  1.38  2.04  1.44  2.55  2.07 
Subject 8  2.27  1.46  2.39  1.63  3.22  2.71 
Subject 9  1.44  0.92  21.11  21.08  2.72  2.44 
Subject 10  1.67  0.90  1.71  0.97  8.39  8.22 
Subject 11  3.76  2.54  3.91  2.74  40.32  40.25 
Subject 12  1.02  0.54  1.02  0.54  1.88  1.63 
Subject 13  1.77  0.94  1.77  0.95  4.10  3.79 
Subject 14  5.10  4.94  6.78  6.67  27.68  27.66 
Subject 15  2.17  1.70  2.18  1.71  2.42  1.97 
Subject 16  2.20  1.54  3.26  2.80  4.66  4.30 
Subject 17  2.82  1.68  2.93  1.83  3.92  3.17 
Subject 18  1.46  1.26  1.46  1.26  4.63  4.57 
Subject 19  1.16  0.89  1.15  0.88  2.14  1.95 
Mean  2.20  1.49  5.52  4.88  9.57  9.22 
Standard Dev.  0.96  0.95  9.79  9.99  12.88  13.02 
Discrepancies between the groundtruth IHR and the IHR estimated from the dsSST can be explained by visually inspecting both signals. In Fig. 4, we show the groundtruth IHR in red, the estimated IHR in blue, and the meanfiltered groundtruth IHR in green. Evidently, the estimated IHR coincides with the trend in the groundtruth IHR, and the error between the estimate and the groundtruth arises from highfrequency fluctuations in the beattobeat interval series.
In Fig. 5, we show the result of the second step of our algorithm: separation of the cardiac and respiratory signals. At the top of each plot, we show a segment of the clarified respiratory signal in blue overlying the original IP in gray. At the bottom of each plot, we show a segment of the cardiac signal in red, aligned with the original electrocardiogram in gray. Note the correspondence between the electrocardiogram’s waveshape and the extracted cardiac signal’s waveshape.
In Table 2, we show the heart rate estimation error for the two additional methods described in 4.1. These methods were derived from the existing lowpass filter approach to clarifying the respiratory signal [15]. The heading SST indicates that the dominant curve was extracted from the timefrequency representation , and the heading HPF indicates that the IHR was estimated by finding a local maximum in the spectrogram of the signal . It is clear that while we can obtain a reasonable estimate for heart rate information in the majority of cases, the quality is not as good as the proposed algorithm. Evidently, the two additional estimates are often distracted by strong respiratory harmonics. Moreover, the HPF method suffers when the heart rate changes significantly during the recording because its estimate is not timevarying. The last observation is that if we only count on the SST but skip the deshape step, the performance is downgraded. The Wilcoxon signedrank test showed a statistically significant improvement in RMSE (respectively RMSE10) by the proposed recycling algorithm with the value (respectively ) in comparison with both additional methods.
See Figure 6 for a qualitative comparison of the extracted hemodynamic signal by the proposed recycling algorithm and the highpass filter described in [15]. It is clear that the oscillatory pattern of the hemodynamic signal extracted by the highpass filter is less regular.
6. Discussion and Conclusion
We recycle the physiological information, the cardiogenic artifact, that is commonly discarded from biomedical measurements such as the IP signal. Our algorithm retrieves IHR information and extracts the oscillatory signal reflecting the hemodynamic activity. The corresponding MATLAB code is made publicly available so that our methods may be reproduced. In our demonstration, which was performed on real signals recorded from 19 subjects, the RMSE and RMSE10 values for the heart rate estimation stage were and , respectively. A comparison with other approaches indicates the superiority of the proposed recycling algorithm. This result indicates that the proposed recycling algorithm can provide an accurate estimate for the heart rate over a 10second time frame. When only a limited number of sensors is available, this algorithm maximizes the amount of information which can be relayed to medical professionals.
It is worth noting that the output of our algorithm should be trusted when the original signal is properly recorded. In general, the accuracy of any algorithm degrades if the signal is not properly recorded. In the case of the IP, when the signal is outofrange, there is no information recorded, and no cardiogenic artifact can be analyzed. Practitioners should ensure that their monitoring devices are properly set up. Additionally, to avoid the outofrange issue, a signal quality index could indicate how much we should trust the result.
We discuss some points specific to the IP. Since the hemodynamic information extracted by the first step of our algorithm tends to align with a meanfiltered version of the groundtruth IHR signal, its suitability for HRV analysis might not be optimal. We are able to account for lowfrequency variability, but highfrequency variability is typically unobserved. Hence, while it could provide extra information about heart rate with accuracy less than , we should use the extracted hemodynamic information as merely an auxiliary resource when doing traditional HRV analysis.
Second, although the respiratory and cardiac signals have wellseparated frequencies, the multiples associated with the nonsinusoidal waveshape function approximating the respiratory signal are not negligible. This is the main reason why the performance of the highpass filter is not comparable with the proposed recycling algorithm – it is nonadaptive to the signal. For patients with higher respiratory rate, the multiples will have a stronger influence on the heart rate spectrum. The same reason explains why the SST performs worse – the existence of multiples associated with the respiratory signal contaminates the heart rate information. On the other hand, when multiples associated with the IP signal are not too strong, the adaptive timevarying bandpass filter strategy in 24 provides a reasonableenough recovery of the respiratory signal, and hence the cardiac oscillation. Note that in this IP example, the possible clinical application of the recovered cardiac waveform needs further exploration; for example, it would be interesting to ask which hemodynamic information could be extracted by analyzing its morphology. For other biomedical signals, particularly when the frequencies of different components are close and/or the waveshape functions are complicated, we may need a more sophisticated approach for the separation. For example, to extract the fetal electrocardiogram from the transabdominal maternal electrocardiogram, due to the overlapping of spectra of the waveshape functions approximating the oscillatory patterns of the fetal and maternal components, the manifold learning approach (e.g. nonlocal Euclidean median) is needed to achieve an efficient separation [21]. A systematic study in this regard will be reported in the near future.
Finally, we mention that the proposed model and our algorithm has the potential to be applied to other oscillatory biomedical signals in the field of health monitoring, where the goal is to transform numerical information (hidden or not) from any of the growing number of measurement devices into a format which can be delivered to and interpreted by medical professionals. However, we also need to mention that the proposed model and algorithm do not take structured noise (in addition to the cardiogenic artifact) into account, which is the main challenge when dealing with mobile devices. Thus, to apply the proposed model and algorithm to mobile devices for mobile health, more environmental information or extra channels are needed to help deal with the inevitable noise artifacts. We will systematically explore this potential in our future work.
Acknowledgements
The authors acknowledge Dr. YuLun Lo for sharing the data. The authors declare no conflicts of interest.
References
 [1] B. Brown, D. Barber, A. Morice, and A. Leathard, Cardiac and respiratory related electrical impedance changes in the human thorax, IEEE Trans. Biomed. Eng. 41 (1994), 729–734.
 [2] G. B. Bucklar, V. Kaplan, and Konrad E. Bloch, Signal processing technique for noninvasive realtime estimation of cardiac output by inductance cardiography (thoracocardiography), Medical and Biological Engineering and Computing 41 (2003), no. 3, 302–309.

[3]
Y.C. Chen, M.Y. Cheng, and H.T. Wu,
Nonparametric and adaptive modeling of dynamic seasonality and trend with heteroscedastic and dependent errors
, J. Roy. Stat. Soc. B 76 (2014), 651–682.  [4] I. Daubechies, J. Lu, and H.T. Wu, Synchrosqueezed wavelet transforms: An empirical mode decompositionlike tool, ACHA 30 (2011), 243–261.
 [5] M. Folke, L. Cernerud, M. Ekström, and B. Hök, Critical review of noninvasive respiratory monitoring in medical care, Medical & Biological Engineering & Computing 41 (2003), no. 4, 377–83.
 [6] Association for the Advancement of Medical Instrumentation et al., American national standard: cardiac monitors, heart rate meters, and alarms, ANSI/AAMI EC13:2002/(R)2007, 2007, http://www.pauljbennett.com/pbennett/work/ec13/ec13.pdf.
 [7] T. Y. Hou and Z. Shi, Extracting a shape function for a signal with intrawave frequency modulation, Phil. Trans. R. Soc. A 374 (2016), 20150194.
 [8] D. Iatsenko, P. V. E. McClintock, and A. Stefanovska, Nonlinear mode decomposition: A noiserobust, adaptive decomposition method, Physical Review E 92 (2015), 032916–1–032916–25.
 [9] C.Y. Lin, S. Li, and H.T. Wu, Waveshape function analysis–when cepstrum meets timefrequency analysis, Journal of Fourier Analysis and Applications 24 (2018), no. 2, 451–505.
 [10] Y.Y. Lin, H.T. Wu, C.A. Hsu, C.W. Wang, P.C. Huang, Y.H. Huang, and Y.L. Lo, Sleep apnea detection based on thoracic and abdominal movement signals of wearable piezoelectric bands, IEEE J. Biomed. Health Inform. 21 (2017), no. 6, 1533–1545.
 [11] Y.L. Lo, H.T. Wu, Y.T. Lin, H.P. Kuo, and T.Y. Lin, Identification of the pattern of hypoventilation by respiratory signal monitor during the induction of bronchoscopic sedation, Anesthesia under review (2018).
 [12] S. Luo, W.J. Tompkins, and J.G. Webster, Cardiogenic artifact cancellation in apnea monitoring, Engineering in Medicine and Biology Society, 1994. Engineering Advances: New Opportunities for Biomedical Engineers. Proceedings of the 16th Annual International Conference of the IEEE, vol. 2, 1994, pp. 968–969.
 [13] A. V. Oppenheim and R. W. Schafer, From frequency to quefrency: A history of the cepstrum, IEEE Signal Processing Magazine 21 (2004), no. 5, 95–106.
 [14] L. Piwek, D. A. Ellis, S. Andrews, and A. Joinson, The rise of consumer health wearables: promises and barriers, PLoS Medicine 13 (2016), no. 2, e1001953.
 [15] L. Poupard, M. Mathieu, R. Sartène, and M. Goldman, Use of thoracic impedance sensors to screen for sleepdisordered breathing in patients with cardiovascular disease, Physiological Measurement 29 (2008), no. 2, 255–267.
 [16] M. A. Sackner, R. A. Hoffman, D. Stroh, and B. P. Krieger, Thoracocardiography: Part 1: Noninvasive measurement of changes in stroke volume comparisons to thermodilution, Chest 99 (1991), no. 3, 613 – 622.
 [17] T. F. Schuessler, S. B. Gottfried, P. Goldberg, R. E. Kearney, and J. H T Bates, An adaptive filter to reduce cardiogenic oscillations on esophageal pressure signals, Annals of Biomedical Engineering 26 (1998), no. 2, 260–267.
 [18] V. P. Seppä, J. Hyttinen, and J. Viik, A method for suppressing cardiogenic oscillations in impedance pneumography, Physiological Measurement 32 (2011), no. 3, 337–345.
 [19] T. C. Smith, A. Green, and P. Hutton, Recognition of cardiogenic artifact in pediatric capnograms, Journal of Clinical Monitoring 10 (1994), no. 4, 270–275.

[20]
Veda C. Storey and IlYeol Song, Big data technologies and management:
What conceptual modeling can do
, Data & Knowledge Engineering
108 (2017), 50 – 67.  [21] L. Su and H.T. Wu, Extract fetal ecg from singlelead abdominal ecg by deshape short time fourier transform and nonlocal median, Frontiers in Applied Mathematics and Statistics 2 (2017), no. 3.
 [22] T. Tamura, Y. Maeda, M. Sekine, and M. Yoshida, Wearable photoplethysmographic sensors – past and present, Electronics 3 (2014), no. 2, 282–302.
 [23] J. A. Uriguen and B. GarciaZapirain, Eeg artifact removal stateoftheart and guidelines, Journal of Neural Engineering 12 (2015), no. 3, 031001.
 [24] D. E. WeeseMayer, R. T. Brouillette, A. S. Morrow, L. P. Conway, L. M. KlemkaWalden, and C. E. Hunt, Assessing validity of infant monitor alarms with event recording, The Journal of Pediatrics 115 (1989), no. 5 PART 1, 702–708.
 [25] A. J. Wilson, C. I. Franks, and I. L. Freeston, Methods of filtering the heartbeat artefact from the breathing waveform of infants obtained by impedance pneumography, Medical and Biological Engineering and Computing 20 (1982), no. 3, 293–298.
 [26] H.T. Wu, Instantaneous frequency and wave shape functions (I), Appl. Comput. Harmon. Anal. 35 (2013), 181–199.
 [27] J. Xu, H. Yang, and I. Daubechies, Recursive diffeomorphismbased regression for shape functions, SIAM Journal on Mathematical Analysis 50 (2018), no. 1, 5–32.
 [28] Y. Yasuda, A. Umezu, S. Horihata, K. Yamamoto, R. Miki, and S. Koike, Modified thoracic impedance plethysmography to monitor sleep apnea syndromes, Sleep Medicine 6 (2005), no. 3, 215–224.
 [29] Z. Zhang, J. Zheng, H. Wu, W. Wang, B. Wang, and H. Liu, Development of a respiratory inductive plethysmography module supporting multiple sensors for wearable systems, Sensors (Switzerland) 12 (2012), no. 10, 13167–13184.
Comments
There are no comments yet.