Recycling cardiogenic artifacts in impedance pneumography

by   Yao Lu, et al.
Duke University

Purpose: We want to capture as much information as possible from biomedical sensors when only a few channels are available. Methods: We apply a nonlinear time-frequency analysis technique, the de-shape synchrosqueezing transform, to adaptively isolate the high- and low-frequency components of a single-channel signal. We demonstrate this technique's effectiveness by deriving hemodynamic information from the cardiogenic artifact in an impedance pneumography (IP). Results: The instantaneous heart rate is extracted, and the cardiac and respiratory signals are reconstructed. Conclusions: The de-shape synchrosqueezing transform is suitable for generating useful information from the cardiogenic artifact in a single-channel IP. We propose that the usefulness of the de-shape synchrosqueezing transform as a recycling tool extends to other biomedical sensors exhibiting cardiogenic artifacts or to any individual sensor which simultaneously monitors two or more distinct oscillatory phenomenae.



There are no comments yet.


page 1

page 2

page 3

page 4


When Ramanujan meets time-frequency analysis in complicated time series analysis

To handle time series with complicated oscillatory structure, we propose...

Over-the-Air Computation via Broadband Channels

Over-the-air computation (AirComp) has been recognized as a low-latency ...

Dyadic aggregated autoregressive (DASAR) model for time-frequency representation of biomedical signals

This paper introduces a new time-frequency representation method for bio...

On the Decomposition of Multivariate Nonstationary Multicomponent Signals

With their ability to handle an increased amount of information, multiva...

Inference of synchrosqueezing transform -- toward a unified statistical analysis of nonlinear-type time-frequency analysis

We provide a statistical analysis of a tool in nonlinear-type time-frequ...

On The Effect of Vibration on Shape Sensing of Continuum Manipulators Using Fiber Bragg Gratings

Fiber Bragg Grating (FBG) has shown great potential in shape and force s...

Parameter Tuning of Time-Frequency Masking Algorithms for Reverberant Artifact Removal within the Cochlear Implant Stimulus

Cochlear implant users struggle to understand speech in reverberant envi...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 low-frequency component of the PPG signal. Furthermore, the cardiac information exists as the high-frequency 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 nonlinear-type time-frequency filtering strategy to extract as much information as possible from a single-channel 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 piezo-electric 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 piezo-electric band to detect central sleep apnea.

In this paper, we apply two recently developed nonlinear-type time-frequency (TF) analysis techniques, namely the synchrosqueezing transform (SST) and the de-shape 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 single-channel IP, and we use this estimate to separate the IP into its respiratory and hemodynamic components. The algorithm is applied to nineteen patients undergoing bronco-endoscopies 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 wave-shape function and the adaptive non-harmonic model [26, 9]. The adaptive non-harmonic 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 time-dependent period and might have a time-varying amplitude and oscillatory pattern.

Clearly, these signals cannot be modeled using simple harmonic functions or without regard to the time-varying nature of their morphologies. In [26], the following adaptive non-harmonic model is proposed to capture signals of this kind:


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 wave-shape function, which captures the oscillatory pattern of the th signal. The model is called “non-harmonic” since the oscillation is non-sinusoidal. For the IP, since it contains not only the respiratory signal, but also the cardiogenic artifact. See Fig. 1 for an example.

Figure 1. An illustration of modeling the IP by the adaptive non-harmonic model. We plot (from top to bottom) the respiratory component, the hemodynamic component, generated noise, and the recorded IP. We highlight the respiratory wave-shape function () and the cardiac wave-shape function () on the right-hand side.

This model is also considered in several other works, for example, [8, 7, 27]. Mild assumptions are needed for the adaptive non-harmonic model to be well-behaved. The model in (1) is further expanded in [9] to capture the time-varying 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 de-shape synchrosqueezing transform

Motivated by biomedical signals and the adaptive non-harmonic 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 well-known short-time 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


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 non-sinusoidal 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 wave-shape functions are non-sinusoidal. To achieve this goal, note that when the wave-shape function of a component is non-sinusoidal, 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 short-time cepstral transform (STCT) for use in our dynamical setup. The STCT is defined by


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 “low-frequency 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


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


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 wave-shape 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


where and means Dirac measure, and the reassignment rule is determined by


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 wave-shape 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. Write

for 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


where . Introduce a parameter so that is the chosen number of pixels in the frequency axis of our time-frequency representation. Evaluate the STFT of , a matrix whose entries are


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


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


The entries of the inverted STCT are calculated by interpolation:


The de-shape STFT of , a matrix , is given by the pointwise product


where and . To sharpen , we first calculate


We choose a threshold and calculate the reassignment operator


The dsSST of , a matrix , is finally given by the formula


We could also use the reassignment operator to obtain the SST of , denoted as , using the formula


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 thquantile of the values


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 shape-preserving piecewise cubic interpolation. Finally, we suggest to remove the low-frequency information below from all of our time-frequency 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 above-mentioned time-frequency analysis techniques is available at

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.


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 beats-per-minute ():


where .

We evaluate the effectiveness of the IHR estimation in the following way. Using the accompanying electrocardiogram, the ground-truth 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 shape-preserving piecewise cubic interpolation over , we recover the ground-truth IHR series, , where the -th entry of is the ground-truth IHR at time in . The following root mean square error (RMSE) metric is used to compare the estimated heart rate, , to the ground-truth heart rate, .


Note that these values appear in units of beats-per-minute (). 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 10-second bi-directional moving average filter to , and we evaluate the 10-second RMSE, defined as

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
Table 1. Parameters for the SST, the dsSST, and curve extraction

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 time-varying 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


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 low-frequency information coming from the spectral linkage of the respiratory signal and noise, we apply a bi-directional highpass third-order Butterworth filter with cutoff frequency as a post-processing step, and we add the low-frequency 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 ground-truth 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 low-frequency 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 one-sided Wilcoxon signed-rank 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.104-0872C). 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 . Thirty-five subjects were enrolled in the study.

The 35 recordings were examined for over-saturation (i.e. out-of-range). We attempted to select a contiguous, in-range, -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 over-saturated 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 ground-truth 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.

Figure 2. We plot short segments of two impedance pneumographies. The first segment is out-of-range and was rejected; the second signal was judged to be of high quality and was included in our analysis.
Figure 3. We plot the squared modulus of the dsSST of a -minute IP. In red, we show the ground-truth 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.
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
Table 2. Heart rate estimation error by different methods

Discrepancies between the ground-truth IHR and the IHR estimated from the dsSST can be explained by visually inspecting both signals. In Fig. 4, we show the ground-truth IHR in red, the estimated IHR in blue, and the mean-filtered ground-truth IHR in green. Evidently, the estimated IHR coincides with the trend in the ground-truth IHR, and the error between the estimate and the ground-truth arises from high-frequency fluctuations in the beat-to-beat interval series.

Figure 4. We visually examine the IHR estimate afforded by the dominant curve in the dsSST. We show the ground-truth IHR determined from the electrocardiogram in red, the estimated IHR in blue, and the mean-filtered ground-truth IHR in green.

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 wave-shape and the extracted cardiac signal’s wave-shape.

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 low-pass filter approach to clarifying the respiratory signal [15]. The heading SST indicates that the dominant curve was extracted from the time-frequency 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 time-varying. The last observation is that if we only count on the SST but skip the de-shape step, the performance is downgraded. The Wilcoxon signed-rank 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 high-pass filter described in [15]. It is clear that the oscillatory pattern of the hemodynamic signal extracted by the high-pass filter is less regular.

Figure 5. We show the separation of the respiratory and cardiac signals. At the top of each plot, we show the respiratory signal in blue and the IP in gray. At the bottom of each plot, we show the extracted hemodynamic signal in red. In the middle of each plot, we show the simultaneously recorded electrocardiogram in gray.

Figure 6. We compare the separation of the respiratory and cardiac signals by different approaches. The IP is shown in gray, over which the extracted respiratory signal by the proposed recycling algorithm is shown in blue. Below the IP signal, the ECG signal is shown in gray. The extracted hemodynamic signal by the proposed recycling algorithm is shown in red, and the extracted hemodynamic signal by the high pass filter of [15] is shown in green.

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 10-second 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 out-of-range, 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 out-of-range 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 mean-filtered version of the ground-truth IHR signal, its suitability for HRV analysis might not be optimal. We are able to account for low-frequency variability, but high-frequency 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 well-separated frequencies, the multiples associated with the non-sinusoidal wave-shape function approximating the respiratory signal are not negligible. This is the main reason why the performance of the high-pass filter is not comparable with the proposed recycling algorithm – it is non-adaptive 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 time-varying bandpass filter strategy in 24 provides a reasonable-enough 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 wave-shape functions are complicated, we may need a more sophisticated approach for the separation. For example, to extract the fetal electrocardiogram from the trans-abdominal maternal electrocardiogram, due to the overlapping of spectra of the wave-shape functions approximating the oscillatory patterns of the fetal and maternal components, the manifold learning approach (e.g. non-local 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.


The authors acknowledge Dr. Yu-Lun Lo for sharing the data. The authors declare no conflicts of interest.


  • [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 non-invasive real-time 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 decomposition-like tool, ACHA 30 (2011), 243–261.
  • [5] M. Folke, L. Cernerud, M. Ekström, and B. Hök, Critical review of non-invasive 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,
  • [7] T. Y. Hou and Z. Shi, Extracting a shape function for a signal with intra-wave frequency modulation, Phil. Trans. R. Soc. A 374 (2016), 20150194.
  • [8] D. Iatsenko, P. V. E. McClintock, and A. Stefanovska, Nonlinear mode decomposition: A noise-robust, adaptive decomposition method, Physical Review E 92 (2015), 032916–1–032916–25.
  • [9] C.-Y. Lin, S. Li, and H.-T. Wu, Wave-shape function analysis–when cepstrum meets time-frequency 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 piezo-electric 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 sleep-disordered 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 Il-Yeol 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 single-lead abdominal ecg by de-shape 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. Garcia-Zapirain, Eeg artifact removal state-of-the-art and guidelines, Journal of Neural Engineering 12 (2015), no. 3, 031001.
  • [24] D. E. Weese-Mayer, R. T. Brouillette, A. S. Morrow, L. P. Conway, L. M. Klemka-Walden, 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 heart-beat 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 diffeomorphism-based 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.