# Efficient Fetal-Maternal ECG Signal Separation from Two Channel Maternal Abdominal ECG via Diffusion-Based Channel Selection

^{1}Department of Mathematics, University of Toronto, Toronto, ON, Canada^{2}Department of Obstetrics and Gynecology, University of Washington, Seattle, USA^{3}Mathematics 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 *F*_{1} 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.

## 1. Introduction

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 **S**mart **A**dapti**V**e **E**cg **R**ecognition. 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. Methods

### 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 ${\text{x}}_{0},{\text{y}}_{0}\in {\mathbb{R}}^{N}$ with the sampling rate ξ_{0}Hz 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 **x**_{0}, **y**_{0}, where the trends are estimated using median filter with window length *L*_{MF} > 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 ${I}\subset (-1,1]$. Define ${\text{z}}_{\theta}=\theta \text{x}+\sqrt{1-{\theta}^{2}}\text{y}$, where $\theta \in {I}$; 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 ${\left\{{\text{z}}_{\theta}\right\}}_{\theta \in {I}}$. 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 ${\text{r}}_{\theta}^{m}=({r}_{\theta ,1}^{m},\dots ,{r}_{\theta ,{k}_{\theta ,m}}^{m})$, 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 ${\text{r}}_{\theta}^{m}$. We use the same notation ${\text{r}}_{\theta}^{m}$ 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 ${\text{r}}_{\theta}^{m}$. Denote the estimated maECG as ${\stackrel{~}{\text{z}}}_{\theta ,m}$.

**Step 2: channel selection**

For each linear combination in ${\left\{{\text{z}}_{\theta}\right\}}_{\theta \in {I}}$, with the estimated maECG, we obtain a rough fECG by a simple subtraction:

Denote ${\left\{{\stackrel{~}{\text{z}}}_{\theta ,f}\right\}}_{\theta \in {I}}$ 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 ${\left\{{\stackrel{~}{\text{z}}}_{\theta ,f}\right\}}_{\theta \in {I}}$ 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 ${\stackrel{~}{\text{z}}}_{\theta ,f}$, 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, *T*_{CS} + 2] seconds of ${\stackrel{~}{\text{z}}}_{\theta ,f}$ into ℝ^{L}, where *T*_{CS} > 0 is chosen by the user and 2 is chosen to avoid the boundary effect associated with the window in the dsSTFT approach. Here *T*_{CS} 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 ${{X}}_{\theta ,f}\subset {\mathbb{R}}^{L}$. Apply the 1-normalization DM to ${{X}}_{\theta ,f}$, 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 $|{\widehat{\varphi}}_{\theta ,1}{|}^{2}$. Denote ξ_{θ, 1}, ξ_{θ, 2}, …, ξ_{θ,nCS} > 0, where *n*_{CS} ∈ ℕ is the number of peaks chosen by the user, to be the frequencies associated with the highest *n*_{CS} peaks in $|{\widehat{\varphi}}_{\theta ,1}{|}^{2}$. Fix *L*_{CS} > 0 and denote ${{J}}_{\theta}:={\cup}_{j=1}^{{n}_{\mathtt{\text{CS}}}}\left[{\xi}_{\theta ,i}-{L}_{\mathtt{\text{CS}}},{\xi}_{\theta ,i}+{L}_{\mathtt{\text{CS}}}\right]$. 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 ${\stackrel{~}{\text{z}}}_{f}^{*}$ 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 ${\stackrel{~}{\text{z}}}_{f}^{*}$ 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 ${\stackrel{~}{\text{z}}}_{f}^{*}$ and extract the dominant curve in the dsSTFT, which represents the estimated fetal IHR.

2. (step 3-2) Compute the polarity of ${\stackrel{~}{\text{z}}}_{f}^{*}$. If the polarity of ${\stackrel{~}{\text{z}}}_{f}^{*}$ is negative, multiply ${\stackrel{~}{\text{z}}}_{f}^{*}$ by −1, and use the same notation ${\text{z}}_{f}^{*}$ to denote the polarity-corrected ECG signal. With the estimated fetal IHR and the polarity-corrected ECG signal, apply the beat tracking algorithm to ${\stackrel{~}{\text{z}}}_{f}^{*}$ to compute the locations of maternal R-peaks. Denote the timestamps of estimated fetal R peaks as ${\text{r}}^{f}=({r}_{1}^{f},\dots ,{r}_{{k}_{f}}^{f})$, where *k*_{f} ∈ ℕ is the number of estimated fetal R peaks.

3. (step 3-3) Adjust the estimated fetal R-peak locations by searching the maximum of ${\stackrel{~}{\text{z}}}_{f}^{*}$ over a small window around **r**^{f}, and use the same notation **r**^{f} 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 **r**^{f}, we could denoise the optimal rough fECG waveform ${\stackrel{~}{\text{z}}}_{f}^{*}$ 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 BSS_{PCA} and the ICA method BSS_{ICA} 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 AM_{LMS} and the ESN method AM_{ESN} 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-AM_{LMS} or ds-AM_{ESN}. 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 TS_{EKF} following the terminologies suggested in Andreotti et al. (2016). The combined algorithm is called the ds-TS_{EKF}.

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 TS_{PCA}. For TS_{PCA}, we take the publicly available code from http://www.fecgsyn.com. We call the combined algorithm ds-TS_{PCA}. 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 Hz^{1}, the input signal to the TS_{EKF} and TS_{PCA} algorithms is resampled to 1000Hz, and we do not change the suggested initialization of the TS_{EKF} code.

### 2.3. Materials

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 *F*_{1} 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 *n*_{TP} is the number of true positive annotations, and ${\stackrel{~}{r}}_{i}^{f}$ and ${r}_{i}^{f}$ 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 *F*_{1} result among all pairs of available channels, denoted as *F*_{1}(1) and report the mean and median of the *F*_{1}(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 *F*_{1} result among all pairs of available channels, called *F*_{1}(0.5), and report the mean and median of the *F*_{1}(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 *F*_{1} 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*.

### 2.5. Parameters

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 ${I}={\left\{-1+k/6\right\}}_{k\text{}=\text{}1}^{12}$. The window length *L*_{MF} 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 *T*_{CS} = 40, *n*_{CS} = 6 and *L*_{CS} = 0.1375 Hz for the adfecgdb database, and *T*_{CS} = 10, *n*_{CS} = 6 and *L*_{CS} = 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.

## 3. Results

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 *F*_{1} 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 *F*_{1}(1) and *F*_{1}(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-AM_{LMS}, 589.83 s for the ds-AM_{ESN}, 308.93 s for ds-TS_{EKF}, 78.30 s for ds-TS_{PCA}, 10.01 s for BSS_{ICA}, and 10.98 s for BSS_{PCA}.

**Table 1. The results of SAVER: F_{1} 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 F_{1} 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 *F*_{1}(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 *F*_{1}(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 *F*_{1}(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 method^{4}. The average running time is 20.29 s for SAVER, 27.26 s for the ds-AM_{LMS}, 100.35 s for the ds-AM_{ESN}, 75.83 s for ds-TS_{EKF}, 17.52 s for ds-TS_{PCA}, 3.29 s for BSS_{ICA}, and 3.20 s for BSS_{PCA}.

**Table 3. The summary statistics of different methods' performance, including F_{1} 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 *F*_{1} 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 estimation^{5}. 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 F_{1} and mean absolute error (MAE), for six pairs of four available channels in the adfecgdb database**.

**Table 5. The summary statistics of SAVER, including F_{1} 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 *F*_{1} 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 ms^{6}. Another finding deserves a discussion is that unlike the adfecgdb database, we can see the discrepancy between the best *F*_{1} out of the 6 pairs reported in Table 3 and the average *F*_{1} of each pair reported in Table 5. This might suggest that the lead system applied in the CinC2013 database is heterogenous across the recordings.

## 4. Discussion

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 *F*_{1} = 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 *F*_{1} 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 TS_{EKF}, and other techniques are applied, and the result with the mean *F*_{1} = 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 BSS_{PCA}, 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 *F*_{1} = 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 *F*_{1} = 89.9% and the MAE = 9.3 ms^{7} 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.

## 5. Conclusion

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.

## Author Contributions

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.

## Acknowledgments

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.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fphys.2017.00277/full#supplementary-material

## Footnotes

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.

2. ^http://www.itam.zabrze.pl/developments-english-version-233/665-komporel.

3. ^It is suggested in Behar et al. (2014b) that the TS_{PCA} algorithm performs better if we applied the 10 Hz high pass filter to remove the baseline wandering before applying the TS_{PCA} 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 *F*_{1}(1) = 93.19±14.19%, MAE(1) = 5.16±3.85 ms, *F*_{1}(0.5) = 74.79±27.74%, and MAE(0.5) = 9.26±6.34 ms. Note that while the results of *F*_{1}(1) and MAE(1) are both better than all algorithms under comparison, the results of *F*_{1}(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 *F*_{1}(1) of all recordings of our method is 94.80±13.17% and the MAE(1) is 5.04±3.88 ms, and the *F*_{1}(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 *F*_{1} 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 *F*_{1} 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 *F*_{1} is the “overall *F*_{1},” which is evaluated by collecting all beats from all recordings, and evaluate the *F*_{1} on all collected beats. If we follow the same procedure and remove a33, a38, a47, a52, a54, a71, and a74, the overall *F*_{1} of SAVER for the combination of channel 1 and channel 4 is 89.77% and the MAE is 4.91 ms.

## References

Akbari, H., Shamsollahi, M. B., and Phlypo, R. (2015). “Fetal ECG extraction using πTucker decomposition,” in *International Conference on Systems, Signals and Image Processing* (London: IEEE).

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., Johnson, A., Clifford, G. D., and Oster, J. (2014a). A comparison of single channel fetal ECG extraction methods. *Ann. Biomed. Eng.* 42, 1340–1353. doi: 10.1007/s10439-014-0993-9

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

Coifman, R. R., and Lafon, S. (2006). Diffusion maps. *Appl. Comput. Harmon. Anal.* 21, 5–30. doi: 10.1016/j.acha.2006.04.006

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

De Lathauwer, L., De Moor, B., and Vandewalle, J. (2000). Fetal electrocardiogram extraction by blind source subspace separation. *IEEE Trans. Biomed. Eng.* 47, 567–572. doi: 10.1109/10.841326

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

El Karoui, N. (2010). On information plus noise kernel random matrices. *Ann. Stat.* 38, 3191–3216. doi: 10.1214/10-AOS801

El Karoui, N., and Wu, H.-T. (2016). Connection graph Laplacian methods can be made robust to noise. *Ann. Stat.* 44, 346–372. doi: 10.1214/14-AOS1275

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

Haghpanahi, M., and Borkholder, D. A. (2013). “Fetal ECG extraction from abdominal recordings using array signal processing,” in *Computing in Cardiology*, Vol. 40 (IEEE), 173–176.

Jenkins, H. M. (1989). Thirty years of electronic intrapartum fetal heart rate monitoring: discussion paper. *J. R. Soc. Med.* 82, 210–214. doi: 10.1177/014107688908200410

Kanjilal, P. P., Palit, S., and Saha, G. (1997). Fetal ECG extraction from single-channel maternal ECG using singular value decomposition. *IEEE Trans. Biomed. Eng.* 44, 51–59. doi: 10.1109/10.553712

Keener, J. (1998). *Mathematical Physiology*. Springer. Available online at: http://www.springer.com/cn/book/9780387227061

Keller, Y., Coifman, R. R., Lafon, S., and Zucker, S. W. (2010). Audio-visual group recognition using diffusion maps. *IEEE Trans. Signal Process.* 58, 403–413. doi: 10.1109/TSP.2009.2030861

Kotas, M., Jezewski, J., Matonia, A., and Kupka, T. (2010). Towards noise immune detection of fetal QRS complexes. *Comput. Methods Prog. Biomed.* 97, 241–256. doi: 10.1016/j.cmpb.2009.09.005

Laguna, P., and Sörnmo, L. (2000). Sampling rate and the estimation of ensemble variability for repetitive signals. *Med. Biol. Eng. Comput.* 38, 540–546. doi: 10.1007/BF02345750

Lederman, R. R., and Talmon, R. (in press). Learning the geometry of common latent variables using alternating-diffusion. *Appl. Comp. Harmon. Anal.* doi: 10.1016/j.acha.2015.09.002

Lee, K., and Lee, B. (2016). Sequential total variation denoising for the extraction of fetal ECG from single-channel maternal abdominal ECG. *Sensors* 16:1020. doi: 10.3390/s16071020

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

Lin, C.-Y., Li, S., and Wu, H.-t. (2016). Wave-shape function analysis–when cepstrum meets time-frequency analysis. arXiv:1605.01805 *J. Four. Anal. Appl.* doi: 10.1007/s00041-017-9523-0

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

Richter, M., Schreiber, T., and Kaplan, D. (1998). Fetal ECG extraction with nonlinear state-space projections. *IEEE Trans. Biomed. Eng.* 45, 133–137. doi: 10.1109/10.650369

Rodrigues, R. (2014). Fetal beat detection in abdominal ECG recordings: global and time adaptive approaches. *Physiol. Meas.* 35, 1699–1711. doi: 10.1088/0967-3334/35/8/1699

Sameni, R. (2008). *Extraction of Fetal Cardiac Signals From An Array of Maternal Abdominal Recordings*. Ph.D. thesis, Sharif University of Technology – Institut National Polytechnique de Grenoble.

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

Singer, A., and Wu, H.-T. (2017). Spectral convergence of the connection laplacian from random samples. *Inform. Inference J. IMA*. 6, 58–123 doi: 10.1093/imaiai/iaw016

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

Stark, J., Broomhead, D., Davies, M., and Huke, J. (2003). Delay embeddings for forced systems. II. Stochastic forcing. *J. Nonlin. Sci.* 13, 519–577. doi: 10.1007/s00332-003-0534-4

Su, L., and Wu, H.-T. (2017). Extract fetal ECG from single-lead abdominal ECG by de-shape short time Fourier transform and nonlocal median. *Front. Appl. Math. Stat.* 3:2. doi: 10.3389/fams.2017.00002

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.

Talmon, R., and Wu, H.-T. (2016). Latent common manifold learning with alternating diffusion: analysis and applications. *arXiv preprint arXiv:1602.00078*.

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 EmiratesReviewed by:

Joachim Behar, Technion - Israel Institute of Technology, IsraelPaolo 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.

*Correspondence: Martin G. Frasch, mfrasch@uw.edu

Hau-Tieng Wu, hauwu@math.toronto.edu