Recycling cardiogenic artifacts in impedance pneumography

RECYCLING CARDIOGENIC AR TIF A CTS IN IMPED ANCE PNEUMOGRAPHY Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 A B S T R AC T . Purpose: Biomedical sensors often exhibit cardiogenic ar- tifacts which, while distorting the signal of interest, carry useful hemo- dynamic information. W e propose an algorithm to remove and extract hemodynamic information from these cardiogenic artifacts. Methods: W e apply a nonlinear time-frequency analysis technique, the de-shape synchrosqueezing transform (dsSST), to adaptiv ely isolate the high- and low-frequenc y components of a single-channel signal. W e demonstrate this technique’ s effecti veness by removing and 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 dsSST is suitable for generating useful hemodynamic information from the cardio- genic artifact in a single-channel IP . W e propose that the usefulness of the dsSST as a recycling tool extends to other biomedical sensors exhibiting cardiogenic artifacts. 1. I N T RO D U C T I O N The popularity and quality of physiological measurement de vices hav e gro wn significantly in recent years [21]. These devices see increasing ap- plicability because of fields like mobile health [14, 23]. In critical clinical scenarios like the intensi ve care unit or the operating room, we are able to af ford multiple sensors, each one optimized to monitor a specific physiologi- cal system. A typical example is the patient monitor commonly seen at the bedside. Ho we ver , this is not alw ays the case. For e xample, in an amb ulance, only a few sensors are av ailable, and for most mobile health applications, only one or two sensors are used. Even in en vironments where patient moni- tors are present, due to technical problems, we may not be able to trust the quality of the information recorded. For e xample, when the respiratory flow channel is dead, we can count on the respiratory information hidden in the 1 Department of Laboratory Medicine and Pathobiology , Uni versity of T oronto, T oronto, ON, Canada. 2 Department of Mathematics, Duke Uni versity , Durham, NC, USA. 3 Department of Statistical Science, Duke Uni versity , Durham, NC, USA. 1 2 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 photoplethysmogram (PPG). Therefore, e very bit of information is precious in any scenario. W e 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 e xample, the respiratory infor- mation exists as the amplitude modulation of the electrocardiogram (ECG), while it exists as the lo w-frequency component of the PPG signal. Further- more, 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 sour ce separation (scBSS) problem. Under some conditions, the problem can be easily solved by applying a bandpass filter . When e xtra channels are a vailable, an adapti ve filtering scheme can be helpful. A more sophisticated approach is required when only one channel is av ailable. 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 irrelev ant 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 e ven turn the artifact into useful physiological information, and this is our target in this paper . T ypical examples include the cardiogenic artifact in the respiratory signal or the respiratory artifact in the PPG signal. W e focus on the scBSS problem when there is only a single recorded channel a vailable. W e use a nonlinear-type time-frequency filtering strate gy to extract as much information as possible from a single-channel recording. 1.1. Cardiogenic artifacts. T ake the commonly seen cardiogenic artifact in biomedical measurements as an example of maximizing information quantity . Cardiogenic artifacts come from the blood v olume mov ement in the thorax [1] and may mask the physiological information in which we hav e interest. Examples of measurements in which cardiogenic artifacts may be observed include the impedance pneumography (IP) [5, 12, 25], the respiratory inductiv e plethysmography (RIP) [30], thoracic or abdominal mov ements recorded via piezo-electric band [10], the capnogram [20], the electroencephalogram [24] and the esophageal pressure signal [18]. Pub- lished algorithms intended for clinical deployment ha ve focused on removing cardiogenic artifacts using techniques such as digital filtering [15, 26], adap- ti ve filtering [12, 29], template subtraction via the electrocardiogram [18, 19], RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 3 the synchrosqueezing transform [11], or blind source separation [24]. Re- moving the cardiogenic artifact enhances and clarifies the information re- layed to medical professionals; it accurately estimates the main information present in the recording and allo ws for its prudent utilization. Filtering out the cardiogenic artifact to enhance the main signal is a strategy suitable for clinical deployment because of the v ariety of resources av ailable in such settings; hemodynamic information may be measured using standard methods such as the electrocardiogram. Howe ver , when only a few channels are av ailable, the cardiogenic artif act is valuable because it provides additional physiological information. Utilization of the cardiogenic artifact has been suggested and applied se veral times in the past. For example, in a thoracocardiography [2, 17], cardiac output is estimated by first separating the cardiac wav eforms from the dominant respiratory signal. In [10], the authors suggest using the cardiogenic artifact in thoracic and abdominal mov ements via piezo-electric band to detect central sleep apnea. 1.2. Our contributions. In this paper , we apply two recently dev eloped nonlinear-type time-frequency (TF) analysis techniques, namely the syn- chrosqueezing transform (SST) and the de-shape SST (dsSST), to achie ve the scBSS task and hence maximize the quantity of physiological informa- tion which may be extracted from an individual signal. W e demonstrate ho w 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 under going bronco-endoscopies for the sake of diagnos- ing pulmonary disease. The MA TLAB code is made publicly a vailable so that our methods may be reproduced and applied to other signals for other purposes. 2. M O D E L I N G T H E C A R D I O G E N I C A RT I FA C T W e model the cardiogenic artifact by the wav e-shape function and the adaptive non-harmonic model [9, 27]. The adapti ve non-harmonic model is moti vated by a need to describe oscillatory physiological phenomena such as respiration and cardiac acti vity . There are se veral 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 . Moreov er , one physiological signal may not contain only one oscillatory component. T ake the IP into account, whose primary purpose is to measure respiration. The amplitude becomes larger and the c ycle 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 acti vity is also recorded as 4 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 another oscillatory component. Since this recorded cardiac acti vity is not the main purpose of the IP , it is commonly vie wed as nuisance and called the cardiogenic artifact. Due to heart rate v ariability , the cardiac acti vity 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 reg ard to the time-varying nature of their morphologies. In [27], the follo wing adaptive non-harmonic model is proposed to capture signals of this kind: (1) f ( t ) = L ∑ l = 1 a l ( t ) s l ( φ l ( t )) + Φ ( t ) , where f is the recorded signal, L describes the number of oscillatory compo- nents, a l ( t ) describes the amplitude of the l th oscillatory component (which is assumed to be positive), s l describes the oscillatory pattern of the l th component (which is assumed to be 1 -periodic), φ l describes the phase of the l th component (which is assumed to be monotonically increasing so that 1 / φ 0 l ( t ) describes the period of the cycle at time t ), and Φ is assumed to be the noise that contaminates the signal. W e call a l ( t ) the amplitude modulation (AM) of the l th component. Note that the in verse of the period of a cycle is in general understood as the frequency , so φ 0 l ( t ) is called the instantaneous fr equency (IF) of the l th component. Finally , s l is called the wave-shape function , which captures the oscillatory pattern of the l th signal. The model is called “non-harmonic” since the oscillation is non-sinusoidal. For the IP , L = 2 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 se veral other works, for example, [7, 8, 28]. Mild assumptions are needed for the adaptiv e non-harmonic model to be well-behav ed. The model in (1) is further expanded in [9] to capture the time-v arying 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 av oid distracting the focus. W e refer readers with interest to [9] for technical details. 3. T H E D E - S H A P E S Y N C H R O S Q U E E Z I N G T R A N S F O R M Moti vated by biomedical signals and the adaptiv e non-harmonic model (1), the dsSST is designed in [9] as a signal processing tool to extract information from f defined in (1). W e summarize the dsSST here. Again, to a void con voluting the main idea of this paper with technical details, the follo wing description will be kept light, and the interested reader can find RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 5 F I G U R E 1 . An illustration of modeling the IP by the adapti ve non-harmonic model. W e plot (from top to bottom) the respi- ratory component, the hemodynamic component, generated noise, and the recorded IP . W e highlight the respiratory wa ve- shape function ( s 1 ) and the cardiac wa ve-shape function ( s 2 ) on the right-hand side. rigorous mathematical details in [3, 4, 9]. Readers familiar with speech processing will recognize a cepstral device which is also employed when calculating mel-frequency cepstral coef ficients for audio signals. [16] W e start by discussing how the well-kno wn short-time Fourier transform (STFT) behav es on a giv en signal f satisfying (1). In general, the STFT aims to capture ho w the signal oscillates at different times by truncating the signal into pieces. Specifically , with a chosen window function h , such as a Gaussian function centered at the origin, the STFT is defined as (2) V ( h ) f ( t , ξ ) = Z f ( τ ) h ( τ − t ) e − i 2 π ξ ( τ − t ) d τ , where t ∈ R indicates time and ξ ∈ R indicates frequency . W e call | V ( h ) f ( t , · ) | 2 the spectrogram of the signal f at time t , since it represents the po wer spec- trum of the truncated signal f ( · ) h ( · − t ) around t . The moti v ation of the dsSST is the following fact. Due to the non- sinusoidal oscillation, at each time t , we will see dominant peaks around the fundamental frequency and its multiples in | V ( h ) f ( t , · ) | 2 . When there are mul- tiple oscillatory components in f , multiples of the fundamental frequency of one component may mask the fundamental frequency of another component. 6 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 Moreov er , it is often difficult to determine whether a ridge in the spectrogram corresponds to a multiple or a fundamental tone. W e 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 frequenc y 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 wav e-shape functions are non- sinusoidal. T o achiev e this goal, note that when the wav e-shape function of a component is non-sinusoidal, there is an oscillatory pattern in the po wer spectrum | V ( h ) f ( t , · ) | 2 whose cycles show at the fundamental frequency and its multiples. A nai ve idea which follo ws is that the frequency of this oscillation provides information about the period of the signal f . W e take into account the old cepstrum idea in signal processing [13] and deri ve the short-time cepstral transform (STCT) for use in our dynamical setup. The STCT is defined by (3) C ( h , γ ) f ( t , q ) : = Z | V ( h ) f ( t , ξ ) | γ e − i 2 π q ξ d ξ , where γ > 0 is suf ficiently small and q ∈ R is called the quefrenc y (its unit is seconds or any feasible unit in the time domain). The reason for taking the γ th po wer of | V ( h ) f ( t , ξ ) | is delicate. While | V ( h ) f ( t , ξ ) | 2 does oscillate, the amplitude of this oscillation changes from one cycle to another . T o remov e the influence of this amplitude modulation, we could take the natural logarithm of | V ( h ) f ( t , ξ ) | so that the amplitude modulation is decoupled as a “lo w-frequency component. ” (The remaining signal consists of cycles of constant amplitude.) Howe ver , taking the natural logarithm might be unstable numerically , so we use the approximation | V ( h ) f ( t , ξ ) | γ , called the “soft logarithm. ” Systematic exploration in this regard can be found in [9]. Ultimately , we obtain the fundamental period and its multiples in C ( h , γ ) f ( t , · ) . The nonlinear mask for the spectrogram is gi ven by (4) U ( h , γ ) f ( t , ξ ) : = C ( h , γ ) f ( t , 1 / ξ ) , where ξ > 0 is gi ven in Hz . This nonlinear mask is designed by taking the fact that the tw o main quantities describing oscillation, namely period and frequency , are in verse to one another . Since C ( h , γ ) f ( t , · ) captures the funda- mental period and its multiples at time t , U ( h , γ ) f ( t , · ) captures the fundamental RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 7 frequency and its di visions. Since the common information in V ( h ) f ( t , · ) and U ( h , γ ) f ( t , · ) is no w the fundamental frequency , we can remove multiples from the STFT by (5) W ( h , γ ) f ( t , ξ ) : = V ( h ) f ( t , ξ ) U ( h , γ ) f ( t , ξ ) , where ξ > 0 is interpreted as frequency . The final TF representation is | W ( h , γ ) f ( t , ξ ) | 2 , which is a nonlinearly filtered spectrogram. This step could be vie wed as applying a nonlinear filter to the signal to remove the influence of the wa ve-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 [3, 4], namely (6) S W ( h , γ , υ ) f ( t , ξ ) = Z | W ( h , γ ) f ( t , η ) | 2 δ | ξ − Ω ( h , υ ) f ( t , η ) | d η , where ξ ≥ 0 and δ means Dirac measure, and the r eassignment rule Ω ( h , υ ) f is determined by (7) Ω ( h , υ ) f ( t , ξ ) : =      − ℑ V ( D h ) f ( t , ξ ) 2 π V ( h ) f ( t , ξ ) when | V ( h ) f ( t , ξ ) | > υ − ∞ when | V ( h ) f ( t , ξ ) | ≤ υ . Here, D h ( t ) is the deriv ati ve of the chosen window function h , ℑ means the imaginary part, and υ > 0 gi ves a threshold so as to a void instability in the computation when | V ( h ) f ( t , ξ ) | is small. S W ( h , γ , υ ) f is what we call the dsSST , and S W ( h , γ , υ ) f ( t , · ) provides a sharpened spectrogram of the oscillatory signal at time t that is free of the influence of the wa ve-shape function. 3.1. Discrete Case. The continuous signal f is uniformly sampled over a discrete set of time points with sampling interv al ∆ t > 0 . The sampling rate is hence f s = ∆ − 1 t . Suppose the recording starts at time t = 0 . Write the uniformly sampled signal as a column vector f ∈ R N , where N is the number of samples and the ` -th entry of f is f ( ` ∆ t ) . W e index our v ectors and matrices beginning with 1 . Choose a discrete windo w function h ∈ R 2 K + 1 , a discretization of a chosen windo w h , which satisfies h ( K + 1 ) = 1 . W e say that 2 K + 1 is the window length. Write h 0 ∈ R 2 K + 1 for the discretization of the deri vati ve of the windo w function. For e xample, a discrete Gaussian windo w (and its deriv ativ e) with standard deviation σ > 0 sampled ov er the 8 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 interv al [ − 0 . 5 , 0 . 5 ] at a sampling interv al of 1 2 K could be h ( k ) = e − ( k − 1 2 K − 0 . 5 ) 2 2 σ 2 (8) h 0 ( k ) = −  k − 1 2 K − 0 . 5  h ( k ) σ 2 , (9) where k = 1 , . . . , 2 K + 1 . Introduce a parameter M so that 2 M is the chosen number of pixels in the frequency axis of our time-frequency representation. Ev aluate the STFT of f , a matrix V f ∈ C N × 2 M whose entries are (10) V f ( n , m ) = 2 K + 1 ∑ k = 1 f ( n + k − K − 1 ) h ( k ) e − i 2 π ( k − 1 )( m − 1 ) 2 M , where f ( l ) : = 0 when l < 1 or l > N , n = 1 , . . . , N is the time index. and m = 1 , . . . , 2 M is the frequency index. Next, e valuate the STCT of f . The STCT is a matrix C f ∈ C N × 2 M whose entries are (11) C f ( n , m 0 ) = 2 M ∑ m = 1 | V f ( n , m ) | γ e − i 2 π ( m − 1 )( m 0 − 1 ) 2 M , where γ > 0 is the chosen power parameter and m 0 = 1 , . . . , 2 M is the que- frency index. W e crop C f and consider only the first M + 1 columns asso- ciated with the positi ve quefrenc y axis. The in verted STCT of f is a matrix U f ∈ R N × ( M + 1 ) . For each time index n , consider the function g n : [ 0 , ∞ ] → R whose kno wn values are g n  1 m − 1  = C f ( n , m ) m = 1 , ..., M + 1 . (12) The entries of the in verted STCT are calculated by interpolation: U f ( n , m ) = g n  m − 1 2 M  . (13) The de-shape STFT of f , a matrix W f ∈ C N × ( M + 1 ) , is gi ven by the pointwise product (14) W f ( n , m ) = V f ( n , m ) U f ( n , m ) , where n = 1 , . . . , N and m = 1 , . . . , M + 1 . T o sharpen W f , we first calculate (15) V 0 f ( n , m ) : = 2 K + 1 ∑ k = 1 f ( n + k − K − 1 ) h 0 ( k ) e − i 2 π ( k − 1 )( m − 1 ) 2 M . W e choose a threshold υ > 0 and calculate the reassignment operator Ω Ω Ω υ f ( n , m ) = ( − ℑ V 0 f ( n , m ) V f ( n , m ) N 2 π ( 2 K + 1 ) when | V f ( n , m ) | > υ − ∞ when | V f ( n , m ) | ≤ υ . (16) RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 9 The dsSST of f , a matrix S W υ f ∈ C N × ( M + 1 ) , is finally gi ven by the formula S W υ f ( n , m ) = ∑ l ; m = l − Ω Ω Ω υ f ( n , m ) W f ( n , l ) . (17) W e could also use the reassignment operator to obtain the SST of f , denoted as S V υ f ∈ C N × ( M + 1 ) , using the formula S V υ f ( n , m ) = ∑ l ; m = l − Ω Ω Ω υ f ( n , m ) V f ( n , l ) . (18) 3.2. Numerical implementation. W e provide numerical details for both the discrete dsSST and the discrete SST . First, to decrease the computational load, when suitable, we suggest to ev aluate the STFT of f ∈ R N at only a sub- set of its N sampling points. For example, choosing a “hop” of s ∈ N samples means we calculate only those rows V f ( n , · ) for which n = s , 2 s , . . . , b N s c s . Here, b x c means the largest inte ger smaller than x > 0 . Second, we take the real part of the STCT and lift negati ve entries to zero. Third, we select a least upper bound u > 0 on the fundamental frequency of all components expected in f . This bound allows us to eliminate en velope information from the STCT below the quefrenc y 1 / u . When performing reassignment in both the dsSST and the SST , the threshold υ is selected by taking the p th quantile of the v alues  | V f ( n , m ) | : n = s , 2 s , ...,  N s  s , m = 1 , ..., M + 1  . (19) Selecting p to be a higher quantile results in the remov al of noise. Since the STFT may e xhibit large changes in norm o ver time, the threshold is deter - mined as a function of each time point m = 1 , ..., M + 1 . When calculating the in verted STCT , we use shape-preserving piecewise cubic interpolation. Finally , we suggest to remove the lo w-frequency information below l > 0 from all of our time-frequency representations because it tends to domi- nate the visualization. For example, in this paper , based on physiological kno wledge, we remove information belo w the frequency 0 . 1 Hz . MA TLAB code for performing all of the abov e-mentioned time-frequency analysis techniques is av ailable at https://github.com/jrvmalik . 4. R E C Y C L I N G A L G O R I T H M The recycling algorithm is demonstrated by e xtracting hemodynamic in- formation 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 pa- rameters for the dsSST are gi ven in T able 1. Parameters were chosen in an ad hoc fashion. The IHR is deemed to be the dominant curv e in the dsSST 10 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 exceeding 50 beats per minute ( bpm ). This dominant curv e is extracted by performing the follo wing optimization [3, (15)]: curve = arg max c ∈ N b N s c    b N s c ∑ k = 1 S W υ f ( ks , c ( k )) − λ b N s c ∑ k = 2 [ c ( k ) − c ( k − 1 )] 2    . (20) The dominant curv e curve is a function of a subset of the sampling points of f and assumes positi ve inte ger values between 1 and M + 1 . The second term in the objecti ve functional penalizes the curve for performing large jumps in the frequency axis between two consecuti ve time points; that is, it controls the regularity of the extracted heart rate. Our implementation relies on a choice of λ = 1 . Applying a linear transformation using the sampling rate f s of f yields an estimate d IHR for the IHR in beats-per-minute (bpm): (21) d IHR ( k ) : = 60 × f s ( curve ( k ) − 1 ) 2 M , where k = 1 , 2 , . . . ,  N s  . In our numerical experiments, the search for the dominant curve includes the follo wing additional steps. First, perform time- av eraging of the deshape matrix to obtain a power spectrum-lik e signal P f ( m ) : = b N s c ∑ k = 1 | S W υ f ( ks , m ) | 2 , (22) where m = 1 , ..., M + 1. Next, detect the maximum τ = arg max  P f ( m ) : 50 ≤ 60 × f s ( m − 1 ) 2 M  . (23) The figure τ plays the role of an initial heart rate range estimator . Finally , optimize (20) while ensuring that the optimizer c ∈ N b N s c satisfies the con- straint 50 ≤ 60 × f s ( c ( k ) − 1 ) 2 M ≤ 60 × f s ( τ − 1 ) 2 M + 30 (24) for all k = 1 , ...,  N s  . In other words, the extracted curve may exceed the initial heart rate estimate by at most 30 bpm . This extra step reduces the risk of detecting a non-vanishing multiple of the cardiac frequency when no human input is pro vided into the curve e xtraction process. Indeed, as is discussed in [9], harmonics may sometimes not be completely eliminated, or artifacts may arise due to issues like a missing fundamental tone, stacked harmonics, or numerical instability when in verting the discrete STCT . W e ev aluate the ef fectiv eness of the IHR estimation in the following way . Using the accompanying electrocardiogram, the ground-truth IHR, which is RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 11 T A B L E 1 . Parameters for the SST , the dsSST , and curve e xtraction Parameter V alue (IHR extraction) V alue (source separation) Input Sampling Rate f s 64 Hz 64 Hz W indow Length 2 K + 1 4001 samples 4001 samples W indow T ype Gaussian Gaussian Frequency Resolution M 15000 15000 Reassignment Quantile p 60 0 Hop s 16 samples 1 sample Power γ 0 . 3 N/A Lower Bound l 5 6 Hz 0 . 1 Hz Upper Bound u 4 . 0 Hz 2 . 0 Hz Heart Rate Search Range 50-240 bpm N/A Curve Penalty λ 1 N/A vie wed as a continuous function on R , is sampled at r i with the value 1 r i − r i − 1 , where r i is the location (in seconds) of the i th R peak in the electrocardiogram. Suppose there are N R detected R peaks. Using shape-preserving piece wise cubic interpolation ov er { r i , 1 r i − r i − 1 } N R i = 1 , we recover the ground-truth IHR series, IHR ∈ R b N s c , where the k -th entry of IHR is the ground-truth IHR at time ks f s in bpm . The following root mean square error (RMSE) metric is used to compare the estimated heart rate, d IHR , to the ground-truth heart rate, IHR . RMSE = v u u t 1  N s  b N s c ∑ k = 1 h IHR ( k ) − d IHR ( k ) i 2 . (25) Note that these values appear in units of beats-per -minute ( bpm ). In clinical practice, physicians make decisions based on the ongoing heart rate , which is an average o ver a maximum delay of 10 seconds of the IHR, according to the ANSI/AAMI standard [6]. For this clinical need, we consider IHR 10 , which is the signal obtained after applying a 10-second bi-directional mo ving av erage filter to IHR , and we ev aluate the 10-second RMSE, defined as RMSE10 = v u u t 1  N s  b N s c ∑ k = 1 h IHR 10 ( k ) − d IHR ( k ) i 2 . (26) 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 vie wed 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 pre viously in se veral places, for example, [10]. The respiratory signal r ∈ R b N s c is reconstructed from the 12 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 SST of f as (27) r ( k ) = Re      ∑ m ; 0 . 1 ≤ ( m − 1 ) f s 2 M ≤ d IHR ( k ) 60 − 0 . 2 S V υ f ( ks , m ) M      , where k = 1 , 2 , . . . ,  N s  and Re denotes taking the real part. Removing the respiratory signal r from the IP f yields an approximation of the cardiac signal. Since the cardiac signal may contain unwanted lo w-frequency infor - mation coming from the spectral linkage of the respiratory signal and noise, we apply a bi-directional highpass third-order Butterworth filter with cutof f frequency 0 . 5 Hz as a post-processing step, and we add the low-frequenc y component back to the separated respiratory signal. 4.1. Comparison with other approaches. T o achiev e 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 achiev ed a reconstruction of the cardiac signal in their attempts to clarify the respiratory signal. In [15], the cardiogenic artifact is e xtracted from the thoracic impedance signal using a smoothing cubic spline filter . W e implement the same smoothing cubic spline filter with a cutof f of 0 . 5 Hz. Since the ground-truth cardiac acti vity present in the IP recording is dif ficult to assess, only a qualitati ve comparison of the tw o algorithms is provided. W e may , howe ver , augment the filtering in [15] with one of two spectral methods to achiev e an estimate for the IHR which riv als the dsSST approach. T o show that the advantages of the dsSST include a cancellation of respiratory harmonics exceeding the threshold 0 . 5 Hz, we do the following. First, remov e the low-frequenc y component of the IP and obtain the signal f hi . Second, estimate the IHR by either: extracting the dominant curve from the SST S V υ f hi in the range 50 - 240 bpm, or determining the highest peak in the po wer spectrum (DFT) of f hi in the range 50 - 240 bpm. The curve extraction method applied to the SST matrix is the same as the curv e extraction method applied to the dsSST matrix. T o e valuate whether the proposed recycling algorithm surpasses these two additional estimates for the IHR, we apply a one-sided W ilcoxon signed-rank test under the null hypothesis that the dif ference between the pairs follo ws a symmetric distribution around zero. W e consider the significance lev el of 0 . 05. 5. R E S U LT W e retrospecti vely analyze a data set containing IP signals from a prospec- ti ve, randomized study conducted at the tertiary medical center , Chang Gung Memorial Hospital, in Linkou, New T aipei, T aiwan. The study protocol RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 13 F I G U R E 2 . W e plot short segments of tw o impedance pneu- mographies. The first segment is out-of-range and was re- jected; the second signal was judged to be of high quality and was included in our analysis. was appro ved by the Chang Gung Medical Foundation Institutional Revie w Board (No.104-0872C). All of the enrolled patients provided written in- formed consent. IP signals were recorded using the Philips Patient Monitor MP60 from subjects under going procedural sedation during bronchoscopy examinations for the sake of pulmonary disease diagnosis. The IP signals were recorded at a sampling rate of 64 Hz . Thirty-fiv e subjects were enrolled in the study . The 35 recordings were e xamined for over -saturation (i.e. out-of-range). W e attempted to select a contiguous, in-range, 10 -minute segment from each recording. Elev en subjects were excluded from our analysis because no such segment existed. Fi ve subjects were excluded because the corresponding electrocardiogram was of lo w quality . Nineteen subjects were left for our analysis. In Fig. 2, we show two impedance pneumographies; the first signal is a typical o ver -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 10 -minute IP is displayed in Fig. 3. In red, we sho w the ground-truth IHR obtained from the corresponding electrocardio- gram. The estimate d IHR 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 T able 2, we show the RMSE and RMSE10 v alues for the 19 se gments selected for analysis. The RMSE and RMSE10 values for all subjects were 2 . 29 ± 0 . 74 bpm (mean ± standard de viation) and 1 . 62 ± 0 . 78 bpm, respecti vely . Discrepancies between the ground-truth IHR and the IHR estimated from the dsSST can be explained by visually inspecting both signals. In Fig. 4, 14 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 F I G U R E 3 . W e plot the squared modulus of the dsSST of a 3 -minute IP . In red, we sho w the ground-truth IHR obtained from the corresponding electrocardiogram. The estimate for the IHR is obtained by e xtracting the dominant curve in the dsSST . In blue, we superimpose the original IP . In gray , we superimpose the corresponding electrocardiogram. T A B L E 2 . Heart rate estimation error by different methods dsSST SST HPF RMSE RMSE10 RMSE RMSE10 RMSE RMSE10 Subject 1 2.26 1.48 2.26 1.48 12.62 12.49 Subject 2 1.69 0.70 1.69 0.70 1.91 1.09 Subject 3 3.11 1.10 3.59 2.05 5.56 4.70 Subject 4 2.41 1.42 2.61 1.74 3.11 2.39 Subject 5 2.59 2.19 42.81 42.76 45.77 45.72 Subject 6 3.13 3.54 1.25 1.06 4.18 4.01 Subject 7 2.03 1.43 2.01 1.40 2.55 2.07 Subject 8 2.27 1.46 2.38 1.63 3.22 2.71 Subject 9 1.78 1.36 2.41 2.11 2.72 2.44 Subject 10 2.02 1.31 2.04 1.32 8.39 8.22 Subject 11 3.77 2.57 3.82 2.63 4.13 3.03 Subject 12 1.16 0.75 1.16 0.76 1.88 1.63 Subject 13 1.77 0.94 1.76 0.93 4.10 3.79 Subject 14 3.73 3.41 3.10 2.80 26.88 26.86 Subject 15 2.17 1.70 2.18 1.71 2.42 1.97 Subject 16 2.20 1.54 3.25 2.79 4.66 4.30 Subject 17 2.82 1.68 2.83 1.67 3.92 3.17 Subject 18 1.45 1.26 1.46 1.26 62.41 62.40 Subject 19 1.14 0.87 1.11 0.90 2.14 1.95 Mean 2.29 1.62 4.41 3.77 10.66 10.26 Standard Dev . 0.74 0.78 9.08 9.21 16.14 16.31 RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 15 F I G U R E 4 . W e visually examine the IHR estimate af forded by the dominant curve in the dsSST . W e sho w the ground- truth IHR determined from the electrocardiogram in red, the estimated IHR in blue, and the mean-filtered ground-truth IHR in green. we sho w 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 interv al series. In Fig. 5, we sho w 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 sho w a segment of the cardiac signal in red, aligned with the original electrocardiogram in gray . Note the correspondence between the electrocardiogram’ s wav e-shape and the extracted cardiac signal’ s wav e-shape. In T able 2, we show the heart rate estimation error for the two additional methods described in 4.1. These methods were deriv ed from the existing lo wpass filter approach to clarifying the respiratory signal [15]. The heading SST indicates that the dominant curve was e xtracted from the time-frequency representation S V υ f hi , and the heading HPF indicates that the IHR was es- timated by finding a local maximum in the po wer spectrum of the filtered signal f hi . It is clear that by using the SST , we can obtain a reasonable estimate for heart rate information in the majority of cases. The deshape step improves the result in the most difficult cases. Finally , the quality of the HPF estimate is not as good as the SST or dsSST estimates. For e xam- ple, subject 1 possessed a heart rate which w as steadily increasing over the duration of the recording, resulting in high RMSE and RMSE10 v alues for the HPF method. The po wer spectrum of subject 11’ s filtered IP signal did not decay quickly enough, resulting in an IHR estimate close to 50 bpm for the HPF method. The power spectrum of subject 14’ s filtered IP signal 16 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 F I G U R E 5 . W e show the separation of the respiratory and cardiac signals. At the top of each plot, we show the respira- tory 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 . contained strong cardiac harmonics which o verpo wered the fundamental car- diac frequency . The synchrosqueezed spectrogram of subject 5 contained a dominant respiratory harmonic, resulting in a poor IHR estimate for both the HPF and SST estimation methods. Evidently , the HPF and SST estimation methods are not immune to detecting respiratory harmonics, while the HPF method suffers when the signal is noisy and when the heart rate changes significantly during the recording. The W ilcoxon signed-rank test showed a statistically significant improv ement in RMSE (respectiv ely RMSE10) by the proposed dsSST algorithm ov er the SST method with the p -v alue 4 . 36 × 10 − 2 (respecti vely 3 . 35 × 10 − 2 ). In comparison with the HPF method, the p -statistics for improvement in RMSE (respecti vely RMSE10) were both less than 10 − 4 . See Figure 6 for a qualitati ve 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 re gular . RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 17 F I G U R E 6 . W e compare the separation of the respiratory and cardiac signals by different approaches. The IP is shown in gray , ov er which the extracted respiratory signal by the proposed recycling algorithm is shown in blue. Belo w the IP signal, the ECG signal is shown in gray . The extracted hemodynamic signal by the proposed recycling algorithm is sho wn in red, and the extracted hemodynamic signal by the high pass filter of [15] is sho wn in green. 6. D I S C U S S I O N A N D C O N C L U S I O N W e recycle the physiological information, the cardiogenic artifact, that is commonly discarded from biomedical measurements such as the IP signal. Our algorithm retrie ves IHR information and extracts the oscillatory signal reflecting the hemodynamic acti vity . The corresponding MA TLAB code is made publicly av ailable so that our methods may be reproduced. The dsSST succeeds by masking respiratory harmonics in the spectrogram of the IP signal, allo wing the fundamental frequency of the cardiac component to be unambiguously estimated. This time-varying estimate is used to guide an adapti ve filtration of the IP signal, leading to a separation of the cardiac and respiratory signals. Algorithms which were pre viously used to remov e cardiogenic artifacts cannot lead to a rob ust IHR estimate because they fail to disregard respiratory harmonics. 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 2 . 29 ± 0 . 74 bpm and 1 . 62 ± 0 . 78 bpm , respecti vely . A comparison with other approaches confirms the superiority of the proposed recycling algorithm. This result indicates that the proposed recycling algorithm can pro vide an accurate estimate for the heart rate over a 10-second time frame. When only a limited number of sensors is av ailable, 18 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 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 de vices are properly set up. Additionally , to av oid the out-of-range issue, a signal quality index could indicate ho w much we should trust the result. W e discuss some points specific to the IP . Since the hemodynamic in- formation 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 HR V analysis might not be optimal. W e are able to account for low-frequenc y v ariability , but high-frequency variability is typically unobserved. Hence, while it could provide e xtra information about heart rate with accuracy less than 2 bpm , we should use the extracted hemodynamic information as merely an auxiliary resource when doing traditional HR V analysis. Second, although the respiratory and cardiac signals have well-separated frequencies, the multiples associated with the non-sinusoidal wav e-shape 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 non-adaptiv e to the signal. For patients with higher respiratory rate, the multiples will ha ve a stronger influence on the heart rate spectrum. The same reason explains why the SST performs worse – the existence of multiples associated with the respi- ratory signal contaminates the heart rate information. On the other hand, when multiples associated with the IP signal are not too strong, the adaptiv e time-v arying bandpass filter strategy in (27) provides a reasonable-enough recov ery of the respiratory signal, and hence the cardiac oscillation. Note that in this IP example, the possible clinical application of the recovered cardiac wav eform needs further exploration; for example, it would be in- teresting 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 wav e-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 o verlapping of spectra of the wa ve-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 ef ficient separation [22]. A systematic study in this regard will be reported in the near future. RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 19 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 an y of the growing number of measurement devices into a format which can be deliv ered to and interpreted by medical professionals. Ho wev er , 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 en vironmental information or extra channels are needed to help deal with the ine vitable noise artifacts. W e will systematically explore this potential in our future work. A C K N O W L E D G E M E N T S The authors ackno wledge Dr . Y u-Lun Lo for sharing the data. The authors declare no conflicts of interest. R E F E R E N C E S 1. B. Brown, D. Barber, A. Morice, and A. Leathard, Cardiac and r espiratory r elated electrical impedance changes in the human thorax , IEEE Trans. Biomed. Eng. 41 (1994), 729–734. 2. G. B. Bucklar , V . Kaplan, and K onrad E. Bloch, Signal pr ocessing technique for non-in vasive real-time estimation of car diac output by inductance cardiogr aphy (thora- cocar diography) , Medical and Biological Engineering and Computing 41 (2003), no. 3, 302–309. 3. Y .-C. Chen, M.-Y . Cheng, and H.-T . W u, Nonparametric and adaptive modeling of dynamic seasonality and tr end with heter oscedastic and dependent err ors , J. Roy . Stat. Soc. B 76 (2014), 651–682. 4. I. Daubechies, J. Lu, and H.-T . W u, Synchr osqueezed wavelet transforms: An empirical mode decomposition-like tool , Appl. Comput. Harmon. Anal. 30 (2011), 243–261. 5. M. Folke, L. Cernerud, M. Ekstr ¨ om, and B. H ¨ ok, Critical r eview of non-in vasive r espi- ratory monitoring in medical car e , Medical & Biological Engineering & Computing 41 (2003), no. 4, 377–83. 6. Association for the Advancement of Medical Instrumentation et al., American na- tional standar d: 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 intr a-wave fr equency modulation , Phil. T rans. R. Soc. A 374 (2016), 20150194. 8. D. Iatsenko, P . V . E. McClintock, and A. Stefanovska, Nonlinear mode decomposition: A noise-r obust, adaptive decomposition method , Physical Re view E 92 (2015), 032916– 1–032916–25. 9. C.-Y . Lin, S. Li, and H.-T . W u, W ave-shape function analysis–when cepstrum meets time-fr equency analysis , Journal of Fourier Analysis and Applications 24 (2018), no. 2, 451–505. 20 Y A O LU 1 , HA U-TIENG WU 2,3 , AND JOHN MALIK 2 10. Y .-Y . Lin, H.-T . W u, C.-A. Hsu, C.-W . W ang, 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 . W u, Y .-T . Lin, H.-P . Kuo, and T .-Y . Lin, Hypoventilation patterns during br onchoscopic sedation and their clinical r elevance based on capnographic and r espiratory impedance analysis , Journal of Clinical Monitoring and Computing In press (2019). 12. S. Luo, W .J. T ompkins, and J.G. W ebster , Cardio genic artifact cancellation in apnea monitoring , Engineering in Medicine and Biology Society , 1994. Engineering Ad- vances: 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 , F r om fr equency 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: pr omises and barriers , PLoS Medicine 13 (2016), no. 2, e1001953. 15. L. Poupard, M. Mathieu, R. Sart ` ene, and M. Goldman, Use of thoracic impedance sensors to scr een for sleep-disor der ed br eathing in patients with car diovascular disease , Physiological Measurement 29 (2008), no. 2, 255–267. 16. L. R. Rabiner and R. W . Schafer , Theory and applications of digital speech pr ocessing , vol. 64, Pearson Upper Saddle Ri ver , NJ, 2011. 17. M. A. Sackner , R. A. Hof fman, D. Stroh, and B. P . Krie ger, Thor acocar diography: P art 1: Nonin vasive measur ement of changes in str oke volume comparisons to thermodilu- tion , Chest 99 (1991), no. 3, 613 – 622. 18. T . F . Schuessler , S. B. Gottfried, P . Goldberg, R. E. K earney , and J. H T Bates, An adaptive filter to r educe car diogenic oscillations on esophageal pr essur e signals , Annals of Biomedical Engineering 26 (1998), no. 2, 260–267. 19. V . P . Sepp ¨ a, J. Hyttinen, and J. V iik, A method for suppr essing cardio genic oscillations in impedance pneumography , Ph ysiological Measurement 32 (2011), no. 3, 337–345. 20. T . C. Smith, A. Green, and P . Hutton, Recognition of car diogenic artifact in pediatric capnograms , Journal of Clinical Monitoring 10 (1994), no. 4, 270–275. 21. V . C. Storey and I.-Y . Song, Big data technologies and manag ement: What conceptual modeling can do , Data & Knowledge Engineering 108 (2017), 50 – 67. 22. L. Su and H.-T . W u, Extract fetal ECG fr om single-lead abdominal ECG by de-shape short time fourier transform and nonlocal median , Frontiers in Applied Mathematics and Statistics 2 (2017), no. 3, 2. 23. T . T amura, Y . Maeda, M. Sekine, and M. Y oshida, W earable photoplethysmogr aphic sensors – past and pr esent , Electronics 3 (2014), no. 2, 282–302. 24. J. A. Uriguen and B. Garcia-Zapirain, Ee g artifact remo val state-of-the-art and guide- lines , Journal of Neural Engineering 12 (2015), no. 3, 031001. 25. D. E. W eese-Mayer , R. T . Brouillette, A. S. Morrow , L. P . Conway , L. M. Klemka- W alden, and C. E. Hunt, Assessing validity of infant monitor alarms with event r ecor d- ing , The Journal of Pediatrics 115 (1989), no. 5 P AR T 1, 702–708. 26. A. J. W ilson, C. I. Franks, and I. L. Freeston, Methods of filtering the heart-beat artefact fr om the br eathing waveform of infants obtained by impedance pneumography , Medical and Biological Engineering and Computing 20 (1982), no. 3, 293–298. 27. H.-T . W u, Instantaneous fr equency and wave shape functions (I) , Appl. Comput. Har- mon. Anal. 35 (2013), 181–199. RECYCLING CARDIOGENIC AR TIF A CTS IN IMPEDANCE PNEUMOGRAPHY 21 28. J. Xu, H. Y ang, and I. Daubechies, Recursive diffeomorphism-based re gression for shape functions , SIAM Journal on Mathematical Analysis 50 (2018), no. 1, 5–32. 29. Y . Y asuda, A. Umezu, S. Horihata, K. Y amamoto, R. Miki, and S. K oike, Modified tho- racic impedance plethysmogr aphy to monitor sleep apnea syndr omes , Sleep Medicine 6 (2005), no. 3, 215–224. 30. Z. Zhang, J. Zheng, H. W u, W . W ang, B. W ang, and H. Liu, Development of a r espi- ratory inductive plethysmography module supporting multiple sensors for wearable systems , Sensors (Switzerland) 12 (2012), no. 10, 13167–13184.