A Singular Spectrum Analysis-Based Data-Driven Technique for the Removal of Cardiogenic Oscillations in Esophageal Pressure Signals

Objective: Assessing the respiratory and lung mechanics of the patients in intensive care units is of utmost need in order to guide the management of ventilation support. The esophageal pressure (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document}) signal is a minimally invasive measure, which portrays the mechanics of the lung and the pattern of breathing. Because of the close proximity of the lung to the beating heart inside the thoracic cavity, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signals always get contaminated with that of the oscillatory-pressure-signal of the heart, which is known as the cardiogenic oscillation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {CGO}$ \end{document}) signal. However, the area of research addressing the removal of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {CGO}$ \end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signal is still lagging behind. Methods and results: This paper presents a singular spectrum analysis-based high-efficient, adaptive and robust technique for the removal of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {CGO}$ \end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signal utilizing the inherent periodicity and morphological property of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signal. The performance of the proposed technique is tested on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signals collected from the patients admitted to the intensive care unit, cadavers, and also on synthetic \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signals. The efficiency of the proposed technique in removing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {CGO}$ \end{document} from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signal is quantified through both qualitative and quantitative measures, and the mean opinion scores of the denoised \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\boldsymbol {P}_{ \boldsymbol {eso}}$ \end{document} signal fall under the categories ‘very good’ as per the subjective measure. Conclusion and clinical impact: The proposed technique: (1) does not follow any predefined mathematical model and hence, it is data-driven, (2) is adaptive to the sampling rate, and (3) can be adapted for denoising other biomedical signals which exhibit periodic or quasi-periodic nature.


I. INTRODUCTION
Pleural pressure is the pressure generated surrounding the lung within the pleural space during respiration. The pleural pressure well reflects the mechanics of respiration and also the work done by the respiratory muscles. However, the pleural pressure is not uniformly distributed throughout the thoracic cavity, and hence, the measured pleural pressure from a single site in the pleural space may not portray the holistic-activity of the respiratory system. Moreover, measurement of the pleural pressure from multiple sites is impractical as it is very difficult to get direct access to the pleural space, and also the process entails great risk, which might result in a collapsed lung [1]. In 1949, Buytendijk showed in his Ph.D. thesis [2] that the esophageal pressure (P eso ) signal can be used as a surrogate of pleural pressure. P eso is measured using an esophageal balloon attached to a long and thin catheter, which is placed at the lower two-thirds of the intra-thoracic esophagus. The balloon is filled with an optimum volume of air [3].
The P eso signal often gets highly contaminated with the oscillatory-pressure-signal of the heart, which is known as the cardiogenic oscillation (CGO) signal. The task of removal of the CGO from P eso is very challenging, as the bandwidths of these two signals are very close. The bandwidth of the P eso signal varies from 0.17 Hz to 0.67 Hz [4], and the bandwidth of the CGO, i.e., the heart-rate signal varies from 0.8 Hz (48 heart-beats per minute) to 4 Hz (240 heart-beats per minute) [5]. Since the upper-band-limit of the P eso signal (0.67 Hz) and the lower-band-limit of the CGO signal (0.8 Hz) are very close (∼0.13 Hz), the direct use of the conventional filters such as bandpass and band-stop filters, having fixed cut-off frequencies on CGO-contaminated P eso signals might not provide a good denoising performance. Use of data-driven or adaptive filtering techniques is always better choices under such circumstances.
Schuessler et al. have proposed an adaptive filtering-based technique for the removal of CGO from P eso signal in [6]. The technique, which is proposed in [6] requires the P eso signal as well as the electrocardiogram (ECG) signal of the same subject, and a linear dynamic filter is used to obtain an artifact-free P eso signal. An enhanced version of the technique, which is proposed in [6] is reported in [7]. In [7], the performance of the technique is tested on P eso signals collected from eight patients only. Moreover, the adaptive filters designed in [6] and [7] require: (i) 1 minute to adapt with that of the heart rate, which suggests that both the techniques are not applicable on short-duration P eso signals and (ii) consecutive 10 stable and clean respiration efforts.
In [4], a modified adaptive noise cancellation (MANC) technique is proposed by Cheng et al. for denoising the P eso signal. The MANC technique utilizes the noisy P eso signal, and an airflow signal as a reference to estimate the CGO. Finally, the estimated CGO signal is subtracted from the noisy P eso signal. The performance of the MANC technique is tested on the P eso signals collected from Brown-Norway rats. The bandwidths of both the P eso and CGO signals are high for rats compared to that of human. Therefore, it is difficult to premise the performance of the MANC technique on P eso signals collected from intensive care unit patients. An enhanced version of the MANC technique is reported in [8].
A template subtraction-based technique is proposed by Grabhoff et al. in [9] for the removal of CGO from P eso signal. Here also, in [9], a reference-signal, namely electromyogram (EMG) is required to generate a template of the CGO signal. The techniqueis designed based on three main steps: detection of the R-peak-indices from EMG signal, template generation and template subtraction.
Zara et al. have recently proposed an ensemble empirical mode decomposition (EEMD)-based technique for the removal of CGO from P eso signal in [10], where a reference signal is not required. In [10], the CGO contaminated P eso signal is decomposed into a number of intrinsic mode functions (IMFs), and the IMFs which are associated with the P eso signal only, are summed up to obtain a CGO-free P eso signal. However, All the above discussed techniques [4], [6]- [9] (except [10]), require an additional signal, be it an ECG or EMG, along with P eso . Acquisition of an extra signal undoubtedly increases the complexity of the system and hinders patients' comfort. In [10], the selection of IMFs is done relying on visual inspection only, and therefore, the technique cannot be considered as a fully-automated one. Moreover, both EMD, and EEMD-based techniques are known to be high time-consuming.
Singular spectrum analysis (SSA) is a model-free and datadriven time-series-decomposition method. It decomposes a time series into three components: trend, seasonal components and noise [27]- [29]. Use of SSA-based methods have been proven to be very efficient in a wide range of applications including climatology [30], electricity consumption forecasting [31], biomedical image and signal processing [14], [14]- [33], gait parameter estimation [12]. However, the potential of denoising CGO contaminated P eso signal using a SSA-based method is not explored to date.
The motivation behind the proposed research work is to design a high-performance data-driven CGO removal technique from P eso signal addressing the hurdles, shortcomings and drawbacks of the aforementioned techniques. The main findings and the novelty of the proposed technique are: (i) it does not require an additional signal as reference, (ii) the proposed technique is fully automated and adaptive, (iii) very high CGO-removal performance, (iv) the proposed technique is applicable on both long and short duration signals, (v) it does not require an adaptation time (vi) the technique is adaptive to the sampling rate of the P eso signal, and (vii) the proposed technique is data-driven, and its performance does not depend on any pre-assumed mathematical function/model unlike the wavelet transform-based ones.
The remainder of this paper is organized as follows. The proposed singular spectrum analysis (SSA)-based data-driven P eso signal denoising technique is presented in Section II. The performance of the proposed technique is analyzed in Section III, and finally, the technique is discussed and conclusions drawn in Section IV.

II. SSA-BASED P eso SIGNAL DENOISING TECHNIQUE
First, the discrete Fourier transform (DFT) of the CGOcontaminated P eso signal is computed. In DFT-domain, the maximum peak-value, which appears at or above 0.8 Hz, is identified, and the corresponding frequency-index, which is denoted as F CGO , is considered as the fundamental-frequency of the CGO signal. Figure 1 exemplifies the operation. Next, the CGO-contaminated P eso signal (which is here denoted as P) of length L is used to form N -lagged vectors each of length M . The lagged-vectors are arranged in the form of a trajectory matrix, which is denoted as P Tra . The first column of the matrix P Tra contains P; the second column of the matrix P Tra contains P shifted by 1 samples, and likewise, the N th column contains P shifted by N − 1 number of samples.
Since the elements of the anti-diagonals of P Tra are the same, i.e., P ij = P ji (i = j), the trajectory matrix P Tra is a Henkel matrix. In the SSA operation, the selection of the proper value of the window-length N is of great importance. Precise separation of two signals from a composite one is possible using SSA technique if their frequency-difference is On the other hand, SSA fails to separate two signals from a composite signal if (1) the amplitudes of both the signals are equal, or (2) their frequency-difference is even less than 1 N Hz. A detail theoretical-explanation and mathematical-exploration of SSA method can be found in [18]. The value of N should be chosen in such a way that all the respiratory-efforts, i.e., the inspiratory and expiratory processes, which are present in P, are also to be present in each of the lagged vectors. A small value of N would not help in denoising the signal. On the other hand, a large value of N is not only redundant, but it also increases the computational burden. Most importantly, a large value of N may cause overlapping of the respiratory events, which in turn leads to a significant loss of the clinical information. In this research work, the window-length N is determined based on the inherent periodicity and morphological properties of the inhalation and exhalation events. The maximum rate of inhalation for humans, which is reported in the literature [4], is 40 breaths per minute. Thus, the window-length N is calculated as shown in Equation 2.
From Equation 2 it can be noted that (i) all the respiratoryefforts which are present in P, are guaranteed to be present in all the columns of P Tra , and (ii) none of the respiratoryefforts, which are present in the adjacent columns of the P Tra would overlap with its preceding respiratory-effort.
Next, the matrix P Tra is decomposed using SVD to represent it as a sum of rank-one biorthogonal elementary matrices. The SVD technique factorizes a matrix into the product of another three matrices: an orthogonal matrix (U), a diagonal matrix (S) and the transpose of an orthogonal matrix (V) [11]. Eigenvalues of the matrix S = P T Tra × P Tra are denoted as λ = λ 1 , λ 2 , . . . , λ N . The Eigenvalues are in decreasing order of magnitude (i.e., λ 1 ≥ λ 2 ≥ · · · ≥ λ N ). The corresponding eigenvectors of S are denoted then it is also possible to write the trajectory matrix as: and v i is the i th right eigenvector [12]. In the SVD operation, the matrix V (= [v 1 , v 2 , . . . , v L ]) is denoted as the right singular matrix, and it contains the information about the most important axis of the data. The vector v 1 reveals the direction having the most variance, and v L reveals the direction having the least variance. The singular values define the variance precisely. The strong inter-sample and inter-cycle correlations that exist in the respiratory signal help enhancing the covariance among the eigenvector of the U matrix. The SVD operation arranges the singular values in its decreasing order of magnitude, and the small singular values barely carry any signal-information. Therefore the small singular values could be discarded [13]. P can be considered as the composition of the P eso signal and the CGO. If the indices I = {i 1 , . . . i p }(p < d) of the eigenvalues, which hold the most of the P eso signal-information, are known, then the matrix corresponding to the P eso -signal-only can be written as follows.
The principal components (PCs) of P are computed by a linear combination of the matrix P Tra and the matrix of eigenvectors V .
The columns of the matrix PC are the principal components of P. Next, a matrix is created by an inverse projection of the principal components and the matrix of eigenvectors V . Finally, the reconstructed components (RCs) of the original P eso signal are computed by averaging along the anti-diagonals of this matrix [14].
is the number of elements in W i . The original P eso signal P can be obtained by summing up all the RC s .
Now, all the RCs are filtered using a zero-phase 4 th order Butterworth bandpass filter having lower and upper cut-off frequencies 0.17 Hz and 1.4 Hz, respectively, and a zero-phase 4 th order Butterworth notch filter, whose center-frequency is set to F CGO in order to remove the VOLUME 8, 2020  interference of the CGO signal. The filtered RCs are denoted as filRCs. Figure 2 shows an example of a few RCs and filRCs, and Figure 3 shows their corresponding frequency-spectras. From Figure 2 it can be observed that the low-order RCs remain almost unaltered, and the high-order RCs are almost diminished after the filtering operation, which suggests that the low-order RCs are signal dominant, and the high-order RCs are noise dominant. Figure 3 corroborates this notion.
A dynamic weight calculation-based method is used to discard the less-significant high-order filRCs. The weights of all the filRCs are calculated as follow: where W (n) is the weight of the n th filRC, N denotes the total number of filRCs, and filRC n (i) denotes the i th coefficients of the n th filRC. Finally, an optimum number of filRCs are summed-up until the cumulative weight reaches a predefined threshold value. The threshold value is denoted as Th W . The value of Th W , which is used in this research work, is explained in Section III. The summed-up filRCs is considered as the denoised P eso signal, i.e., the CGO-free P eso signal. The filRC-summing algorithm is given below.
Step 1: n = 1, rec = 0, weigh = 0 Step 2: rec = rec + filRC n Step Step 4: if weigh ≥ Th W dP eso = rec n opt = n else n = n + 1 go to Step2 end where n opt is the optimum number filRCs which is required to attain Th W .

III. PERFORMANCE OF THE PROPOSED TECHNIQUE
The performance of the proposed technique is quantified through qualitative as well as quantitative measures. The technique is implemented on MATLAB platform with a computer having 64-bit Windows 7 operating system, 16GB RAM and Intel Xeon CPU E3-1225 v3 3.20 GHz.

A. DATA ACQUISITION
The P eso signals are collected in three different settings as shown in Table 1 breathing, ii) breathing using a t-piece, and iii) passive breathing. T-piece is a T-shape tubing that can be connected to the endotracheal tube while the patient is disconnected from the ventilator allowing for oxygen supplementation normally used during the so called ''spontaneous breathing trials'' to test the patient's ability to breathe without the ventilator's assistance prior to removal of the endotracheal tube. Spontaneous breathing and breathing using a T-piece are considered to be 'active breathing' conditions, which signifies that the patients are able to recruit their respiratory muscles to provide some degree of effort during inspiration. In this dataset, the spontaneously breathing patients received pressure support ventilation (PSV). During PSV, the patients' inspiratory efforts triggered the ventilator to provide a pre-set positive pressure to supplement the patients' breath. On the other hand, when the patients were breathing using a T-piece, they were provided an external oxygen supply but were temporarily disconnected from the mechanical ventilator. The use of PSV and/or T-piece are common methods of weaning patients off mechanical ventilation by encouraging them to take on a greater proportion of the ventilatory effort. The other condition in this dataset is 'passive breathing', which means that the patients did not provide any inspiratory effort during inspiration. Under this condition, respiration is achieved solely by the mechanical ventilator. To achieve the passive breathing condition in this dataset, the patients were sedated, and a paralysis-inducing drug was administered to their respiratory muscles. The distinction between 'active' and 'passive' breathing conditions is important, because the morphology of the P eso signal changes depending on the condition. During 'active' breathing, the P eso signal will exhibit a negative deflection during inspiration. The reason for the negative deflection is because during inspiration, the respiratory muscles contract and the volume of the chest cavity increases. An increase in the chest cavity volume generates a pressure gradient and allows for air to flow into the lungs. Since pressure is inversely proportional to volume according to Boyle's Law [19], a negative deflection in the P eso signal is observed. Figure 4 illustrates an example of a P eso signal from a spontaneously breathing patient. In passive-breathing conditions, the patients' respiratory muscles do not contract during inspiration. Thus, the P eso signal will deflect upwards during inspiration as a result of mechanical insufflations and the chest-wall elastance. Figure 5 shows an example of a P eso signal from a passivelybreathing patient.  Setting #2 represents synthetic data under an 'active' breathing condition. Nine P eso signals were simulated using the IngMar Medical ASL 5000 (Active Servo Lung) High Fidelity Breathing Simulator and are normalized to ±1. The settings that are used on ASL5000 to generate the nine simulated P eso signals are shown in Table 2. The last set of data, under setting #3, was captured during cadaver studies. The P eso signals captured from cadavers takes into account the complex characteristics of the chest wall, respiratory muscles, and the pulmonary system during inspiration and expiration. The cadaver P eso signals provided another approach to testing and validating the proposed technique since no CGO is present in the signal. VOLUME 8, 2020 FIGURE 6. (A) Black→CGO-contaminated P eso signal, red→denoised P eso signals, (B) Black→Frequency spectra of CGO-contaminated P eso signal, red→Frequency spectra of denoised P eso signals.
Two types of catheters were used in the clinical study and in the cadavers: i) Adult Esophageal Balloon Catheter (Cooper Surgical, Inc., USA) consisting of 86 cm closed-end catheter with multiple perforations surrounded by a balloon of 9.5 cm length (polyethylene) [20] and ii) Nutrivent Multifunction Naso-Gastric Catheter (Sidam, Italy) consisting of 108 cm closed-end catheter with multiple perforations surrounded by a balloon of 10 cm length (polyethylene) [21].

B. QUALITATIVE MEASURE
First, all the P eso signals, which are collected in setting #1, are denoised using the proposed technique, and the performance of the technique is assessed through qualitative assessment. Semi-blind mean opinion score (MOS) test [15] of thirteen clinicians (clinicians having extensive expertise and experience; 3-years to 25-years, on P eso signal analysis) from around the globe (8 different countries) were carried out. A web-based survey-form has been created containing the CGO-contaminated P eso signals and the corresponding denoised P eso signals. Three P eso features: (i) end-expiratory pressure (P ee ), (ii) P eso -pressure swing ( P), (iii) pressurechanging-instance (t PC ), and the overall P eso denoising performance were considered for the MOS test. The evaluators were asked to quantify whether the features in the denoised P eso signal contained sufficient clinical information for monitoring and diagnosis by providing a quality rating. The quality ratings were: 1 (does not contain any diagnostic information), 2, 3, 4, and 5 (contain sufficient diagnostic information). MOS of the denoised P eso signals were calculated using the following equation.
where N eva = total number of evaluators, N feat = total number of features, Q = quality rating of the f th feature given by the e th evaluator. The MOS value of a particular P eso feature is expressed as MOS error of each feature and the overall denoised P eso signal are calculated using Equation 8 and 9, respectively.
According to the MOS error criteria, the quality of the reconstructed biosignal is considered to be (i) 'very good' if the MOS error lies in between 0 and 15%, and (ii) 'good' if it lies in between 15% and 35% [16]. Table 3 shows the MOS errors of various P eso features and overall denoised P eso signals. According to MOS error criteria the reconstructed P eso signals (overall) and also all of its features fall under the 'very good' category.  Here it is worth mentioning that the minimum heart rate, which has been encountered while denoising the 75 signals collected in setting #1 is 53 beats per minute (∼0.88 Hz).

C. QUANTITATIVE MEASURE
In doing the quantitative assessment, the CGO signals are extracted from all 75 P eso signals (which are collected in setting #1) using Equation 10.

CGO = CGO contaminated P eso signal
−denoised P eso signal (10) Next, the extracted CGO signals are added with the P eso signals that are acquired in settings #2 and #3, using Equation 11.
where, CGO = CGO max|CGO| , F is the scaling factor, which varies from 0.1 to 1, and P 2,3 eso is the P eso signal collected in settings #2 and #3.
In total, 1800 SynP eso signals are generated using Equation 11. The SynP eso signals are then denoised using the proposed technique. Next, the quality of the denoised-SynP eso signals are quantified through quantitative measures, such as percent root mean square (PRD). PRD is a widely accepted numerical measure for assessing the quality or acceptability of processed biosignals [15]. Hence, PRD is used in this research work to gauge the quality of the denoised SynP eso signal. The PRD values are calculated using Equation 12. (12) where, and SynP eso is the denoised-SynP eso signal. According to the globally considered standard, the quality of the reconstructed biosignal is considered to be (i) 'very good' if the PRD value lies in between 0 and 2%, and (ii) 'good' if the PRD value lies in between 2% and 9% [16]. Figure 9 shows the variation of PRD with F on SynP eso signal. From this figure it can be noted that the quality of the reconstructed signals fall under the category 'good' at F ≤ 0.3. Figures 10 to 12 show the denoising performance of the proposed technique on SynP eso signals at different values of F.
The variation of PRD with F and Th W is shown in Figure 13. From this figure it can be seen that the PRD values are reduced with increasing Th W for all the values of F, and the minimum PRD value is obtained at Th W = 99.99%. As the last few filRCs are dominated by the components of the CGO signal, the PRD value becomes high when all the filRCs are summed-up, i.e., at Th W = 100%. Therefore, in this research, the value of Th W is set to 99.99%.
Time-complexity is a figure of merit, which is often used to gauge the runtime of a technique. The runtime of the proposed technique as a function of F and the length of the signal is show in Figure 14. From Figure 14 it can be noted that the runtime of the proposed P eso signal denoising technique varies linearly with both the length of the input signal and F. In this proposed P eso signal denoising VOLUME 8, 2020    technique, the reconstructed components are filtered through Butterworth bandpass and notch filters. What would be the result if the bandpass and notch filters are directly applied on the CGO contaminated P eso signal, is shown in Figure 15. From this figure it can be noted that proposed SSA-based P eso  signal denoising technique performs much better than that of when a bandpass and notch filter are used alone.

IV. PERFORMANCE COMPARISON
A direct comparison of performance of the proposed P eso signal denoising technique with that of other techniques which are reported in [4] and [6]- [9] is not possible as the performance these techniques are evaluated with different settings and P eso databases. However, in [6], the simulated P eso signals were added with artificial CGO and white Gaussian noise (WGN) to have a signal-to-noise ratio of 10 dB, and then the noisy signals are denoised. The proposed technique is also tested likewise. A P 2,3 eso signal is added with CGO (at an F value of 0.5) and WGN (the input-SNR = −1.92 dB), and then the signal is denoised using the proposed technique. The performance of the proposed technique under the presence of both CGO and WGN is shown in Figure 16. Now, in order to make a fair comparison between the proposed P eso signal denoising technique and the one which is proposed in [10], the performance of [10] is evaluated on the SynP eso signals, which are generated using Equation 11. Then the performances of these two techniques are compared in terms of PRD and runtime at different values of F. The runtime of the technique [10] is evaluated in the same computational environment as mentioned before. The result is shown in Table 4.
From Table 4 it can be seen that the runtime of technique, which is proposed in [10] is much higher compared to the proposed one, and also the PRD values of the proposed P eso  signal denoising technique is smaller compared to [10] for all the values of F. Here it is worth mentioning that the run time of [10], which is mentioned in Table 4 comprises only the time required to decompose the P eso signal using EEMD technique. After decomposition, IMFs are selected based on visual inspection, i.e., manually, and then, those selected IMFs are summed-up to obtain a CGO-free P eso signal.

A. CLINICAL VALIDATION
A number of parameters that define the mechanics of the patient's respiratory system such as the elastance of the lung and chest wall, transpulmonary driving pressure, and those, which are needed to quantify the strength of the respiratory effort such as pressure-time product of the respiratory muscles (PTP) and work of breathing (WOB) are calculated based on the esophageal pressure signal [22], [23]). Values of these parameters can significantly be altered by the presence of the cardiac oscillations which modify the amplitude and slope of the esophageal pressure signal. The result can be either an overestimation or underestimation of the real magnitudes of the measured parameters, which can lead to an erroneous clinical judgment. In order to evaluate the performance of the proposed technique on the measurement of clinically relevant parameters, five tracings were selected randomly from patients during unassisted breathing on T-piece, including 425 breaths. These tracings were then FIGURE 17. Agreement between parameters of inspiratory effort measured using the noisy and denoised P eso signal. A and E: Bland-Altman plots showing agreement between the pressure-time product per breath and work of breathing per liter measured from the denoised signal (PTP/br filtered and WOB filtered ) compared to that calculated based on the noisy signal (PTP/br noise and WOB noise ). The difference between and the PTP/br filtered or WOB filtered and corresponding PTP/br noise or WOB noise is plotted against the average of the two variables. Black horizontal continuous lines represent mean biases, and dashed lines represent the upper and lower limits of agreement. On the right side of the figure, representative examples of breaths where the PTP/br filtered is higher (B), equal (C) and lower (D) than the PTP/br noise are shown. denoised using the proposed technique and three parameters namely, PTP per breath, PTP per minute, and WOB were calculated from both the noisy and denoised tracings. PTP per breath was calculated as the integral of the P eso signal from the beginning of inspiratory effort until the end of inspiratory flow limited by the chest wall recoil pressure (product of tidal volume and chest wall elastance), PTP per breath was calculated as the average PTP per breath for one tracing times the respiratory rate. WOB was calculated as the area enclosed in the P eso -volume loop divided by the tidal volume. The mean (standard deviation) differences between these parameters derived from the original and the denoised tracings are −1.0 (0.8) cm H2O.sec for the PTP per breath, −18.3 (10.6) cm H2O.sec/min for the PTP per minute and −0.1 (0.1) J/L for the WOB representing a relative difference of −15% (11%), −15% (8%), and −11% (11%), respectively. Figure 17 exemplifies the operation. VOLUME 8, 2020

V. CONCLUSION AND DISCUSSION
A high efficient, robust and data-driven cardiogenic oscillation removal technique from the esophageal-pressure signal is proposed in this research work. Though the performance evaluation metrics of the proposed techniques are rigorously analyzed in Section III, Figure 16 perhaps better convey the robustness of the technique in removing not only the cardiogenic oscillation, but also the white Gaussian types of noises. A very important factor about the proposed technique is that its performance does not depend on any predefined mathematical-model or function unlike the wavelet transform. Another advantage of the proposed technique is that the singular spectrum analysis parameters are made adaptive to the sampling rate of the signal, and therefore it is not required to change the parameters manually if the sampling rate of the signal alters. The proposed technique can also be easily adapted for denoising other biomedical signals such as the electrocardiogram and photoplethysmogram, which exhibit periodic or quasi-periodic nature.
The performance of the proposed technique is tested on 75 esophageal-pressure signals, which are collected from the intensive care unit of the St. Michael's Hospital, Toronto, Ontario, Canada, and also on 1800 signals which are generated by adding the synthetic esophageal-pressure signals with real cardiogenic-oscillation noises. Both the quantitative and qualitative distortion measure metrics show that the proposed denoising technique is robust enough to expel out the cardiogenic-oscillation noises efficiently. The main reasons for achieving such an attractive denoising performance are: (i) choosing the optimum value of the windowlength, (ii) enhanced covariance among the eigenvector, (iii) implementing the bandpass and notch filtering operations at the reconstruction-component levels, and (iv) the singular spectrum analysis technique captures the periodicity of the oscillatory modes of the esophageal-pressure signal better than that of the fixed filtering-based approaches. Figure 15 shows a comparison of the denoising performance between the proposed technique and fixed-filtering approaches. The resulting signals suggest that the proposed technique outperforms fixed-filtering techniques in removing CGO interference. Figures 6 to 8 show that the proposed technique is efficient enough to process the esophageal-pressure signals of different morphologies.
At its present setting, the proposed technique is applied on a full-length esophageal-pressure signal at once. However, the technique can also be adapted to use in real-time applications. Formation of a M × N Henkel matrix is the foremost criteria of a singular spectrum analysis-based method. In order to process an esophageal pressure signal using the proposed technique, following Equation 2, the number of columns of the Henkel matrix should be N = 1.5 × sampling rate (Hz). Hence, at least N + 1 numbers of samples are required in order to form the Henkel matrix, and in real-time application the proposed technique can be applied iteratively on N + 1 number of samples. For an example, if the sampling rate of the esophageal pressure signal is 200 Hz, then the minimum number of samples, which is required to process using the proposed technique, is (1.5 × 200) + 1 = 301. It is observed that the average processing time of 301 samples is ∼0.025 seconds. A real-time implementation of this proposed denoising technique can help providing a precise estimation of the respiratory mechanics and breathing effort. Moreover, it can also be applied in existing monitoring devices, such as ventilators, for the clinicians to calculate these parameters correctly in real-time allowing a personalized management of ventilator settings and sedation to avoid harm.
From a clinical perspective, we showed that the denoised signal allows for calculation of clinically important parameters derived from P eso . PTP and WOB correlate with energy expenditure of the respiratory muscles [24] and allow for estimation of the risk of patient self-inflicted lung injury and myotrauma [25], [26]. The relative and absolute differences between the parameters calculated with the noisy and denoised signal are small and clinically acceptable; the precision (i.e. distribution of the difference) is also acceptable. The results show that, in average, the parameters calculated by the denoised signal are slightly lower than that of the noisy signal. However, upper limit of agreement (mean + 1.96 SD) crosses zero suggesting that for some breaths the use of our denoising technique results estimation of higher measures of inspiratory effort and for others in a lower effort compared to the noisy signal. To tease out if there is an overall slight overestimation of the real measurement of inspiratory effort by the noisy signal or slight underestimation by the denoised signal, an independent comparison with other measures of muscular activity such as energy consumption should be performed in future studies. Only one level of pressure support; 5 cm H 2 O, was tested in the research which is considered as a relatively low level of support in the clinical setting.
The result of the subjective assessment i.e., the qualitative measure provides a better insight into the quality of the denoised esophageal-pressure signals as the assessment is conducted by the field-experts. As per the mean opinion score error criteria the denoised esophageal-pressure signals (overall) and also all of its features fall under the category 'very good'.