Fetal Electrocardiogram Extraction from the Mother’s Abdominal Signal Using the Ensemble Kalman Filter

Fetal electrocardiogram (fECG) assessment is essential throughout pregnancy to monitor the wellbeing and development of the fetus, and to possibly diagnose potential congenital heart defects. Due to the high noise incorporated in the abdominal ECG (aECG) signals, the extraction of fECG has been challenging. And it is even a lot more difficult for fECG extraction if only one channel of aECG is provided, i.e., in a compact patch device. In this paper, we propose a novel algorithm based on the Ensemble Kalman filter (EnKF) for non-invasive fECG extraction from a single-channel aECG signal. To assess the performance of the proposed algorithm, we used our own clinical data, obtained from a pilot study with 10 subjects each of 20 min recording, and data from the PhysioNet 2013 Challenge bank with labeled QRS complex annotations. The proposed methodology shows the average positive predictive value (PPV) of 97.59%, sensitivity (SE) of 96.91%, and F1-score of 97.25% from the PhysioNet 2013 Challenge bank. Our results also indicate that the proposed algorithm is reliable and effective, and it outperforms the recently proposed extended Kalman filter (EKF) based algorithm.


Introduction
A national study reported by the Centers for Disease Control (CDC) showed that the U.S. fetal mortality rate remained unchanged from 2006 through 2013 from the rate of 6.05 to 5.96 per 1000 births [1]. Thus, fetal monitoring is essential throughout pregnancy for recognition of elements that might imperil the life of the fetus and mother. However, the 2015 Cochrane review of gold standard antenatal cardiotocography (CTG) for fetal assessment showed no clear evidence that it improves perinatal outcome [2].
The current global COVID-19 pandemic has led to critical demand for person-centered healthcare instead of hospital-centered healthcare system. Expectant moms usually come to the obstetric clinics many times during pregnancy for checkups. According to the CDC, a study on 461,825 women with COVID-19 showed that pregnant women with COVID-19 are more likely admitted to an intensive care unit, receiving invasive ventilation, extracorporeal membrane oxygenation, or die compared with non-pregnant women [3]. This would stop pregnant women from having physiological measurements, ultrasound examination, or non-stress test, which are critical for both the mom and the unborn baby. This calls for an urgent need for novel home-based tools and systems for reliable prenatal monitoring.
There are already existing ultrasound-based consumer products for fHR assessment in the home setting. Those require active scanning over the abdomen coupling with ultrasound gel to locate the fetal heart to obtain the fHR, which is highly technique dependent, easily make difficult by maternal of fetal movement and larger body size. Thus, it is especially challenging for non-medical persons to administer it. Further, the U.S. Food and Drug Administration (FDA) issued a warning in 2014 against the use such ultrasound fHR home monitors [4], citing safety concerns over their repeated use.
The fetal electrocardiogram (fECG) can also be used to assess fetal wellbeing and even enable further diagnoses. Precisely-extracted fECG will provide vital information about fHR, fetal development, fetal maturity, and the existence of abnormalities or distress during pregnancy [5]. However, continuous and noninvasive fECG monitoring has remained challenging in the research community in terms of both signal acquisition and processing. There are several methods proposed for extracting fECG from aECG. These methods can be generally classified into three groups: blind source separation (BSS), template subtraction, and filtering techniques. The BSS methods include methods such as parallel linear predictor (PLP) filter, principal component analysis (PCA), independent component analysis (ICA), and periodic component analysis (πCA) [6][7][8][9]. The BSS methods consider that the abdominal signal is a combination of fECG, mECG, and noises [10]. Although these methods perform greatly for fECG extraction, they need multiple-channel aECG signals, which makes them unsuitable for continuous non-invasive fECG monitoring. In addition, after extraction, the order of the separated independent component could not be determined. Thus, it is challenging to identify the fECG component for further processing [11]. Therefore, the BSS methods usually require the determination of other parameters (e.g., t-test, correlation coefficient, heart rate) to automatically identify the extracted components [12][13][14]. Template subtraction (TS) is another widely used approach. The method involves subtracting a synthetic mECG generated by estimating the QRS complex waveform (mQRS) of mECG from the abdomen signal [15][16][17][18][19][20]. The main challenge of this method involves mQRS detection [21], which becomes more challenging if the fetal R waves overlap with the maternal R waves. This drawback degrades the effectiveness of the template subtraction method for fECG extraction. The popular filtering techniques include adaptive filtering [22][23][24][25], Kalman filtering [26][27][28], and wavelet transform [29,30]. These filtering techniques are mostly and effectively applied for denoising of single-channel ECG signals. Adaptive filtering-based algorithms have been proposed for fECG extraction [31]. Such methods, however, require additional reference signals for separating the different components of the aECG.
Our group has been developing abdominal patch devices and systems for the acquisition and extraction of fECG [12,31,32]. As we aim to make the device compact and unobtrusive, a single recording channel is desired for saving space; consequently, fECG extraction is almost impossible, especially with the presence of motion noise [33]. In this case, there are multiple signals buried in the single-channel aECG including maternal ECG (mECG), fECG, maternal muscle activity, fetal movement activity, and noise. Here, we propose a novel algorithm based on the Ensemble Kalman filter (EnKF) to extract fECG from a single-channel aECG signal. The EnKF is an approximate filtering method that is used for state estimation of large-scale nonlinear systems. Specifically, the EnKF is a Sequential Monte Carlo (SMC) method, and it shows a better performance, especially in systems with strong nonlinearity, than the popular Extended Kalman filter (EKF) that applies analytic linearization. We studied and compared the performance of EKF and EnKF using our aECG data obtained from 10 pregnant subjects, and data from the PhysioNet database. The results demonstrated that our EnKF algorithm is more effective and outperforms the regular EKF in accurately extracting fetal ECG.

Ensemble Kalman Filter (EnKF)
EnKF is a variant of the celebrated Kalman filter used to estimate time-varying parameters in problems that arise in various disciplines [34]. It is applicable for problems that can be represented as dynamic systems and formulated in a state-space model with unknown time-varying state parameters. When the state-space model of the system is linear and Gaussian, the optimal estimate of the state parameters can be obtained using the Kalman filter [35]. However, when the problem is nonlinear or non-Gaussian, other variants of the Kalman filter, such as the EnKF, are used to obtain close-to-optimal solutions [36]. Our EnKF algorithm is developed by considering a Bayesian filtering framework and formulating the fECG extraction problem as a dynamic system whose state and measurement equations are represented in a state-space form. The dynamic model is adopted from the models proposed by McSharry et al. [37] and later discretized by Sameni [38].
Suppose the unknown time-varying state vector of a dynamic state-space model is denoted by x n ∈ R D x where n = 1, 2, . . . , N represents time instants and D x represents the dimension of x n . We assume that x n has a Markovian property, and its evolution is given by: where f (.) represents a state function which, in general, is nonlinear, and u n denotes the state noise vector with a known probability density function (pdf). Furthermore, the observation equation of the state-space model is given by: where y n ∈ R D y denotes the measurement vector obtained at time n, D y represents the dimension of the vector y n , and w n denotes the measurement noise vector whose pdf is assumed known. Given the state-space model (1) and (2), our objective is to make a sequential estimate of the evolution of the state vector x 1:n = {x 1 , . . . , x n } in real-time as the measurement vector denoted by y 1:n = {y 1 , . . . , y n } becomes available.
The Ensemble Kalman filter is a variant of the Kalman filter where the state error statistics are approximated using the Monte Carlo method. Recall that if the state Equation (1) and the measurement Equation (2) are linear, and the state noise u n and measurement noise w n are Gaussian, the optimal estimate of the state vector can be analytically obtained using Kalman filter. To understand the EnKF, let us review the Kalman filter algorithm's two steps: the time update and the measurement update. For convenience, we rewrite the state and measurement equations for a linear and Gaussian systems as follows: y n = Hx n + w n x n = Fx n−1 + u n , where F is a D x × D x matrix, H is a D y × D x matrix, and u n and w n are zero-mean Gaussian probability densities with covariances Q u and Q w , respectively. The time update step of the Kalman filter algorithm obtains the predicted value of the state vector and the covariance of its error. The predicted state vector and the covariance of its error at time n are computed from propagating their corresponding values at time n − 1 through the state dynamics as follows: wherex n|n−1 is the predicted value of the state vector, P n|n−1 is the covariance of the error of the predicted state vector given by is the expectation operator. Once the measurement vector is received, the measurement update step is used to compute the estimate of the state vector,x n and the covariance of the error of the state vector, P n , by applying corrections to the corresponding predicted values based on the measurement obtained as follows: where K n is the Kalman gain, and I is a diagonal matrix. We note that the Kalman gain can be expressed as K n = P xy,n P yy,n −1 , where P xy,n is the cross-covariance of the error of x n|n−1 and y n , and P yy,n is the covariance of the error of y n . These covariances are given by: When the state equation and measurement equation are not linear, the above equations of the covariances of the errors, P n|n−1 and P n , as well as the equation of the Kalman gain, K n , are not valid, and, therefore, generally, the Kalman gain cannot be analytically determined. A common approach of circumventing such problems in nonlinear state-space models is to apply extended Kalman filter (EKF), which is a variant of the Kalman filter, that approximates the state and measurement equations by linearizing them using the Taylor series. Consequently, the EKF approximates the posterior probability density of the state vector by a Gaussian distribution. When the true posterior distribution of the state vector is not 'close' to Guassian, such approximations may not be valid and the EKF may diverge. In such cases, sequential Monte Carlo methods show superior performance over EKF [39].
EnKF is the Monte Carlo-based Kalman filter which can be used for nonlinear and non-Gaussian models. The underlying idea of the method is to approximate the Kalman gain and state vector propagations using Monte Carlo technique. EnKF computes the Kalman gain by approximating P xy,n and P yy,n using their corresponding sample covariances, P xy,n andP yy,n . To do so, N number of ensembles, , are first drawn from the prior probability density of the state vector, p(x n|n−1 ), which has the same probability distribution function as the state noise with a mean of f (x (i) n−1 ). Once the ensembles are generated, the sample covariances of the errors are computed as follows: where Then, the Kalman gain is approximated by: and the ensembles of the state vector, {x (i) n|n−1 }, are computed as where v (i) n are samples obtained from Gaussian distribution with mean y n and covariance Q w . Once the ensembles of the state vector are computed, the estimate of the state vector is obtained by taking the averages of ensembles as follows:

State-Space Model of a Synthetic ECG
In the work reported in, McSharry et al. proposed a dynamic model which consists of a set of nonlinear state equations to generate synthetic ECG signals in the Cartesian coordinate system [37]. Further, Sameni et al. transformed the model to a polar coordinate system and provided a convenient discrete-time mathematical model [38]. The model represents an ECG signal by a sum of five Gaussian functions, each corresponding to the five waves of an ECG signal, namely P, Q, R, S, T waves. The state vector of the dynamic model is defined by x k = [θ k , z k ] T , and the state equation is given by: where ∆θ i = (θ k − θ i )mod 2π is the phase increment, ∆ is the sampling period, η k is the state noise, ω is the angular velocity of the trajectory as it moves around the limit cycle, and α i , b i and θ i represent the amplitude, width, and center of the Gaussian functions of the five PQRST waves, respectively. The measurement vector is defined by where φ k is the observed phase representing the linear time wrapping of the R-R time interval into [0, 2π], and s k is the observed amplitude. The measurement equation is given by where u k and v k denote the measurement noises.

EnKF Based fECG Extraction Algorithm
Given the state-space model (13) and (14), the EnKF algorithm in Figure 1 is applied to filter out the mECG from aECG, assuming the remaining signal composed of fECG and noise is Gaussian distributed. The extracted mECG signal is then subtracted from aECG to obtain a noisy fECG signal. Finally, the EnKF algorithm is applied to the residual signal to denoise the fECG signal. Before running the extraction algorithm, the acquired aECG signals are processed to remove the baseline wander, the powerline, and ambient interferences. The baseline wander is removed using a lowpass filter, and a notch filter is used to suppress the powerline interference noise. We also applied the wavelet filtering and thresholding technique and compare [40]. This will be discussed in detail in later sections. Finally, the fetal QRS complex (fQRS) is detected using the Pan-Tompkin algorithm [32,38].

Data for Testing
Three different datasets were used to test the performance of the proposed algorithm and to compare with the EKF. They are the data obtained from the PhysioNet 2013 Challenge database, the same data with motion noise added, and our own data obtained in a recent pilot study with pregnant subjects.

Challenge Databank
The PhysioNet 2013 Challenge databank used in this work consists of 75 recordings, excluding a number of recordings (a33, a38, a47, a52, a71, and a74) that had inaccurate reference annotations. Each recording comprises four different abdominal signals. All sig- nals have been sampled at 1 kHz and recorded for 60 s. In each case, reference annotations marking the locations of each fetal QRS complex were produced, usually with reference to a direct fECG signal, acquired from a fetal scalp electrode. The reference annotations are produced by a team of experts manually [41]. . It propagates and updates the ensemble to track the distribution of . The state estimation is conducted by calculating the sample mean (red five-pointed-star) and covariance (red ellipse) of the ensemble.

Data for Testing
Three different datasets were used to test the performance of the proposed algorithm and to compare with the EKF. They are the data obtained from the PhysioNet 2013 Challenge database, the same data with motion noise added, and our own data obtained in a recent pilot study with pregnant subjects.

Challenge Databank
The PhysioNet 2013 Challenge databank used in this work consists of 75 recordings, excluding a number of recordings (a33, a38, a47, a52, a71, and a74) that had inaccurate reference annotations. Each recording comprises four different abdominal signals. All signals have been sampled at 1 kHz and recorded for 60 s. In each case, reference annotations marking the locations of each fetal QRS complex were produced, usually with reference to a direct fECG signal, acquired from a fetal scalp electrode. The reference annotations are produced by a team of experts manually [41].

Modified Signals with Motional Artifacts Added
The aECG recordings acquired in real-life settings would possess a variety of interferences, including motion artifacts. The online databank, however, was obtained in the

Modified Signals with Motional Artifacts Added
The aECG recordings acquired in real-life settings would possess a variety of interferences, including motion artifacts. The online databank, however, was obtained in the clinical setting, where motion noise was mostly non-existent since the subjects were in a resting position. We added realistic motion artifacts to the online databank. The realistic motion noise generation process is described in our previous paper [32], to mimic real-life scenarios. In this experiment, the ECG data were recorded from a healthy subject during different types of activities. Then, the normalization had been used to reinsure aECG and the motion noise have realistic amplitudes. The recorded data were normalized between −1 and 1. Then, the motion noise was extracted by using the EKF. Before adding the motion noise to the new aECG, the aECG data should be normalized with the same threshold. In the last step, the generated motion noises are added to the PhysioNet 2013 Challenge databank. Figure 2 illustrates the modified data with motion noise process. different types of activities. Then, the normalization had been used to reinsure aECG and the motion noise have realistic amplitudes. The recorded data were normalized between −1 and 1. Then, the motion noise was extracted by using the EKF. Before adding the motion noise to the new aECG, the aECG data should be normalized with the same threshold. In the last step, the generated motion noises are added to the PhysioNet 2013 Challenge databank. Figure 2 illustrates the modified data with motion noise process.

Our Human Data
We developed the gen-2 'fetal monitoring patch' containing non-contact electrodes (NCEs), electronics, and secure communication with a smart device via Bluetooth Low Energy (BLE) [12]. The compact patch is ~4 inches long and unobtrusive, thus, it can be integrated on or inside maternity garment. The gen-2 patch has one single NCE channel to collect the aECG of the pregnant subject. The collected data are sent to an Android app connected to a cloud server for analytics. The system was validated on 10 pregnant women between 28 and 34 weeks of gestation in the UCI Medical Center. Each subject was resting on a chair during recording, and a maternity belt with the fECG patch attached is worn so that the NCEs are located on the abdominal area below the navel. The data were then collected in 5 min for each posture. More details are described in our previous paper [42]. All clinical recordings were collected anonymously under Institutional Review Board (IRB) approval #2020-6342 at the University of California, Irvine (UCI).

Comparison Criteria
The performance of the extraction methods is assessed by comparing the beat-to-beat length of the extracted fECG QRS complex and the corresponding annotated data. According to the American National Standards Institute/Association for the Advancement of Medical Instrumentation (ANSI/AAMI) guideline, sensitivity (SE), positive predictive value (PPV), and the accuracy measure (F1 score) which is the harmonic mean of PPV and SE, were used for assessment. These statistical indices are computed as follows where TP, FP, and FN are the number of true positive, false positive, and false negative, respectively.

Our Human Data
We developed the gen-2 'fetal monitoring patch' containing non-contact electrodes (NCEs), electronics, and secure communication with a smart device via Bluetooth Low Energy (BLE) [12]. The compact patch is~4 inches long and unobtrusive, thus, it can be integrated on or inside maternity garment. The gen-2 patch has one single NCE channel to collect the aECG of the pregnant subject. The collected data are sent to an Android app connected to a cloud server for analytics. The system was validated on 10 pregnant women between 28 and 34 weeks of gestation in the UCI Medical Center. Each subject was resting on a chair during recording, and a maternity belt with the fECG patch attached is worn so that the NCEs are located on the abdominal area below the navel. The data were then collected in 5 min for each posture. More details are described in our previous paper [42]. All clinical recordings were collected anonymously under Institutional Review Board (IRB) approval #2020-6342 at the University of California, Irvine (UCI).

Comparison Criteria
The performance of the extraction methods is assessed by comparing the beat-tobeat length of the extracted fECG QRS complex and the corresponding annotated data.  (15) where TP, FP, and FN are the number of true positive, false positive, and false negative, respectively.

Results
The execution time of each algorithm is calculated from the start of the pre-processing to the end of the fECG R-peak detection. The EnKF with ensemble size of 70 has the same execution time such as the EKF. Hence, the proposed EnKF-based algorithm was run with an ensemble size of N = 70 for single-channel signals from the aforementioned datasets. For comparison purposes, the EKF-based algorithm using the same synthetic ECG parameters and datasets was also carried out. The three-dimensional trajectory generated from (13) consists of a unit-radius circular limit cycle that goes up and down when it approaches one of the P, Q, R, S, or T points (Figure 3). The projection of these trajectory points on the z-axis gives a synthetic ECG signal. Figure 3 shows plots of the ECG signals versus the assigned phases in polar coordinates on the unit-radius circle. The figure depicts a typical phase-wrapped aECG signal (a), EKF extracted mECG signal (b), and EnKF extracted mECG signal (c) plotted using a sample signal taken from the PhysioNet database. an ensemble size of N = 70 for single-channel signals from the aforementioned datasets. For comparison purposes, the EKF-based algorithm using the same synthetic ECG parameters and datasets was also carried out. The three-dimensional trajectory generated from (13) consists of a unit-radius circular limit cycle that goes up and down when it approaches one of the P, Q, R, S, or T points (Figure 3). The projection of these trajectory points on the z-axis gives a synthetic ECG signal. Figure 3 shows plots of the ECG signals versus the assigned phases in polar coordinates on the unit-radius circle. The figure depicts a typical phase-wrapped aECG signal (a), EKF extracted mECG signal (b), and EnKF extracted mECG signal (c) plotted using a sample signal taken from the PhysioNet database. Examples of two abdominal ECG signals of records "a01" and "a03", and fetal signals extracted using the EKF and EnKF, are shown in Figure 4. Panels a and d are the original signals, while b-e and c-f are fQRS extracted using the EKF and EnKF approaches, respectively. The fQRS annotation is shown in an orange asterisk (*). The red arrows show the places that fetal QRS was wrongly detected. The blue arrows show the missing fetal QRS. It can be seen that the detected fQRS (fetal R-peaks) follow the annotated ones with high accuracy. It is worth noting that in "a01" (Figure 4d), the fetal QRS complexes are reversed due to electrode placement, but it did not affect the extraction algorithms. It should be emphasized that in the case of overlapping of the fetal QRS and maternal QRS, the EnKF algorithm still gives favorable results. Comparison between Figure 4b,c show that the EKF failed when maternal and fetal QRS complexes overlap in time (e.g., at t ~= 21.8 s and t ~= 23.7 s) for 'a03', while our EnKF method still performed successfully. It can be also seen that the EnKF functions reasonably well with the presence of noise ((e.g., at t ~= 19.7 s for 'a01', Figure 4d-f). Examples of two abdominal ECG signals of records "a01" and "a03", and fetal signals extracted using the EKF and EnKF, are shown in Figure 4. Panels a and d are the original signals, while b-e and c-f are fQRS extracted using the EKF and EnKF approaches, respectively. The fQRS annotation is shown in an orange asterisk (*). The red arrows show the places that fetal QRS was wrongly detected. The blue arrows show the missing fetal QRS. It can be seen that the detected fQRS (fetal R-peaks) follow the annotated ones with high accuracy. It is worth noting that in "a01" (Figure 4d), the fetal QRS complexes are reversed due to electrode placement, but it did not affect the extraction algorithms. It should be emphasized that in the case of overlapping of the fetal QRS and maternal QRS, the EnKF algorithm still gives favorable results. Comparison between Figure 4b,c show that the EKF failed when maternal and fetal QRS complexes overlap in time (e.g., at t ∼ = 21.8 s and t ∼ = 23.7 s) for 'a03', while our EnKF method still performed successfully. It can be also seen that the EnKF functions reasonably well with the presence of noise ((e.g., at t ∼ = 19.7 s for 'a01', Figure 4d-f).  fECG extracted using EKF of record "a03"; (c) fECG extracted using EnKF of record "a03"; (d) Abdominal ECG (raw data) of record "a01" with reversed fetal QRS complexes; (e) fECG extracted using EKF of record "a01"; (f) fECG extracted using EnKF of record "a01". Figure 5 depicts the fECG extraction results from the modified PhysioNet data with added motion artifacts. As seen, the fECG extracted by the EKF was incorrect and its F1 score was reduced significantly (average F1 = 78, see Table 1). fQRS complexes extracted by the EnKF, however, are still visible in most time points, yielding a favorable F1 score of 89. . fECG extraction using the EKF and EnKF with the PhysioNet data. The fQRS annotation is shown in orange asterisk (*). The red arrows show the places that fetal QRS was wrongly detected. The blue arrows show the missing fetal QRS. (a) Abdominal ECG (raw data) of record "a03"; (b) fECG extracted using EKF of record "a03"; (c) fECG extracted using EnKF of record "a03"; (d) Abdominal ECG (raw data) of record "a01" with reversed fetal QRS complexes; (e) fECG extracted using EKF of record "a01"; (f) fECG extracted using EnKF of record "a01". Figure 5 depicts the fECG extraction results from the modified PhysioNet data with added motion artifacts. As seen, the fECG extracted by the EKF was incorrect and its F1 score was reduced significantly (average F1 = 78, see Table 1). fQRS complexes extracted by the EnKF, however, are still visible in most time points, yielding a favorable F1 score of 89. Figure 4. fECG extraction using the EKF and EnKF with the PhysioNet data. The fQRS annotation is shown in orange asterisk (*). The red arrows show the places that fetal QRS was wrongly detected. The blue arrows show the missing fetal QRS. (a) Abdominal ECG (raw data) of record "a03"; (b) fECG extracted using EKF of record "a03"; (c) fECG extracted using EnKF of record "a03"; (d) Abdominal ECG (raw data) of record "a01" with reversed fetal QRS complexes; (e) fECG extracted using EKF of record "a01"; (f) fECG extracted using EnKF of record "a01". Figure 5 depicts the fECG extraction results from the modified PhysioNet data with added motion artifacts. As seen, the fECG extracted by the EKF was incorrect and its F1 score was reduced significantly (average F1 = 78, see Table 1). fQRS complexes extracted by the EnKF, however, are still visible in most time points, yielding a favorable F1 score of 89. Figure 5. fECG extraction using the EKF and EnKF with the motion artifacts added PhysioNet data. The fQRS annotation is shown in orange asterisk (*). The red arrows show the places that fetal QRS was wrongly detected. The blue arrows show the missing fetal QRS. (a) Abdominal ECG (raw data); (b) fECG extracted using EKF; and (c) fECG extracted using EnKF. Figure 6 shows the performance of the proposed algorithm, the EnKF, and the EKF on one sample from our clinical data. Panels a-c show the original and extracted signals, while panels d-f are signals with extensive preprocessing. Specifically, in Figure 6d, a lowpass filter (with a cut-off frequency of 1 Hz), a notch filter, and a Wavelet filter (a 10 level 1-D stationary wavelet decomposition with Coiflet mother wavelet) were used to suppress the background noise and artifacts. The signal extracted from the preprocessed data with the EKF and EnKF algorithms are shown in Figure 6e,f. In Figure 6f, we also further notice that fetal characteristic waves, such as P and T waves, may be conserved. This judgment is also strengthened by the EnKF extraction carried out on the original   Figure 6 shows the performance of the proposed algorithm, the EnKF, and the EKF on one sample from our clinical data. Panels a-c show the original and extracted signals, while panels d-f are signals with extensive preprocessing. Specifically, in Figure 6d, a lowpass filter (with a cut-off frequency of 1 Hz), a notch filter, and a Wavelet filter (a 10 level 1-D stationary wavelet decomposition with Coiflet mother wavelet) were used to suppress the background noise and artifacts. The signal extracted from the preprocessed data with the EKF and EnKF algorithms are shown in Figure 6e,f. In Figure 6f, we also further notice that fetal characteristic waves, such as P and T waves, may be conserved. This judgment is also strengthened by the EnKF extraction carried out on the original PhysioNet database in Figure 4f, where the conserved features are probably fECG waves. Table 1 presents the average F1, average PPV, and average SE results in our own clinical aECG records. The average F1, PPV, and SE indices of our proposed EnKF are 94.3%, 100%, and 89.2%; while those using the EKF method are 82.3%, 71.4%, and 100%, respectively. Table 1 shows the average F1 scores, the average PPE, and the average SV of the performance of the EnKF and EKF algorithms. These statistical indices are computed by determining the accuracy, TP, FP and FN, of the locations of the R-peaks obtained by the EKF and EnKF algorithms against the reference annotations. The F1 score, PPV, and SE are computed using 68 one-minute aECG records from PhysioNet 2013 Challenge databank and our own clinical data. The results shown in Table 1 indicate that the EnKF method is reliable on its own. In all cases, the EnKF outperforms the EKF. Figure 4f, where the conserved features are probably fECG waves. Table 1 presents the average F1, average PPV, and average SE results in our own clinical aECG records. The average F1, PPV, and SE indices of our proposed EnKF are 94.3%, 100%, and 89.2%; while those using the EKF method are 82.3%, 71.4%, and 100%, respectively. Figure 6. fECG extraction using the EKF and EnKF with our own clinical data. The fQRS annotation is shown in orange asterisk (*). The red arrows show the places that fetal QRS was wrongly detected. The blue arrows show the missing fetal QRS. (a) Abdominal ECG (raw data), (b) fECG extracted from EKF and (c) fECG extracted from EnKF". (d) Abdominal ECG after Wavelet preprocessing; (e) fECG extracted using EKF of preprocessed data; and (f) fECG extracted using EnKF of preprocessed data. The inset shows additional peaks, likely P and T waves; they were conserved after EnKF extraction. Table 1 shows the average F1 scores, the average PPE, and the average SV of the performance of the EnKF and EKF algorithms. These statistical indices are computed by determining the accuracy, TP, FP and FN, of the locations of the R-peaks obtained by the EKF and EnKF algorithms against the reference annotations. The F1 score, PPV, and SE are computed using 68 one-minute aECG records from PhysioNet 2013 Challenge databank and our own clinical data. The results shown in Table 1 indicate that the EnKF method is reliable on its own. In all cases, the EnKF outperforms the EKF.

PhysioNet database in
Another parameter that should be taken into account is the computational complexity. The computational complexity of the EnKF algorithm is proportional to the number of ensembles used. In our simulations, we have observed that increasing the size of the ensemble by more than 70 does not improve the performance of the algorithm significantly. The F1 score obtained when the algorithm was run for ensemble sizes between 5 to 350 always remained in the range between 94.5% and 98.6% for all of the 68 aECG records obtained from the PhysioNet database. Another parameter that should be taken into account is the computational complexity. The computational complexity of the EnKF algorithm is proportional to the number of ensembles used. In our simulations, we have observed that increasing the size of the ensemble by more than 70 does not improve the performance of the algorithm significantly. The F1 score obtained when the algorithm was run for ensemble sizes between 5 to 350 always remained in the range between 94.5% and 98.6% for all of the 68 aECG records obtained from the PhysioNet database.

Discussion and Conclusions
The current fECG extraction methods from a single-channel signal are not robust when particularly (i) the fECG and mECG waveforms temporally overlap, and (ii) the amplitude of fECG is low compared to the noise level. Recently, an EKF-based algorithm was proposed for the extraction fECG from a single-channel aECG signal [43]. Generally, those EKF-based algorithms are less effective as at every instant of time, they approximate the posterior probability density of the parameter of interest by a Gaussian distribution. When the true posterior density is not Gaussian, sequential Monte Carlo (SMC) filtering methods show superior performance over EKF methods. Therefore, here, our proposed EnKF is an SMC-based method, for noninvasive extraction of fECG from the single-channel aECG. As described above, our EnKF-based algorithm exhibited robust performance when tested using public online data as well as our own clinical data. Ten records, including 20-min aECG from our own clinical database, 68 records, including 1-min aECGs in Phys-ioNet Challenge 2013, and the online data bank added realistic motion artifacts were used to assess the performance of the proposed method.
As can be seen in Figure 3, the EnKF-based algorithm follows the dynamics of ECG and thus suppresses the noise better than EKF. In the EKF method, since some of the fECG peaks were still incorporated in the estimated mECG (3b), the subtraction between the aECG and mECG will not produce the correct fECG. In contrast, the mECG estimated by EnKF is significantly clearer (3c), which explains the better performance in fECG extraction (Figures 4-6 and Table 1).
Our experiments proved that the EnKF-based algorithm is a robust method and has superior performance over the EKF for extracting fECG in various scenarios (Figures 4-6). Results support the expectations that the EnKF not only extracts the fetal signal when the fECG and mECG waveforms temporally overlap, but it also potentially extracts accurate fECG signal with characteristic waves such as P and T waves. Currently, there is no approach available to extract full-feature fECG, especially with aECG in daily life. Hence, our solution may hold the potential to revolutionize fetal monitoring as it can be used to diagnose potential congenital effects. Currently, this is carried out with genetic test and echocardiogram [44][45][46]. However, those need to be carried out in the clinics and they cannot provide continuous information of the fetal heart over a long period of time. For example, with our fetal ECG patch and this EnKF method, full-feature fECG can be obtained 24/7 in the home setting in the second and third trimesters of pregnancy, which provides information about not only fetal wellbeing and development but also any fetal cardiac anomalies.
Due to the nature of aECG acquisition, the signal is always accompanied by other bioelectric potentials, such as maternal muscle activity, fetal movement activity, and noise [33]. The results shown in Figures 4-6 suggested that, prior to fECG extraction, preprocessing is a critical step which needs to be optimized. Specifically, preprocessing may help extraction significantly as it removes noise components, but it may also eliminate low-amplitude precious components such as P and T waves. In Figure 6d, we can see that, when the aECG signal is somewhat over-processed, both EKF and EnKF perform reasonably well; the signal, however, may have lost its intrinsic peaks. The filtering scheme used here is a low-pass filter, a digital Notch filter, and the wavelet filtering and thresholding technique described in [40]. As P and T waves have similar frequency and amplitude with noise components, it is almost impossible to set up optimal thresholds to eliminate the noise and keep the desired waves in all cases. Therefore, we usually carry out this manually until the best performance is obtained. The challenge is the noise level varies from person to person, setting to setting, and even time to time within one person in one setting. To address this, we are planning to add an accelerometer to collect the motion noise as carried out in [12], and explore further algorithms to extract uterine contraction and fetal movement.
The fECG would provide valuable information that could help deliver better prenatal care as well as assisting clinicians in making appropriate and timely decisions during labor. In this context, our EnKF presented in this work may bring significant impacts as a robust and efficient algorithm deployed in compact wearables for fECG assessment in the home setting. Future work includes assessing the fECG patch and the EnKF algorithm in out-of-clinic settings. We are planning to first validate the performance with pregnant subjects during sleep overnight in their home first and then move forward to real-life trials for several days. We also plan to correlate with other physiological signals using machine learning and cloud computing to link the mother's activity with fetal development, such as the mother's sleep quality or diet [42,47]. In the future work, the fetal ECG monitoring in twins and multi-fetal prenatal should also be studied.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board of UC Irvine, protocol #2020-6342.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: Not applicable.