Efficient Fetal-Maternal ECG Signal Separation from Two Channel Maternal Abdominal ECG via Diffusion-Based Channel Selection
- 1Department of Mathematics, University of Toronto, Toronto, ON, Canada
- 2Department of Obstetrics and Gynecology, University of Washington, Seattle, USA
- 3Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan
There is a need for affordable, widely deployable maternal-fetal ECG monitors to improve maternal and fetal health during pregnancy and delivery. Based on the diffusion-based channel selection, here we present the mathematical formalism and clinical validation of an algorithm capable of accurate separation of maternal and fetal ECG from a two channel signal acquired over maternal abdomen. The proposed algorithm is the first algorithm, to the best of the authors' knowledge, focusing on the fetal ECG analysis based on two channel maternal abdominal ECG signal, and we apply it to two publicly available databases, the PhysioNet non-invasive fECG database (adfecgdb) and the 2013 PhysioNet/Computing in Cardiology Challenge (CinC2013), to validate the algorithm. The state-of-the-art results are achieved when compared with other available algorithms. Particularly, the F1 score for the R peak detection achieves 99.3% for the adfecgdb and 87.93% for the CinC2013, and the mean absolute error for the estimated R peak locations is 4.53 ms for the adfecgdb and 6.21 ms for the CinC2013. The method has the potential to be applied to other fetal cardiogenic signals, including cardiac doppler signals.
Fetal electrocardiogram (ECG) and the fetal heart rate (HR) provide enormous information about fetal health. For example, the fetal distress monitoring (Jenkins, 1989) or the potential risk for fetal hypoxia detection and alert by the ST analysis monitor (Belfort et al., 2015). Moreover, from clinical studies and animal models, evidence is accumulating that perinatal brain injury originates in utero, yet no means exist to detect its onset early, reliably and with simple, widely accessible means (Anblagan et al., 2016). A harbinger of brain injury is the fetal inflammatory response (Hagberg et al., 2015). There is an urgent need for early antenatal detection of fetal inflammatory response to prevent or at least mitigate the developing perinatal brain injury. In adults and neonates, complex mathematical features of heart rate fluctuations have proven promising as early diagnostic tools (Bravi et al., 2013; Fairchild et al., 2014). For the fetal monitoring, our team addressed the challenge by developing a series of biomarkers relying on non-invasively obtainable fetal HR. Our fetal inflammatory index tracks inflammation along with the fetal plasma IL-6 temporal profile in a fetal sheep model of subclinical chorioamnionitis (Durosier et al., 2015). We also derived a set of fetal HR features that is specific to brain or gut inflammation (Liu et al., 2016). Such systemic and organ-specific tracking of inflammation via fetal HR is possible due to the brain-innate immune system communication reflected in the fetal HR fluctuations, commonly referred to as the cholinergic anti-inflammatory pathway (Fairchild et al., 2011; Olofsson et al., 2012; Garzoni et al., 2013).
In spite of its broad usefulness in the fetal health, it is fair to state that in the fetal HR monitoring realm, the technological progress has been coming more gradually. This has been not due to the plethora of studies attempting and testing various approaches, but, rather, due to the intrinsic limitations of the currently used fetal HR monitoring technology. This technology is outdated, as it deploys the traditionally set low sampling rate of heart rate or ECG signal. In animal model and human cohorts, we showed that such sampling rate is bound to miss the faster temporal fluctuations of vagal modulations of fetal HR variability and leads to inaccuracies in detection of early fetal acidemia (Durosier et al., 2014; Li et al., 2015). A sampling rate of the ECG signal around 1,000 Hz is required to capture these vagal influences and this is the commonly used sampling rate for the postnatal studies and our above-cited studies on the fetal inflammatory index.
Postnatal clinical studies are typically based on multi-lead ECG recordings which, even in newborns, and certainly in adults, poses no technical challenge to attach and record from. In fetuses, however, this is not the case. Since the fetal cardiac electric field strength is order of magnitude weaker than maternal ECG's, and the lack of clinical motivation in higher quality fetal HR data, little development had been done to focus on fetal ECG (fECG) signal in the clinical monitoring until today, except the Doppler-based fetal HR extraction techniques that dominate the market. The Doppler-based fetal HR extraction techniques, however, suffer from low fetal HR sampling rates, largely due to the auto-correlation algorithms deployed in the devices (Durosier et al., 2014). Transabdominal ECG (aECG) machines overcome this limitation by capturing the actual cardiac electric field and have returned to the market during the last decade. However, their arrival has been slower than we would have hoped. Perhaps this is in part due to the general acceptance speed of new technology in medicine (related to regulatory and safety testing as well as the specific cultures), due to the high cost for each device to upgrade a hospital's delivery unit, or, more likely, the technical limitation of the fetal ECG extraction from the aECG signals. To make the technology of high quality and low-cost fetal ECG widely accessible, we need algorithms for fetal ECG extraction from easily deployable aECG devices. We refer the readers to the up-to-date guidance (Behar et al., 2016) for more information about the non-invasive fetal ECG extraction and analysis.
The current study addresses this challenge by proposing an algorithm capable of working with only two composite (maternal and fetal) aECG channels to derive the fetal signal from it. It is based on the currently developed single-lead fECG algorithm based on the modern time-frequency analysis and manifold learning technique (Su and Wu, 2017) and a novel proposed diffusion-based channel selection criteria. All the proposed methods have rigorous mathematical backups, and numerically they can be efficiently implemented to handle long signal. We call the proposed algorithm SAVER, which stands for Smart AdaptiVe Ecg Recognition. To validate SAVER, we report the analysis results of two publicly available databases, and compare the algorithm with other available algorithms in the literature.
The paper is organized in the following way. In Section 2, we detail our proposed algorithm, describe the algorithms we will compare, and describe the databases we validate the algorithm. The results are shown in Section 3, and the discussions with the future works are provided in Section 4. The paper closes with the conclusion shown in Section 5. The necessary theoretical background is provided in SAVER Section SI.1 in appendix, particularly the diffusion-based channel selection criteria. We refer the readers to Su and Wu (2017) for the details of the de-shape short time Fourier transform (dsSTFT), beat tracking and the nonlocal median.
2.1. Two-Lead fECG Algorithm – SAVER
We now describe the proposed two-channel fECG algorithm, which the authors coined as SAVER. The overall algorithm is illustrated in Figure 1.
Figure 1. The flow chart of the proposed two-channel fECG algorithm, SAVER. The x-axis of all figures are of the unit second. The data is the a2 recording from the database used in the 2013 PhysioNet/Computing in Cardiology Challenge, and channel 1 and channel 4 are shown in this illustration. Only three linear combinations are shown for the illustration purpose. The signal quality index for the channel selection is shown on the third block.
Denote two simultaneously recorded aECG signals as with the sampling rate ξ0Hz over the interval from the 0-th second to the N/ξ0-th second. If the signal is sampled more slowly than 1,000 Hz, to enhance the R peak detection and the the nonlocal median (Su and Wu, 2017), the signal is upsampled to 1,000 Hz (Laguna and Sörnmo, 2000). We use the same notations to denote the upsampled signal.
Step 0: pre-processing
To suppress the noise, the signal is low-pass filtered below 100 Hz. Then, subtract the estimated trend from x0, y0, where the trends are estimated using median filter with window length LMF > 0 s. If needed, the power-line interference is suppressed by two notch-filters at 50 and 60 Hz, since the origin of the tested database in this paper is unknown (if the resource of the database is known, the notch-filter will be designed according to the power system of that region). Denote the pre-processed signal as x and y. Take a discrete finite subset . Define , where ; that is, zθ is a linear combination of two aECG signals. This linear combination could be viewed as a generalization of the augmentation technique considered in Andreotti et al. (2014, Section 2.3.3).
Step 1: maternal ECG estimation
We iterate the dsSTFT and nonlocal median algorithms proposed in Su and Wu (2017) to decompose the maECG from each linear combination in . The algorithm is summarized below. For each θ we run the following three sub-steps.
1. (step 1-1) Apply the dsSTFT to zθ and extract the dominant curve in the dsSTFT (Su and Wu, 2017, Section 3.1.2), which represents the estimated maternal IHR.
2. (step 1-2) Compute the polarity of zθ, where the polarity is either positive or negative. If the polarity of zθ is negative, multiply zθ by −1; that is, flip the sign of zθ. We use the same notation zθ to denote the polarity-corrected ECG signal. With the estimated maternal IHR and the polarity-corrected ECG signal, apply the beat tracking algorithm (Su and Wu, 2017, Section 3.1.3) to zθ to compute the locations of maternal R-peaks. Denote the timestamps of estimated maternal R peaks as , where kθ, m∈ ℕ is the number of estimated maternal R peaks.
3. (step 1-3) Adjust the estimated maternal R-peak locations by searching the maximum of zθ over a small window around . We use the same notation to denote the adjusted estimated maternal R-peak locations. Apply the nonlocal median (Su and Wu, 2017, Section 3.1.4) to estimate the maECG in zθ based on the estimated R-peak locations . Denote the estimated maECG as .
Step 2: channel selection
For each linear combination in , with the estimated maECG, we obtain a rough fECG by a simple subtraction:
Denote to be the set of rough fECG signals estimated from Step 1. We apply the lag map and the diffusion map (DM) to each rough fECG in and select the optimal linear combination by the following procedure. See Section SI.1 in the Appendix for the theoretical background of this approach.
For each rough fECG, say , we evaluate the signal quality index (SQI) for the channel selection purpose in the following way. Apply the L-step lag map to embed the interval [2, TCS + 2] seconds of into ℝL, where TCS > 0 is chosen by the user and 2 is chosen to avoid the boundary effect associated with the window in the dsSTFT approach. Here TCS is chosen to be short enough to guarantee the computational efficiency and to avoid the possibility non-stationarity inherited in the fECG signal, and long enough to capture the periodicity of the fECG. Denote the embedded point cloud as . Apply the 1-normalization DM to , where the bandwidth of the kernel is chosen in the following way suggested in Keller et al. (2010). We first set ϵ0 to be the smallest value such that each data point has at least one neighbor within the distance ϵ0. Then we set the bandwidth to be 2ϵ0. Denote ϕθ, 1 be the first nontrivial eigenvector of the corresponding graph Laplacian. Compute the power spectrum of ϕθ, 1, denoted as . Denote ξθ, 1, ξθ, 2, …, ξθ,nCS > 0, where nCS ∈ ℕ is the number of peaks chosen by the user, to be the frequencies associated with the highest nCS peaks in . Fix LCS > 0 and denote . The SQI for the channel selection purpose is thus defined as
Under the assumption that the better the quality of the rough fECG is, the closer the embedded point cloud is to the one-dimensional circle, we know that the higher the SQI, the better the rough fECG is. More precisely, if the embedded point cloud is close to the one-dimensional circle, the first non-trivial eigenvector should behave like an oscillatory function. With the designed SQI, we could choose the optimal rough fECG as the one with the highest SQI. Denote to be the optimal rough fECG with the highest signal quality index we can obtain from the given two channels.
Step 3: fetal R peaks estimation
With the rough fECG obtained from the optimal linear combination, we finish the algorithm by estimating the fetal R peaks and fECG by again applying the dsSTFT and the nonlocal median algorithm. This part of the algorithm is essentially the same as that for the maternal ECG estimation, and we repeat the three sub-steps below for the sake of completeness.
1. (step 3-1) Apply the dsSTFT to and extract the dominant curve in the dsSTFT, which represents the estimated fetal IHR.
2. (step 3-2) Compute the polarity of . If the polarity of is negative, multiply by −1, and use the same notation to denote the polarity-corrected ECG signal. With the estimated fetal IHR and the polarity-corrected ECG signal, apply the beat tracking algorithm to to compute the locations of maternal R-peaks. Denote the timestamps of estimated fetal R peaks as , where kf ∈ ℕ is the number of estimated fetal R peaks.
3. (step 3-3) Adjust the estimated fetal R-peak locations by searching the maximum of over a small window around rf, and use the same notation rf to denote the adjusted estimated fetal R-peak locations. Finally, output the fetal R peaks.
Remark 2.1. We mention that by applying the nonlocal median again based on rf, we could denoise the optimal rough fECG waveform and obtain a clean fetal waveform. However, since the result is similar to that shown in Su and Wu (2017), and the focus of this paper is the fetal R peak detection, we skip the details of the fECG reconstruction in this study, and leave the fetal waveform reconstruction in the future work.
2.2. Comparison with Benchmark Algorithms
There have been several algorithms proposed in the field suitable for analyzing fECG from multiple channel aECG signals. Note that the two-channel aECG signals fall in the category of the blind source separation (BSS) (De Lathauwer et al., 2000; Akhbari et al., 2013; Di Maria et al., 2014; Varanini et al., 2014) and its variations (Sameni et al., 2008; Haghpanahi and Borkholder, 2013; Akbari et al., 2015). It is well known that usually we need more than 4 channels to have a reasonable result (Andreotti et al., 2016). Due to the stationarity assumption of the ICA, the input signal should be truncated to be short enough, like 30 s long. An important step in the BSS approach is channel selection, which is critical to identify the decomposed channel that contains the maternal or fetal ECG. Although we only have two channels, for the comparison purpose, we still show the results of the BSS approaches, including the joint approximation diagonalization of eigen-matrices (JADE) for the independent component analysis (ICA) and the principal component analysis (PCA). Since there are only two decomposed signals, we do not carry out the channel selection algorithms proposed in, for example, Andreotti et al. (2014); instead, we take the ground truth annotation to select the optimal channel that is more likely to be the fECG, and report the detected R peaks from this detected channel. Note that we do not take the ground truth annotation into account in any other algorithms considered in this paper except this BSS approach, due to the limited number of channels. We apply the publicly available codes provided in http://www.fecgsyn.com, and call the PCA method BSSPCA and the ICA method BSSICA following the terminologies suggested in Andreotti et al. (2016).
Another set of algorithms allow us to take only single mECG signal, but need to simultaneously acquire the maternal thoracic-lead ECG signal (tECG). Examples include adaptive methods (AM) based on the least mean square (LMS) (Widrow et al., 1975) or the recursive least square (RLS) (Behar et al., 2014a) and its variations, like the echo state neural network (ESN) (Behar et al., 2014a), blind adaptive filtering (Graupe et al., 2008), extended Kalman filter (EKF) (Sameni, 2008; Niknazar et al., 2013; Andreotti et al., 2014), etc. In these algorithms, the maternal thoracic ECG signal (mtECG) is needed and is viewed as the reference channel. The mtECG contains the maternal cardiac activity information that we want to remove from the aECG. Based on the assumption that the mtECG and the maternal cardiac activity in the aECG are linearly related, the LMS or RLS helps to extract the fECG from the aECG by removing the maternal cardiac activity in the aECG. If the relationship between the tECG and the maternal cardiac activity in the aECG is nonlinear, then ESN could help. However, it is not always the case that we could get the mtECG, particularly in our setup, so these algorithms could not be directly applied for our purpose. Since it has been shown in Su and Wu (2017) that by combining the dsSTFT and nonlocal median, we are able to estimate the maECG signal accurately, we could thus view the estimated maECG signal as the reference. This consideration can also be found in, for example, Rodrigues (2014). For the LMS or ESN, we take the publicly available code from http://www.fecgsyn.com, and call the LMS method AMLMS and the ESN method AMESN following the terminologies suggested in Andreotti et al. (2016). We thus consider the following combinations of the proposed two channel fECG algorithm and the LMS or ESN. Precisely, in our proposed algorithm, we replace the direct subtraction (0.0.1) in Step 2 by the LMS or ESN, by taking the estimated maECG as the reference channel to get the rough fECG. We call the combined algorithm ds-AMLMS or ds-AMESN. Note that under the assumption that the nonlocal median does a good job to recover the maECG, the reference channel should be the same as, or linearly related to, the maternal cardiac activity in the aECG, so the LMS could be applied. The same idea could be applied to other algorithms, like RLS, but to keep the discussion simple, we focus on the above-mentioned two typical algorithms, LMS and ESN.
We could also consider the EKF algorithm. In the EKF algorithm, the information of the maternal R peak location is needed to cancel the maternal cardiac activity. Again, since it has been shown that by combining the dsSTFT and nonlocal median, we are able to estimate the maternal R peaks location accurately (Su and Wu, 2017), we could use the estimated maternal R peaks as the input to the EKF algorithm, and replace the direct subtraction (0.0.1) in Step 2 by the EKF. For the EKF, we take the publicly available code from http://www.fecgsyn.com, and call the EKF method TSEKF following the terminologies suggested in Andreotti et al. (2016). The combined algorithm is called the ds-TSEKF.
To have a complete comparison, we also consider the template subtraction (TS) algorithm, which is suitable for the single lead mECG signal, and replace the direct subtraction (0.0.1) in Step 2 by the TS algorithm. In this work, the TS method we apply is the singular value decomposition approach proposed in Kanjilal et al. (1997) and nominated in Behar et al. (2014b) and Andreotti et al. (2016) as TSPCA. For TSPCA, we take the publicly available code from http://www.fecgsyn.com. We call the combined algorithm ds-TSPCA. Other TS methods could be combined in the same way and we do not report the results to simplify the discussion.
We follow the suggested optimized parameters accompanying the code without any modification; for example, the input signal to the AM algorithms, like LMS or ESN, is resampled to 250 Hz1, the input signal to the TSEKF and TSPCA algorithms is resampled to 1000Hz, and we do not change the suggested initialization of the TSEKF code.
We validate the proposed two-channel algorithm on two publicly available databases of aECG signals.
The first database is the PhysioNet non-invasive fECG database (adfecgdb), where the aECG signals with the annotation provided by experts are publicly available https://www.physionet.org/physiobank/database/adfecgdb/ (Goldberger et al., 2000; Kotas et al., 2010). There are five pregnant women between 38 and 40 weeks of pregnancy in this database. Each has 4 aECG channels and one direct fECG signal recorded from the Komporel system (ITAM Institute, Zabrze, Poland)2. The four abdominal leads are placed around the navel, a reference lead is placed above the pubic symphysis, and a common mode reference electrode with active-ground signal is placed on the left leg. See Figure 2 for an illustration of the leads placement. The signal lasts for 5 min and is sampled at a fixed rate 1,000 Hz with the 16 bit resolution. The R peak annotation is determined from the direct fECG recorded from the fetal scalp lead.
The second database is the 2013 PhysioNet/Computing in Cardiology Challenge (https://physionet.org/challenge/2013/#data-sets), abbreviated as CinC2013. We focus on the set A composed of 75 recordings for an assessment of our proposed algorithm since it is the only one with the provided the R peak annotation with reference to a direct FECG signal, acquired from a fetal scalp electrode. Each recording includes four noninvasive mECG channels that were obtained from multiple sources using a variety of instrumentations with differing frequency response, resolution, and configurations. Although they are from different resources, all recordings are resampled at the sampling rate 1,000 Hz and last for 1 min. There is no publicly available information about where the leads are placed on the maternal abdomen. Note that some recordings come from the adfecgdb database, but no detail is available publicly. More details about these two databases can be found on the website. We follow the suggestion in Andreotti et al. (2014) to disregard the recording a54 since it was discarded by the Challenge's organizers, and focus on the remaining 74 recordings.
2.4. Evaluation Metrics
In the whole analysis, the R peak detection result is evaluated by beat-to-beat comparisons between the detected beats and the provided annotations. We follow the criterion in Guerrero-Martinez et al. (2006) and choose a matching window of 50 ms. Denote TP, FP, and FN to be true positive rate, false positive rate, and false negative rate, where TP means correctly detected peaks, FP means nonexistent peaks that were falsely detected, and FN means existing peaks that were not detected.
We report the sensitivity (SE) and the positive predictive value (PPV) defined as
and the F1 score, which is the harmonic mean of PPV and SE,
We also report the mean absolute error (MAE) of the estimated R peak locations. We follow the suggestion in Andreotti et al. (2016) to report the MAE only on true positive annotations to make the evaluation independent of the detection accuracy. Thus, the MAE is defined as
where nTP is the number of true positive annotations, and and are the temporal location of the i-th true positive reference R-peak and temporal location of the i-th true positive detected R peak.
For each database, we will report two sets of statistics. First, for each subject, we record the best F1 result among all pairs of available channels, denoted as F1(1) and report the mean and median of the F1(1) of all subjects, and the corresponding summary statistics of the MAE, denoted as MAE(1). To see how stable the algorithm is, we also record the median F1 result among all pairs of available channels, called F1(0.5), and report the mean and median of the F1(0.5) of all subjects, as well as the corresponding summary statistics of the MAE, denoted as MAE(0.5). Second, to evaluate the lead placement issue, for each pair of available channels, we report the the mean and median of the F1 of all subjects, and the corresponding summary statistics of the MAE. To avoid the boundary effect inevitable in the dsSTFT algorithm due to the window length, the first and last 2 s in every recording are not evaluated. The notation a±b indicates the mean a with the standard deviation b.
For a fair comparison and the reproducibility purposes, here we summarize the parameters for SAVER. The parameters are fixed for all signals throughout the paper unless otherwise stated. For the linear combination of two channels, we fix . The window length LMF of the median filter for the baseline wandering removal is chosen to be 0.1 s. For the dsSTFT, the beat tracking, and the nonlocal median, the parameters are set to be the same as those reported in Su and Wu (2017).
For the channel selection, we set the lag to L = 7 for the lag map; we choose the Gaussian kernel and α = 1 normalization for the DM; we choose TCS = 40, nCS = 6 and LCS = 0.1375 Hz for the adfecgdb database, and TCS = 10, nCS = 6 and LCS = 0.25 Hz for the CinC database. We mention that the above parameters are chosen in the ad-hoc fashion without any optimization pursue. Those parameters could be optimized based on the application field and the environment.
The algorithms are tested on MacBook Air (13-inch, Mid 2013) with Processor 1.3GHz Intel Core i5, Memory 4 GB1600MHz DDR3, Mac OS Sierra (Version 10.12.2), and Matlab R2015b without implementing the parallel computation.
For the adfecgdb database, the direct fECG measurement was lost between 187 and 191 s and between 203 and 211 s in the r10 record, and these two segments were discarded in the evaluation. The evaluation results of our proposed algorithm for each combination of two channels out of four available channels of all subjects in the adfecgdb database are shown in Table 1 for a clear comparison purpose. Except the combination of Channel 2 and Channel 3 in r01 and r08, all the other combinations have the F1 consistently greater than 94%. For the MAE, the result is always smaller than 9ms except the combination of Channel 2 and Channel 3 in r08. Table 2 shows the comparison of the proposed method with other available algorithms. The F1(1) and F1(0.5) of all 6 pairs for each subject are recorded, and the summary statistics of all subjects are shown. It is clear that SAVER is consistently better than the other algorithms. The average running time is 141.55 s for SAVER, 194.44 s for the ds-AMLMS, 589.83 s for the ds-AMESN, 308.93 s for ds-TSEKF, 78.30 s for ds-TSPCA, 10.01 s for BSSICA, and 10.98 s for BSSPCA.
Table 1. The results of SAVER: F1 score, mean absolute error (MAE), positive predictive value (PPV), and sensitivity (SE) of all pairs of two channels out of four available channels and all subjects over the whole 5 min signals in the adfecgdb database.
Table 2. The summary statistics of different methods' performance, including F1 and mean absolute error (MAE), evaluated in the adfecgdb database.
For the CinC2013 database, in Table 3 we compare SAVER with the other available algorithms in the CinC2013 database. The F1(1) of all recordings of our method is 92.99±16.0% and the corresponding MAE(1) is 5.38±4.52 ms, which are both better than the other compared methods. The median F1(0.5) of all recordings of our method is 85.44±22.42% and the MAE(0.5) of our method is 6.54±4.92 ms, which are both better than the best result determined by other methods. It should be noted that the median of F1(0.5) over 6 pairs of our proposed algorithm is still as high as 96.32%, while other methods decline dramatically to less than 60%3. This result suggests the stability of the proposed method4. The average running time is 20.29 s for SAVER, 27.26 s for the ds-AMLMS, 100.35 s for the ds-AMESN, 75.83 s for ds-TSEKF, 17.52 s for ds-TSPCA, 3.29 s for BSSICA, and 3.20 s for BSSPCA.
Table 3. The summary statistics of different methods' performance, including F1 and mean absolute error (MAE), evaluated in the CinC2013 database.
To further evaluate the influence of the lead placement, or to answer if we could design the best lead placement scheme for the proposed two-channel algorithm, we report the summary statistics of all pairs of two channels for the adfecgdb database in Table 4 and the CinC2013 database in Table 5. It is interesting to see that for the adfecgdb database, except for the combination of channel 2 and channel 3, the mean F1 accuracy is great than 97%. The outlier of the combination of channel 2 and channel 3 comes from the fact that the fECG is strong in case r08, which confuses the channel selection step. As a result, SAVER extracts the maternal ECG as the fECG, which leads to a wrong fECG estimation5. While determining the role of each component is a common issue for the fetal-maternal ECG separation algorithms and commonly we need more information to handle it, we leave this open problem for the future work.
Table 4. The summary statistics of SAVER, including F1 and mean absolute error (MAE), for six pairs of four available channels in the adfecgdb database.
Table 5. The summary statistics of SAVER, including F1 and mean absolute error (MAE), for six pairs of four available channels in the CinC2013 database.
Compared with the result of the adfecgdb database, the performance of SAVER in the CinC2013 database is not uniform cross different combinations of channels. Note that the lead placement scheme is unknown for the CinC2013 database, so it is not possible to conclude which pair of channels is the best. However, if we assume that the lead placement scheme for all recordings in the CinC2013 database is the same as the lead placement scheme shown in Figure 2, then the CinC2013 database results suggest that the best combination is channel 1 and channel 4; the F1 has the mean of 87.93% with the standard deviation 22.64%, and the median 97.60% with the interquartile range 6.92%; the MAE has the mean of 6.21 ms with the standard deviation 6.03 ms, and the median 4.34 ms with the interquartile range 5.62 ms6. Another finding deserves a discussion is that unlike the adfecgdb database, we can see the discrepancy between the best F1 out of the 6 pairs reported in Table 3 and the average F1 of each pair reported in Table 5. This might suggest that the lead system applied in the CinC2013 database is heterogenous across the recordings.
The encouraging results of SAVER indicate the possibility to design a “two-lead system” for the noninvasive, and long term fECG monitoring purpose. To the best of our knowledge, less is published about two aECG channels approach (for example, in Rodrigues, 2014, the considered algorithm can be applied to the two channel aECG), and our proposed method focuses on this direction. The main innovation of our approach, compared with other methods, is twofold. First, based on the geometry of the inherited oscillatory structure of the cardiac activity, the diffusion-based manifold learning technique is applied to do the channel section. While other channel selection criteria mainly are based on the power spectral distribution, wave morphology entropy, root mean square error, etc, to find the clearest and most enhanced QRS complexes (Di Maria et al., 2014; Ghaffari et al., 2015), our approach is different since we carefully examine the nontrivial underlying geometric structure hosting the cardiac activity by the DM and look for the linear combination that is most like a simple closed curve. Second, we apply the modern time-frequency analysis technique, the dsSTFT, and the beat tracking algorithms detailed in Su and Wu (2017) to obtain an accurate R peak locations, and the nonlocal median, to better estimate the maternal ECG morphology and fetal ECG morphology. Compared with other available algorithms, we use more information hidden in the aECG, including decomposing the non-sinusoidal oscillatory pattern from the time-varying frequency, and the low dimensional parametrization of all possible cardiac oscillations. We mention that an important advantage of the approach in Su and Wu (2017) is the ability to separate mECG and fECG with temporal overlap by the nonlocal median. Furthermore, due to its nonlocal nature, it can directly handle a long signal without dividing it into small fragments. Notice that unlike the traditional AF-like methods, SAVER does not cancel the maternal ECG in one channel by designing a filter from another channel; instead, it directly cancels the maternal ECG in a single linear combination, as is mentioned in Step 1.
Our results deserve a discussion and comparison with the previous reported findings. For the adfecgdb database, our result is overall compatible with, or better than, the state-of-art result reported in the field. For example, if we choose the pair of channel 1 and channel 2, our result is better than the best channel result based on the continuous wavelet transform based single-channel algorithm (Table 5, Castillo et al., 2013). However, it is not a fair comparison since the algorithm used in Castillo et al. (2013) is a single-channel algorithm. On the other hand, if we compare with the methods based on ICA on four channels (Table 1, Poian et al., 2015), our result is compatible. The MAE, which is less reported in the literature, is as small as 10 ms, which indicates the potential of applying the SAVER to do the fetal heart rate variability (HRV) analysis.
For the CinC2013 database, our result is compatible, or better than, the reported results. At the first glance, it is not the case, since by the ICA-based algorithms (Andreotti et al., 2014; Behar et al., 2014b), the accuracy could be as high as have the mean F1 = 96%, under the same setup that a detected R-peak was labeled as TP if within 50 ms of a reference R-peak. However, we mention that unlike SAVER, these algorithms are ICA-based and four channels are simultaneously used. Specifically, in Behar et al. (2014b, Table 3), among different combinations of different algorithms, the algorithm FUSE-SMOOTH achieved the best result – the mean F1 over all recordings is 96%, after removing a33, a38, a47, a52, a54, a71, and a74; in Andreotti et al. (2014, Table 1), the augmentation, the ICA, the template adaptation or TSEKF, and other techniques are applied, and the result with the mean F1 = 97.3% over all recordings with the standard deviation 0.108 is reported based on the template adaptation, after removing a54. Our proposed algorithm, on the other hand, outperforms the algorithm based on four channels and the BSSPCA, for example, Di Maria et al. (2014). In Di Maria et al. (2014, Section 3.2), the accuracy of the proposed algorithm in detecting the fetal heart beats gives the mean F1 = 89.8% over all recordings, under the setup that a detected R-peak was labeled as TP if within 100 ms of a reference R-peak and removing 9 recordings, including a29, a38, a54, a56, a33, a47, a52, a71, and a74. Another novel method based on the channel selection over 4 channels followed by the sequential total variation denoising (Lee and Lee, 2016, Table 5) leads to the accuracy with F1 = 89.9% and the MAE = 9.3 ms7 under the setup that a detected R-peak was labeled as TP if within 50 ms of a reference R-peak and removing a33, a38, a47, a52, a54, a71, and a74. We emphasize that while our algorithm does not outperform some of the above-mentioned algorithms, based on two channels, SAVER leads to the MAE as small as 6.21 ms in channel 1 and channel 4 combination in the CinC2013 database, which again indicates the potential of applying the SAVER to do the fetal HRV analysis.
As discussed above, theoretically, the chance is low that the fetal cardiac axis orientation would be so much orthogonal to the 2-dim affine subspace spanned by the two leads that no fECG shape can be reconstructed. This is a big advantage compared with the single-lead system, as the chance that the fetal cardiac axis orientation is orthogonal to the 1-dim affine subspace spanned by the single lead is much higher. Thus, while there have been several successful algorithms for the one aECG channel, like (Castillo et al., 2013; Behar et al., 2014a; Su and Wu, 2017) and the citations inside, if the recorded one channel signal does not have fECG information, there is nothing the algorithm can do. From the practical viewpoint, since only two leads are needed, the corresponding hardware could be lighter and more deployable than the currently available four-lead or multiple-lead systems. While it is certainly possible to generalize our algorithm to a three-lead or four-lead system (and the algorithm can be changed directly according to the setup), to have a better balance between the prediction accuracy, the hardware design, and practical purposes, we focus on the two-lead system in our research.
Despite of the above-mentioned benefits, there are several challenges we need to solve until this possible system is clinically usable. As is shown above, the performance of SAVER depends on how the two leads are put on the abdomen. The fECG situation is clearly different from the adult ECG system, like the widely applied 12 lead ECG system. Since fetus does move and rotate inside the uterus, the uterus differs from female to female, and the maternal body profile varies, we may not expect to have a two-lead system universal for all women. Therefore, for the practical purpose, particularly for the long term monitoring purpose and the future digital health, like the wearable biosensors (Li et al., 2017), it is important to ask if we could adaptively find the best lead placement scheme for different females. For the practical purpose, due to the inevitable non-stationary noise of different types, like the motion artifact and uterine contraction, an automatic system providing a SQI to alarm/warn the low quality of the lead system, and hence improve the overall fECG extraction quality, is urgently needed. We leave this important engineering problem to the future work. Another interesting question naturally raises from the current work is if we could generalize the current algorithm to study the twin dataset. Theoretically it is possible, if we take the fact that geometrically the twin will locate in different positions. We would expect to study this problem when the dataset is available.
From the algorithmic viewpoint, there are several directions we could improve the proposed two-channel fECG algorithm. The main ingredient in SAVER is the diffusion geometry. Since we have more than one aECG channel, we could consider modern diffusion-based manifold learning techniques to extract information common in two channels, like the alternating diffusion (Papyan and Talmon, 2016; Talmon and Wu, 2016; Lederman and Talmon, in press). The non-stationary nature of the fECG signal, which often presents itself as a time-varying frequency, might jeopardize the diffusion-based approach. We could consider to entangle the nontrivial time-varying frequency nature of the signal by further applying the modern nonlinear-type time-frequency analysis technique, like the synchrosqueezing transform or concentration of frequency and time (see Daubechies et al., 2016 and the citations inside). In this work, the parameters for the channel selection are chosen in the ad hoc fashion and are fixed across different algorithms for a fair comparison. For the practical purpose, we may optimize these parameters to improve the results. A systematic survey of this issue will be reported in the future work.
Another important algorithmic question left unanswered in this paper is how to improve the nonlocal median algorithm so that the reconstructed fECG could provide more accurate electrophysiological information about the heart, for example, the ECG morphology like the Q wave and ST-segment section information (Amer-Wåhlin et al., 2001). The main difficulty encountered in this problem is the lack of the “ground truth,” and a careful design of the clinical trial to acquire a reliable ground truth for the morphological study of the fetal cardiac activity is needed. As important as this clinical information could be, we will focus on it as an independent research and report the result in the future work.
Last but not the least, the databases we tested are small and not specifically designed for our purpose. We thus need a large scale and well designed prospective study to confirm the result.
A novel two-channel fetal-maternal ECG signal separation algorithm, SAVER, is proposed. The potential of the proposed algorithm is supported by the positive validation results on two publicly available databases. The algorithm is both computationally efficient and is supported by the underlying rigorous mathematical model and theory. Its clinical applicability will be evaluated in the future work.
RL and HW carried out the data analysis. MF and HW carried out the paper writeup.
Conflict of Interest Statement
The authors filed a patent for this method.
HW research is partially supported by Sloan Research Fellow FR-2015-65363. The authors acknowledge the help of Jan Hamanishi to prepare the illustration in Figure 2. The authors thank the reviewers for their constructive comments to improve the manuscript.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fphys.2017.00277/full#supplementary-material
1. ^We mention that the publicly available code for AM algorithms is optimized for the 250 Hz signal, and it is likely that the results would be improved if the corresponding parameters had been optimized for 1,000 Hz.
3. ^It is suggested in Behar et al. (2014b) that the TSPCA algorithm performs better if we applied the 10 Hz high pass filter to remove the baseline wandering before applying the TSPCA algorithm. If we replace the median filter in Step 0 by the 10 Hz high pass filter to remove the baseline wandering, over all subjects except a54, we have F1(1) = 93.19±14.19%, MAE(1) = 5.16±3.85 ms, F1(0.5) = 74.79±27.74%, and MAE(0.5) = 9.26±6.34 ms. Note that while the results of F1(1) and MAE(1) are both better than all algorithms under comparison, the results of F1(0.5) and MAE(0.5) are worse. This discrepancy might be caused by the carried out optimization in Behar et al. (2014b). Since the baseline wandering algorithm is out of the scope of this paper, in all tables we only report the results with the median filter as mentioned in Step 0.
4. ^It is suggested in Behar et al. (2014b, p. 1569) to remove six more recordings, a33, a38, a47, a52, a71, and a74, in addition to a54, because of some inaccurate reference annotations identified by the visual inspection of authors in Behar et al. (2014b). The F1(1) of all recordings of our method is 94.80±13.17% and the MAE(1) is 5.04±3.88 ms, and the F1(0.5) of all recordings of our method is 87.04±21.27% and the MAE(0.5) of our method is 6.29±4.56 ms.
5. ^If we are allowed to use the physiological information that both the fetus and the mother are healthy so that the fetal IHR is on average higher than maternal IHR, then we could correct this confusion by swapping the fetal IHR and maternal IHR. This leads to the mean F1 of the combination of channel 2 and channel 3 93.44% with the standard deviation 6.99% and the mean MAE 5.57 ms with the standard deviation 2.41 ms, and the results of other combinations unchanged.
6. ^If we remove a33, a38, a47, a52, a54, a71, and a74 from the CinC2013 database (Behar et al., 2014b), for the combination of channel 1 and channel 4, the F1 has the mean 89.81% with the standard deviation 20.84%, and the median becomes 98.41% with the interquartile range 5.10%; the MAE has the mean of 5.74 ms with the standard deviation 5.33 ms, and the median 4.20 ms with the interquartile range 5.48 ms.
7. ^In a private communication, the authors confirmed that this F1 is the “overall F1,” which is evaluated by collecting all beats from all recordings, and evaluate the F1 on all collected beats. If we follow the same procedure and remove a33, a38, a47, a52, a54, a71, and a74, the overall F1 of SAVER for the combination of channel 1 and channel 4 is 89.77% and the MAE is 4.91 ms.
Akhbari, M., Niknazar, M., Jutten, C., Shamsollahi, M. B., and Rivet, B. (2013). Fetal electrocardiogram R-peak detection using robust tensor decomposition and extended Kalman filtering. Comput. Cardiol. 40, 189–192. Available online at: http://ieeexplore.ieee.org/document/6712443/
Amer-Wåhlin, I., Hellsten, C., Norén, H., Hagberg, H., Herbst, A., Kjellmer, I., et al. (2001). Cardiotocography only versus cardiotocography plus st analysis of fetal electrocardiogram for intrapartum fetal monitoring: a swedish randomised controlled trial. Lancet 358, 534–538. doi: 10.1016/S0140-6736(01)05703-8
Anblagan, D., Pataky, R., Evans, M. J., Telford, E. J., Serag, A., Sparrow, S., et al. (2016). Association between preterm brain injury and exposure to chorioamnionitis during fetal life. Sci. Rep. 6:37932. doi: 10.1038/srep37932
Andreotti, F., Behar, J., Zaunseder, S., Oster, J., and Clifford, G. D. (2016). An open-source framework for stress-testing non-invasive foetal ECG extraction algorithms. Physiol. Meas. 37, 627–648. doi: 10.1088/0967-3334/37/5/627
Andreotti, F., Riedl, M., Himmelsbach, T., Wedekind, D., Wessel, N., Stepan, H., et al. (2014). Robust fetal ECG extraction and detection from abdominal leads. Physiol. Meas. 35, 1551–1567. doi: 10.1088/0967-3334/35/8/1551
Behar, J., Andreotti, F., Zaunseder, S., Oster, J., and Clifford, G. D. (2016). A practical guide to non-invasive foetal electrocardiogram extraction and analysis. Physiol. Meas. 37:R1. doi: 10.1088/0967-3334/37/5/R1
Behar, J., Oster, J., and Clifford, G. D. (2014b). Combining and benchmarking methods of foetal ECG extraction without maternal or scalp electrode data. Physiol. Meas. 35, 1569–1589. doi: 10.1088/0967-3334/35/8/1569
Belfort, M. A., Saade, G. R., Thom, E., Blackwell, S. C., Reddy, U. M., Thorp, J. M. Jr., et al. (2015). A randomized trial of intrapartum fetal ECG ST-segment analysis. N. Engl. J. Med. 373, 632–641. doi: 10.1056/NEJMoa1500600
Bravi, A., Green, G., Herry, C., Wright, H., Longtin, A., Kenny, G., et al. (2013). Do physiological and pathological stresses produce different changes in heart rate variability? Front. Physiol. 4:197. doi: 10.3389/fphys.2013.00197
Castillo, E., Morales, D. P., Botella, G., Garcia, A., Parrilla, L., and Palma, A. J. (2013). Efficient wavelet-based ECG processing for single-lead FHR extraction. Digital Signal Process. 23, 1897–1909. doi: 10.1016/j.dsp.2013.07.010
Daubechies, I., Wang, Y. G., and Wu, H.-t. (2016). ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform. Philos. Trans. A Math. Phys. Eng. Sci. 374:20150193. doi: 10.1098/rsta.2015.0193
Di Maria, C., Liu, C., Zheng, D., Murray, A., and Langley, P. (2014). Extracting fetal heart beats from maternal abdominal recordings: selection of the optimal principal components. Physiol. Meas. 35, 1649–1664. doi: 10.1088/0967-3334/35/8/1649
Durosier, L. D., Herry, C. L., Cortes, M., Cao, M., Burns, P., Desrochers, A., et al. (2015). Does heart rate variability reflect the systemic inflammatory response in a fetal sheep model of lipopolysaccharide-induced sepsis? Physiol. Meas. 36:2089. doi: 10.1088/0967-3334/36/10/2089
Durosier, L. D., Green, G., Batkin, I., Seely, A. J., Ross, M. G., Richardson, B. S., et al. (2014). Sampling rate of heart rate variability impacts the ability to detect acidemia in ovine fetuses near-term. Front. Pediatr. 2:38. doi: 10.3389/fped.2014.00038
Fairchild, K. D., Sinkin, R. A., Davalian, F., Blackman, A. E., Swanson, J. R., Matsumoto, J. A., et al. (2014). Abnormal heart rate characteristics are associated with abnormal neuroimaging and outcomes in extremely low birth weight infants. J. Perinatol. 34, 375–379. doi: 10.1038/jp.2014.18
Fairchild, K. D., Srinivasan, V., Randall Moorman, J., Gaykema, R. P. A., and Goehler, L. E. (2011). Pathogen-induced heart rate changes associated with cholinergic nervous system activation. Am. J. Physiol. 300, R330–R339. doi: 10.1152/ajpregu.00487.2010
Garzoni, L., Faure, C., and Frasch, M. (2013). Fetal cholinergic anti-inflammatory pathway and necrotizing enterocolitis: the brain-gut connection begins in utero. Front. Integr. Neurosci. 7:57. doi: 10.3389/fnint.2013.00057
Ghaffari, A., Mollakazemi, M. J., Atyabi, S. A., and Niknazar, M. (2015). Robust fetal QRS detection from noninvasive abdominal electrocardiogram based on channel selection and simultaneous multichannel processing. Australas Phys. Eng. Sci. Med. 38, 581–592. doi: 10.1007/s13246-015-0381-2
Goldberger, A., Amaral, L., Glass, L., Hausdorff, J., Ivanov, P., Mark, R., et al. (2000). Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals. Circulation 101, e215–e220. doi: 10.1161/01.CIR.101.23.e215
Graupe, D., Graupe, M. H., Zhong, Y., and Jackson, R. K. (2008). Blind adaptive filtering for non-invasive extraction of the fetal electrocardiogram and its non-stationarities. Proc. Inst. Mech. Eng. H J. Eng. Med. 222, 1221–1234. doi: 10.1243/09544119JEIM417
Guerrero-Martinez, J. F., Martinez-Sober, M., Bataller-Mompean, M., and Magdalena-Benedito, J. R. (2006). “New algorithm for fetal QRS detection in surface abdominal records,” in Computers in Cardiology (IEEE), 441–444.
Hagberg, H., Mallard, C., Ferriero, D. M., Vannucci, S. J., Levison, S. W., Vexler, Z. S., et al. (2015). The role of inflammation in perinatal brain injury. Nat. Rev. Neurol. 11, 192–208. doi: 10.1038/nrneurol.2015.13
Keener, J. (1998). Mathematical Physiology. Springer. Available online at: http://www.springer.com/cn/book/9780387227061
Li, X., Dunn, J., Salins, D., Zhou, G., Zhou, W., Schüssler-Fiorenza Rose, S. M., et al. (2017). Digital health: tracking physiomes and activity using wearable biosensors reveals useful health-related information. PLoS Biol. 15:e2001402. doi: 10.1371/journal.pbio.2001402
Li, X., Xu, Y., Herry, C., Durosier, L. D., Casati, D., Stampalija, T., et al. (2015). Sampling frequency of fetal heart rate impacts the ability to predict pH and BE at birth: a retrospective multi-cohort study. Physiol. Meas. 36:L1. doi: 10.1088/0967-3334/36/5/L1
Liu, H. L., Garzoni, L., Herry, C., Durosier, L. D., Cao, M., Burns, P., et al. (2016). Can monitoring fetal intestinal inflammation using heart rate variability analysis signal incipient necrotizing enterocolitis of the neonate? Pediatr. Crit. Care Med. 17, 165–176. doi: 10.1097/PCC.0000000000000643
Niknazar, M., Rivet, B., and Jutten, C. (2013). Fetal ECG extraction by extended state Kalman filtering based on single-channel recordings. IEEE Trans. Biomed. Eng. 60, 1345–1352. doi: 10.1109/TBME.2012.2234456
Olofsson, P. S., Rosas-Ballina, M., Levine, Y. A., and Tracey, K. J. (2012). Rethinking inflammation: neural circuits in the regulation of immunity. Immunol. Rev. 248, 188–204. doi: 10.1111/j.1600-065X.2012.01138.x
Poian, G. D., Bernardini, R., and Rinaldo, R. (2015). “Sparse representation for fetal QRS detection in abdominal ECG recordings,” in The 5th IEEE International Conference on E-Health and Bioengineering - EHB 2015 (Iasi), 1–4. doi: 10.1109/ehb.2015.7391501
Sameni, R., Jutten, C., and Shamsollahi, M. (2008). Multichannel electrocardiogram decomposition using periodic component analysis. IEEE Trans. Biomed. Eng. 55, 1935–1940. doi: 10.1109/TBME.2008.919714
Stark, J., Broomhead, D., Davies, M., and Huke, J. (1997). Takens embedding theorems for forced and stochastic systems. Nonlin. Anal. Theory Methods Appl. 30, 5303–5314. doi: 10.1016/S0362-546X(96)00149-6
Takens, F. (1981). “Detecting strange attractors in turbulence,” in Dynamical Systems and Turbulence, Vol. 898 of Lecture Notes in Mathematics, eds D. Rand and L.-S. Young (Berlin; Heidelberg: Springer), 366–381.
Varanini, M., Tartarisco, G., Billeci, L., Macerata, A., Pioggia, G., and Balocchi, R. (2014). An efficient unsupervised fetal QRS complex detection from abdominal maternal ECG. Physiol. Meas. 35, 1607–1619. doi: 10.1088/0967-3334/35/8/1607
Keywords: de-shape short time Fourier transform, fetal electrocardiogram, maternal abdominal electrocardiogram, nonlocal median, diffusion maps
Citation: Li R, Frasch MG and Wu H-T (2017) Efficient Fetal-Maternal ECG Signal Separation from Two Channel Maternal Abdominal ECG via Diffusion-Based Channel Selection. Front. Physiol. 8:277. doi: 10.3389/fphys.2017.00277
Received: 04 February 2017; Accepted: 18 April 2017;
Published: 16 May 2017.
Edited by:Ahsan H. Khandoker, Khalifa University, United Arab Emirates
Reviewed by:Joachim Behar, Technion - Israel Institute of Technology, Israel
Paolo Melillo, Università degli Studi della Campania “Luigi Vanvitelli” Caserta, Italy
Copyright © 2017 Li, Frasch and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.