Quantiﬁcation of Feto-Maternal Heart Rate from Abdominal ECG Signal Using Empirical Mode Decomposition for Heart Rate Variability Analysis

: In this paper, a robust method of feto-maternal heart rate extraction from the non-invasive composite abdominal Electrocardiogram (aECG) signal is presented. The proposed method is based on the Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) method, in which a composite aECG signal is decomposed into its constituent frequency components called Intrinsic Mode Functions (IMFs) or simply “modes”, with better spectral separation. Decomposed IMFs are then selected manually according to probable maternal and fetal heart rate information and are processed further for quantiﬁcation of maternal and fetal heart rate and variability analysis. The proposed method was applied to aECG recordings collected from three different sources: (i) the PhysioNet (adfecgdb) database; (ii) the PhysioNet (nifecgdb) database; and (iii) synthetic aECG signal generated from mathematical modeling in the LabVIEW software environment. An overall sensitivity of 98.83%, positive diagnostic value of 97.97%, accuracy of 96.93% and performance index of 96.75% were obtained in the case of Maternal Heart Rate (MHR) quantiﬁcation, and an overall sensitivity of 98.13%, positive diagnostic value of 97.62%, accuracy of 95.91% and performance index of 95.69% were obtained in case of Fetal Heart Rate (FHR) quantiﬁcation. The obtained results conﬁrm that CEEMDAN is a very robust and accurate method for extraction of feto-maternal heart rate components from aECG signals. We also conclude that non-invasive aECG is an effective and reliable method for long-term FHR and MHR monitoring during pregnancy and labor. The requirement of manual intervention while selecting the probable maternal and fetal components from “ n ” number of decomposed modes limits the real-time application of the proposed methodology. This is due to the fact that the number of modes “ n ” produced by the CEEMDAN decomposition is unpredictable. However, the proposed methodology is well suited for applications where a small time-delay or offset in feto-maternal monitoring can be acceptable. In future, application-speciﬁc modiﬁcation of the CEEMDAN algorithm can be implemented to eliminate manual intervention completely and will be suitable for long-term feto-maternal monitoring.


Introduction
Electronic Fetal Monitoring (EFM) is an essential requirement when assessing the health state of the fetus and its overall growth and development inside the uterus during pregnancy and labor [1].In high-risk conditions like an Intra-Uterine Growth Retarded (IUGR) fetus or pregnancies complicated by pathologies like diabetes, infections, pre-eclampsia or placental abruption, etc., long-term fetal monitoring is required [2].Fetal Heart Rate (FHR) analysis using the Cardiotocograph (CTG) or by direct fetal scalp ECG are the most effective techniques for fetal well-being assessment [3].However, CTG is not suitable for long-term monitoring and requires highly experienced medical personnel for the interpretation of FHR traces.Fetal scalp Electrocardiogram (ECG) recordings are considered to be the gold standard in FHR monitoring, but can only be applied when the membrane has been ruptured [4].Considering the above limitations, a non-invasive alternative for FHR monitoring is abdominal ECG (aECG) recordings.In the aECG method, electrodes are placed on the mother's abdomen, which acquire a composite ECG signal, which consists of both fetal ECG (fECG) and maternal ECG (mECG) [5] along with other artifacts like muscle contraction noise, power line interferences, various electrical noise, fetal movement and uterine contractions [6].The interpretation of morphology, magnitude and frequency of fECG helps physicians to identify overall fetal health state [7]; however, the magnitude of fECG signals is usually overwhelmed by the magnitude of mECG signals, and significant frequency overlaps with other artifacts make the Signal to Noise Ratio (SNR) of the fECG signal very low.Therefore, advanced signal processing techniques are required for the extraction of fECG from the aECG signal [8].
Researchers have been working in this direction and successfully evaluated various advanced signal processing methods.Kanjilal et al. proposed a method based on Singular Value Decomposition (SVD) for extraction of fECG from the single channel aECG signal by decomposing it into orthogonal modes and then selecting the most dominant periodic component by Singular Value Ratio (SVR) spectrum analysis [9].Based on the Artificial Neural Network (ANN) technique, Hasan et al. used the supervised Multilayer Perception (MLP) network [10]; Jia et al. used the adaptive linear neural network [11]; and Golzan et al. used the multi-layer perceptron neural network model for extraction of fECG [12].Based on the Adaptive Noise Cancellation (ANC) technique, Zhang et al. used a combination of SVD and Smooth Window (SW) techniques to build a reference signal in an ANC, which is done to eliminate the limitation of having a similar thoracic mECG waveform as that present in acquired the aECG waveform [13], and Zeng et al. proposed a Recursive Least Squares (RLS)-based ANC approach for fECG extraction by elimination of the mECG signal [14].Based on Principal Component Analysis (PCA), Rahmati et al. detected morphological ECG features like FHR, fetal R-peaks and RR interval from maternal aECG recordings collected from the PhysioNet/Computing in Cardiology Challenge 2013 (CinC2013) [15], Lipponen et al. also used the principal component model for removing the maternal ECG from aECG recordings for the purpose of fetal QRS detection [16], and Petrolis et al. used multistage PCA for the subtraction of mECG from aECG signals, leaving behind a signal with a higher energy of fECG components; they also tested their methodology on the database collected from PhysioNet/CinC2013 [17].Based on the Blind Source Separation (BSS) technique, He et al. proposed a method of mECG and fECG separation from a single channel abdominal signal.They also used the Empirical Mode Decomposition (EMD) method to decompose the collected single channel signal into multiple Intrinsic Mode Function (IMF) mapping single channels [18].Ghazdali et al. proposed an approach for denoising of the observed ECG signals using a Bilateral Total Variation (BTV) filter and then minimizing the Kullback-Leibler divergence [19] between copula densities to separate the fECG signal from the mECG [20].Sato et al. proposed a modified blind source separation with the reference signal (BSSR) for fECG extraction by cancellation of the mECG component [21].Zarzoso et al. tested three methods for fECG extraction PCA, Higher-Order Singular-Value Decomposition (HOSVD) and Higher-Order Eigenvalue Decomposition (HOEVD) and found that HOEVD was the best amongst them [22]; in a consecutive study, they implemented higher-order statistics and Widrow's multireference ANC approach for more robust BSS performance in the fECG extraction problem [23].Based on the Adaptive Neuro-Fuzzy Inference System (ANFIS), Assaleh et al. demonstrated the use of the ANFIS network to identify the nonlinear relationship that exists between the mECG signal acquired from thorax and the mECG component derived from the aECG signal, thus helping in the fECG extraction process by subtracting the aligned version of the mECG signal from the aECG signal [24].Panigrahy et al. implemented ANFIS along with the Differential Evolution (DE) algorithm and the Extended Kalman Smoother (EKS) for fECG extraction, requiring some degree of operator interaction for assigning the parameters and the EKS state equations [25].Su et al. applied a nonlinear time-frequency analysis technique called de-shape Short Time Fourier Transform (STFT) on PhysioNet (adfecgdb) database and CinC2013 database for fECG extraction [26].Almeida et al. validated the Wavelet Transform (WT)-based QRS algorithm for localization of the fetal QRS complex in aECG signals acquired from the PhysioNet database [27].
Recently, researchers have been involved in techniques like Empirical Mode Decomposition (EMD) for denoising of ECG signals [28,29], baseline wander removal [30], Power Line Interference (PLI) removal [31][32][33], QRS complex detection [34,35], extraction of the respiratory signal from ECG recordings [36,37], etc. EMD has emerged to be an essential tool for decomposing a nonlinear and non-stationary signal into its constituent frequency components with higher temporal and spatial resolution as compared to Fourier-based and wavelet-based time-frequency methods.Several improvements in basic EMD have been proposed like Ensemble Empirical Mode Decomposition (EEMD) [38] and Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) [39] for minimizing the limitations of basic EMD like the "end-effect" and "mode-mixing".However, very limited research has been done on the topic of the extraction of fetal heart rate information from the aECG signal using EMD [40][41][42].
In this paper, we propose the CEEMDAN-based method for reliable extraction of both maternal and fetal heart rate information from the composite aECG signal.The aECG signal datasets used in this research are acquired from three different publicly available online sources.The CEEMDAN technique decomposes the composite aECG signal into a finite number of components separated by instantaneous frequencies called modes [43,44].Careful selection of the modes according to both maternal and fetal components and post-processing for valid R-peak detection results in the extraction of FHR, MHR and their Heart Rate Variability (HRV).Feto-maternal HRV analysis helps the clinician to predict the overall health state of the fetus along with the cardiac autonomic regulation of the fetal heart.
The remainder of the paper is structured as follows: In Section 2, the review of the EMD technique and its improved versions like EEMD and CEEMDAN is presented.The datasets and proposed methodology are described in Section 3. In Section 4, experimental results are presented.In Section 5, the discussion about the results is presented, and finally, in Section 6, the main conclusions are pointed out.

Empirical Mode Decomposition
EMD is an adaptive signal processing method developed by Huang et al. in 1998 [45], which is used for decomposing a non-linear and non-stationary temporal signal into a finite set of Intrinsic Mode Functions (IMFs).IMFs or "modes" are a set of Amplitude-and Frequency-Modulated (AM-FM) oscillating components intrinsic to the signal, as shown in Equation (1), and are the bases of the decomposition.
where the phase φ k (t) is a non-decreasing function, φ k (t) ≥ 0, the envelope is non-negative A k (t) ≥ 0 and both the envelope A k (t) and the instantaneous frequency ω k (t) = φ k (t) vary much slower than the phase φ k (t) [46,47].
An IMF should satisfy two basic conditions: 1.The number of extrema (minima and maxima) and the number of zero crossings must be the same or differ at most by one 2.At any point, the mean value of the envelope defined by local extrema (minima and maxima) should be zero The basic principle (also known as the shifting process) required for decomposing a signal x (t) into its intrinsic oscillatory components (IMFs) is as follows: 1.The first step is to find all the local maxima and minima of the original signal x (t). 2. Use cubic spline interpolation to connect all the identified maxima and minima of the original signal x (t) to form an upper envelope e up (t) and a lower envelope e low (t).3. The mean m 1 (t) of the upper and lower envelope is calculated: 4. The difference d 1 (t) between the original signal x (t) and the mean m 1 (t) is calculated: 5.If the difference d 1 (t) satisfies the two criteria of IMF, then it is considered as the first IMF component of the original signal x (t) and defined as c 1 (t).The first IMF c 1 (t) contains the highest frequency component present in the original signal x (t).6.If the difference d 1 (t) does not satisfies the two criteria of IMF, then the shifting process described in Steps ( 1)-( 4) is repeated on d 1 (t) considering it as the original signal.7. The residue is calculated as r 1 (t) = x (t) − c 1 (t), which is treated as the original signal for the next iteration for higher order IMF calculation.8. Repeat the above process "n" number of times to extract "n" IMFs along with the final residue r n (t).9.The shifting process is terminated when r n (t) becomes a monotonic function from which no more IMFs can be extracted.Certain boundary condition like the Normalized Standard Deviation (NSD) can be imposed for terminating the shifting process.
when the value of NSD is less than a particular threshold "T", then the decomposition process is terminated.In this scenario, the residual r n (t) is either a constant, or a monotonic slope, or a function with only one extremum [48].10.At the end of the decomposition process, the original signal x (t) is represented as: where c n (t) is the "n-th" order IMF and r n (t) is the final residue or mean trend of the original signal x (t).By this convention, as the order of the IMF increases, the frequency content of the IMF decreases.A flowchart of the EMD decomposition process is shown in Figure 1.

Ensemble Empirical Mode Decomposition
A major drawback with EMD is frequent "mode-mixing", which occurs when a specific signal may not be separated into the same IMFs every time or when a signal of the same frequency is obtained from different modes.To solve this problem, Wu et al. in 2009 proposed a noise-assisted data analysis method, called EEMD [38].In EEMD, a white noise series is added to the original data.Then, this signal is decomposed into IMFs repeatedly with different white noise series each time and obtaining the (ensemble) means of corresponding IMFs of the decompositions as the final result.

Complete Ensemble Empirical Mode Decomposition with Adaptive Noise
EEMD successfully addresses the problem of mode-mixing by adding white noise into the signal, but this also leads to a problem that the added noise cannot be eliminated fully and may produce additional spurious modes.To resolve this problem, Torres et al. in 2011 proposed the CEEMDAN algorithm, which provides good spectral separation of the modes [39].Hence, it gives an exact reconstruction of the original signal with a lower computational cost.The following are the steps involved in the decomposition process by the CEEMDAN algorithm: 1. Decompose signal x(t) + ω 0 ε i (t) to obtain the first mode by using EMD algorithm.
where ω 0 is the ratio of the standard deviation of the added white noise, ε i (t) is the white noise with unit variance under the condition of the i-th ensemble number, N is the total ensemble number and c i 1 (t) is the first mode obtained after EMD decomposition.2. Compute the difference signal: to obtain the first mode, and define the second mode by: here, E 1 and ε i (t) stand for a function to extract the first IMF decomposed by EMD and the white noise with unit variance, respectively. 4. For k = 2, ..., K, calculate the k-th residue, and obtain the first mode.Define the (k + 1)-th mode as follows: where E k (•) is a function to extract the k-th IMF decomposed by EMD. 5. Repeat Step (4) until the residue contains no more than two extremes.The residual mode is then defined as: therefore, the signal x (t) can be expressed as follows:

Example
A typical noisy heart pulse signal and its resulting IMFs after CEEMDAN decomposition are shown in Figure 2.
This noisy heart pulse signal was obtained via an infra-red Photoplethysmograph (PPG) sensor (Heartbeat Pulse Sensor-1260, Sunrom Electronics, Ahmedabad, Gujarat, India) placed on the finger tip of a healthy subject, and the sampling rate was fixed at 1000 Hz.The PPG sensor was connected with a National Instruments (NI) myRIO-1900 data acquisition device (myRIO Student Embedded Device Model 1900 with bundled LabVIEW 2015, National Instruments Corporation, Austin, TX, USA) which captures the PPG recording and saves it in a PC (running on Windows 7 operating system with intel core i3 processor and 4GB of RAM memory) installed with NI LabVIEW software (2015 licensed version, National Instruments Corporation, Austin, TX, USA).
A segment of recorded PPG data with 4000 samples is shown in Figure 3 (Top), and it can be seen that the recorded PPG signal is full of noise due to power line interference and other ambient noise.Applying CEEMDAN to the noisy PPG signal and upon observing the obtained IMFs in Figure 2, it can be ascertained that most of the high frequency noise components are present in lower order IMFs, and the clinically important heart pulse information is present in higher order IMF 9 and IMF 10 (approximately).Therefore, the addition of the oscillating components present in IMF 9 and IMF 10 results in a denoised version of the PPG signal, as shown in Figure 3 (Bottom).Heart rate information from the denoised PPG signal can easily be extracted by the application of a peak detection algorithm.This proves that CEEMDAN is a very robust method for denoising a physiological signal without much loss of clinical information.

Data Acquisition
To evaluate the performance of the proposed method for quantification of fetal and maternal heart rate information contained within the composite abdominal ECG signal, we have collected data from three different sources as mentioned below: The first dataset was collected from the PhysioNet [49] website entitled "Abdominal and Direct Fetal Electrocardiogram Database" (adfecgdb), which contains multichannel fECG recordings obtained from 5 different women in labor, between 38 and 41 weeks of gestation [50][51][52].Each recording consists of a reference direct fECG signal registered directly from the fetal scalp and four differential signals acquired by placing Ag-AgCl electrodes around the navel region of the maternal abdomen surface.The researchers who contributed to this database state that these recordings are an excellent material for testing the credibility of new signal processing algorithms like maternal ECG suppression, extraction of fetal QRS complexes from abdominal recordings and evaluating the fetal heart rate variability for detecting overall fetal health state.Data were collected at a sampling rate of 1000 Hz with a resolution of 16 bits and were already digitally filtered for removal of power line interference of 50 Hz and baseline drift.Therefore, these datasets comprise an excellent choice for evaluating the application of the CEEMDAN technique for the separation of fetal and maternal heart beat information.Data from each subject consist of a 5-min recording of the direct fECG signal and four single channel recordings of 5 min each from the abdomen at a sampling rate of 1000 Hz, resulting in total of 3,00,000 samples from each subject per channel.In our experiment, we have selected the 4th channel from all the recordings (r0104, r0404, r0704, r0804 and r1004) and taken only one-tenth of the samples from the beginning, i.e., 30,000 samples (30 s of the recording), to reduce the computational complexity of the decomposition algorithm.
The second dataset was also collected from the PhysioNet [49] website entitled "Non-Invasive Fetal Electrocardiogram Database" (nifecgdb), which contains a series of 55 multichannel abdominal non-invasive fECG recordings, taken from a single subject from 21st week to 40th week of pregnancy.The researchers who contributed to this database state that these recordings may be very useful for testing signal separation algorithms.Each recording consists of 2 thoracic mECG signals and 3 or 4 aECG signals.The database states that the recordings were collected using a Ag-AgCl transducer, pre-amplified and filtered with a 50-Hz notch filter having a bandwidth of 0.01-100 Hz and a sampling rate of 1000 Hz at 16-bit resolution.We collected this database in CSV (Comma Separated Value) format via the PhysioBank ATM services.The whole database consisted of 55 records, but none of them contained fetal QRS annotations; only maternal QRS annotations were present.To determine the fetal QRS complex manually, we applied bandpass filtering (low cut-off frequency of approximately 60 Hz and high cut-off frequency of approximately 100 Hz), envelope generation and peak detection process on all of the recordings.We found clear and distinguished fetal R-peaks in only 6 of the records (nif0403, nif0404, nif0904, nif1003, nif1703, nif3403).Figure 4 shows the process of finding the fetal R-peaks in one of the records (nif0404).In the rest of the records, either the noise level was too high or the fetal QRS amplitude was too small.Although it was mentioned in the preamble of the "nifecgdb" database that the analog amplifier used to collect the recordings also includes a 50-Hz notch filter, in our analysis, we encountered the significant presence of power line interference at 50 Hz, which eventually distorts the fetal QRS complex.The fetal R-peaks obtained from the above process were considered as reference fetal QRS annotations in our research, and their locations were used as reference to calculate the true positive cases after application of the CEEMDAN-based methodology.We have taken 30,000 samples from the beginning of the recording data (30 s of the signal) to reduce the computational complexity of the CEEMDAN decomposition algorithm.The third dataset was generated from a synthetic Bio-Physiological Signal Generator (BPSG), which is a type of Virtual Instrument (VI) created by the authors in the LabVIEW software environment and was inspired by McSharry's dynamical ECG model [53].The VI generates a synthetic abdominal ECG signal from a combination of synthetic maternal and fetal ECG signals with added power line interference, baseline wander and Gaussian white noise.Signals were synthesized at a sampling rate of 500 Hz with respiration (12 respiration per minute with an amplitude of 0.15 mV) induced baseline wander.A comparison between a real aECG signal acquired from PhysioNet and a synthetically-simulated aECG signal using mathematical equations given in (10)-( 13) is shown in Figure 5.The BPSG VI front panel contains various user-defined controls for varying the model parameters like MHR, FHR, the amplitude of white noise, maternal respiration rate, RR interval and heart rate standard deviation.The best selected model parameters for mECG and fECG modeling are shown in Tables 1 and 2, respectively.The maternal QRS annotations and fetal QRS annotations were known from the synthetic maternal and fetal ECG signals used to synthesize the abdominal ECG signal.We have taken 20,000 samples from the beginning of the synthetic signal (40 s of the signal) to reduce the computational complexity of the CEEMDAN decomposition algorithm.To check the dependency of the CEEMDAN decomposition algorithm on feto-maternal heart rates, we have generated 5 synthetic abdominal signals (syn0101, syn0201, syn0301, syn0401, syn0501) with variable maternal and fetal heart rates, as shown in Table 3, and used these signals in our research.
FHR(GA) = −0.065GA+ 152.503,R 2 = 0.0016 ( 12) where a i , b i and c i are the amplitude, width and center of the Gaussian term, respectively.f 2 is the respiratory frequency, which is used to create baseline wander, and A is the amplitude of the respiratory sinusoid.The effect of gestation age (GA) and uterine contraction (UC) on FHR are also implemented in the simulator.

Proposed Methodology
The steps present in the proposed methodology are discussed below, and a flowchart of the process is shown in Figure 6.

CEEMDAN
Abdominal ECG datasets acquired from PhysioNet (sampled at 1000 Hz), as well as datasets synthetically generated by BPSG (sampled at 500 Hz) were decomposed into their respective IMF components using the CEEMDAN method in the LabVIEW software environment.Figure 7 shows the result of the CEEMDAN process.

Probable IMF Selection
The CEEMDAN process decomposes the aECG signal into various high and low frequency components.It can be observed from Figure 7 that most of the high frequency components are present in lower order IMFs and most of the low frequency components are present in higher order IMFs.Application of CEEMDAN to the PhysioNet dataset yields 12 IMF components, while that of synthetic dataset yields 11 IMF components.After careful examination, we have selected some lower order IMFs as probable candidates for containing most of the maternal heart rate information and selected some higher order IMFs as probable candidates for containing most of the fetal heart rate information from the decomposed datasets.In our research, the step of manual intervention for the selection of IMFs is required because we have performed our research on recordings collected from a variety of databases.Some recordings were collected from different subjects (adfecgdb), some at different times from a single subject (nifecgdb) and some at variable maternal and fetal heart rates (synthetic signals).Hence, to deal with this, we have introduced an additional step of IMF selection.The summation of the selected IMFs for maternal and fetal heart rate information was done to form a composite signal for further processing as shown in Figures 8 and 9.

Downsampling
Downsampling of the selected maternal and fetal IMF components was done to reduce the number of samples and to speed up the computation of heart rate.Proper precaution must be taken before downsampling the signal to avoid the aliasing affect.The addition of an anti-aliasing filter before downsampling prevents undesired signals or noise to alias into the desired signal band.Downsampling is an optional step in our proposed method.

Bandpass Filtering
Filtering a signal introduces a delay, i.e., the output signal is shifted in time with respect to the input signal.Therefore, a smart filter subVI is constructed in LabVIEW, which consists of two Butterworth filters, which completely eliminates the time shift in the output filtered signal.This smart filter is used for various low-pass and bandpass filtering operations required in our proposed work.The frequency of the aECG signal ranges approximately from 0.1-150 Hz [54,55].In order to remove unwanted frequencies that lie outside this range and preserve the frequencies of interest, we have applied bandpass filtering.

Envelope Generation
The envelope of the band passed signal is extracted to find out the actual periodicity of the heart beat.We have examined two different methods of envelope generation: the first uses a squaring and low pass fliter (LPF) method, and the second uses Hilbert transform and the LPF method.We found squaring and the LPF method to be slightly more accurate and to contain less ripples.Hence, we have selected this method for envelope generation.Generation of the envelope is essential because the envelope represents the energy of the signal at that instant of time.

Peak Detection
Peak detection is one of the essential time-domain analyses performed in signal processing.In LabVIEW, we have examined both the wavelet-based peak detection method and the quadratic fit and interpolation-based peak detection method.Both methods serve this purpose very efficiently by defining the proper threshold value.We have selected the quadratic fit and interpolation-based method because its subVI contains output for the number of peaks, location of peaks and amplitude of peaks.Accurate information about the location of the detected R-peaks is necessary for the identification of true positive cases.

Heart Rate Calculation
PhysioNet "adfecgdb" dataset contains the direct fECG signal acquired from fetal scalp along with its respective fetal QRS annotations.As the PhysioNet "nifecgdb" dataset does not contain any information about fetal QRS, we have manually extracted fetal QRS annotations by the process already mentioned in Section 3.1.These fetal QRS annotations are considered as a reference in our proposed method for identifying the true positive cases in heart rate detection.The number of peaks detected by our proposed method was directly compared with reference fetal QRS annotations for calculating various performance measures of the proposed methodology.

Heart Rate Variability Analysis
Heart Rate Variability (HRV) is a bio-physiological phenomenon of variation in the time interval between successive heart beats (RR interval) primarily due to two fluctuations, Respiratory Sinus Arrhythmia (RSA wave) and baroreceptor reflexes (Mayer wave).HRV analysis is a method to estimate autonomic regulation of the heart through quantification of these fluctuations.The frequency domain analysis shows that the low frequency band (0.04-0.15 Hz) includes physiologic oscillations associated with Mayer waves, and the high frequency band (0.15-0.40 Hz) is associated with RSA.The powers in these bands reflects the overall autonomic function of that individual.

Experiment and Results
We have used the PhysioNet aECG dataset, as well as the synthetic aECG dataset to evaluate the performance of CEEMDAN for the process of maternal and fetal heart rate extraction from a composite aECG signal.The performance of the proposed method was evaluated by the following measures: where FN (False Negative) indicates that the proposed method failed to detect a real heart beat, FP (False Positive) indicates that the proposed method detects a heart beat when their is no actual beat present, TP (True Positive) indicates that the real heart beat was detected accurately and N is the total number of the actual QRS complexes present in the reference signal.The F1 score is calculated based on sensitivity and positive diagnostic value and varies between [0,1].A value of "1" represents the optimal recognition capability of the proposed method, and a value of "0" represents that the proposed method is incapable of recognition.Tables 4 and 5 present the values of performance measures obtained after applying the proposed CEEMDAN-based feto-maternal heart rate extraction algorithm on all collected aECG datasets.For PhysioNet datasets, we assumed that a maternal R-wave and a fetal R-wave are considered to be True Positives (TP) if they are detected within ±10 ms from the reference peak.At a sampling rate of 1000 Hz, ±10 ms corresponds to a difference of 10 samples between the reference peak and the detected peaks.For the synthetic database, we have considered a TP case when a peak is detected within a difference of five samples between the reference peak and the detected peaks.The quality of feto-maternal heart rate separation depends on the overlap between their frequency spectrum in the composite abdominal ECG signal.Less frequency overlap means that it is easier to extract feto-maternal heart rate information from aECG signal using decomposition techniques like CEEMDAN.A comparison between the power spectral density of the aECG, mECG and fECG signals is shown in Figure 10, and the difference between the magnitude of power contained within the mECG signal and within the fECG signal can clearly be observed.In the case of fetal heart rate detection, from among a total number of 727 maternal QRS complexes recognized in the reference signals, 718 complexes were detected accurately, only nine complexes (1.24%) were missed, whereas the number of detected false complexes was only 16 (2.2%).The mean sensitivity of 98.83%, positive diagnostic value of 97.97%, accuracy of 96.93%, performance index of 96.75% and F1 score of 0.98 were obtained while detecting the maternal QRS complexes after the application of the proposed methodology as shown in Table 4.
In the case of fetal heart rate detection, from among a total number of 1231 fetal QRS complexes recognized in the reference signals, 1210 complexes were detected accurately, only 21 complexes (1.71%) were missed, whereas the number of detected false complexes was only 28 (2.28%).The mean sensitivity of 98.13%, positive diagnostic value of 97.62%, accuracy of 95.91%, performance index of 95.69% and F1 score of 0.98 were obtained as shown in Table 5.
We further studied the HRV analysis of the corresponding maternal and fetal heart rate intervals (RR intervals) determined from reference R-peaks and detected R-peaks after the CEEMDAN process.Table 6, presents the time domain results, while Tables 7 and 8 present the frequency domain results.Figure 11 shows the RR intervals plot (tachogram) for both MHR signals and FHR signals derived from the reference signal and the signal obtained after the CEEMDAN process.Figure 11A shows the RR-tachogram obtained from the reference mECG signal and extracted mECG signal and their corresponding Fast Fourier transform (FFT) spectrum and Autoregressive (AR) power spectrum plots.Various studies in the past have been done to evaluate the relationship between maternal HRV and fetal HRV during pregnancy-related pathologies like hypertension and pre-eclampsia [56,57].The influence of maternal Respiratory Sinus Arrhythmia (RSA) on maternal HRV and fetal autonomic nervous regulation in normal gestation is also investigated by researchers [57].
Figure 11B shows the RR-tachogram obtained from the reference fECG signal and extracted fECG signal and their corresponding FFT spectrum and Autoregressive (AR) power spectrum plots.The energy content in the frequency range of up to 0.2 Hz indicates the long-term variability, and the frequency range up to 1.5 Hz indicates the short-term variability.Long-term and short-term FHR variability are often related to fetal conditions [58,59].Fetuses also exhibit HRV analogous to an adult subject [60]; therefore, detailed analysis of fetal HRV can give sensitive information about the autonomic nervous system and cardiovascular system development of the fetus [60].

Discussion
From the results listed in Tables 4 and 5, the following important points can be made: 1.The selection of IMF for containing probable maternal and fetal frequency components depends on the subject, the recording conditions and the placement of the abdominal electrodes.In our research, manual intervention for IMF selection is required because we have applied the CEEMDAN methodology on the database acquired from various sources and acquired in different recording conditions.The proposed methodology can be automated for long-term patient-specific monitoring when there is no significant variation in recording conditions during the procedure.2. The selection of IMF does not depend on the variation of maternal or fetal heart rates.As seen from the results of the CEEMDAN decomposition of the synthetic database, where a combination of variable maternal and fetal heart rates (given in Table 3) were examined, still the selected IMF were similar (Mode 7 + 8 for the maternal component and Mode 4 for the fetal component).
3. In the PhysioNet "nifecgdb" dataset, the location of the fetal R-peaks obtained via the alternate process of band pass filtering (explained in Section 3.1) was assumed as reference fetal R-peaks for identifying TP cases.A mean accuracy of approximately 95% was obtained between the reference peaks and detected peaks, which indicates that the alternate method can be implemented for remaining records present in "nifecgdb" for which fetal QRS annotations are not defined.
Some limitations and challenges were also encountered during this research work and are listed below: 1.Although the CEEMDAN-based decomposition method provides exact reconstruction of the original signal with better frequency separation between the decomposed IMFs [44], like EMD and EEMD, it is a very time-consuming process due to its iterative nature [61,62] and is very difficult to implement in a real-time environment [63].2. The requirement of manual intervention for the purpose of feto-maternal IMF selection after CEEMDAN decomposition restricts our methodology to be implemented in real time.In the future, this limitation can be eliminated by introducing an automated IMF selection step in the algorithm.3. Implementing the methodology in the LabVIEW software environment instead of MATLAB makes the research a bit tedious.However, the graphical programming approach was found to be more suitable for our application.

Conclusions
Monitoring of fetal and maternal heart rate variability provides critical diagnostic information about the health state of the mother and her fetus.Long-term FHR monitoring is crucial in high-risk pregnancies, which can be accomplished only by a non-invasive method like aECG.Although CTG and ultrasound are extensively employed diagnostic tools for non-invasive monitoring of various physiological and morphological parameters of the fetus inside the uterus, it is not preferred for long-term bedside monitoring due to some degree of inconvenience to the patient and the need for frequent attendance by an experienced medical professional.
A CEEMDAN-based approach for extraction of both fetal and maternal heart rates from composite aECG recordings is presented in this research.The proposed method is based on the fact that in most of the cases, there is a frequency difference between maternal ECG components and fetal ECG components.CEEMDAN is a robust technique to decompose a composite signal into its respective frequency components called IMFs.CEEMDAN decomposition of the aECG signal produces IMFs in which the lower order IMFs contain higher frequency components (fetal ECG) and higher order IMFs contain lower frequency components (maternal ECG).Manually selecting the IMFs according to probable maternal and fetal components and taking the summation of theses components produces a signal for further post-processing.Post-processing involves downsampling, band pass filtering, envelope generation and peak detection.The detected peaks for both maternal and fetal cases were compared with the reference peaks for identifying TP, FP and FN cases.A TP case is considered when a detected R-peak is within ±10 ms from the reference R-peak.An overall sensitivity of 98.83%, a positive diagnostic value of 97.97%, an accuracy of 96.93% and a performance index of 96.75% were obtained in the case of MHR quantification, and an overall sensitivity of 98.13%, a positive diagnostic value of 97.62%, an accuracy of 95.91% and a performance index of 95.69% were obtained in the case of FHR quantification.The obtained results prove that the proposed CEEMDAN-based methodology is a very robust way of separating maternal and fetal components with good accuracy and sensitivity.It also confirms that an abdominal ECG recording is indeed a very useful technique for fetal heart rate extraction when long-term monitoring of a pregnant woman is required.Furthermore, the CEEMDAN method produces the best results when there is less overlap between the spectral components to be separated.Some limitations like manual selection of IMFs and time consumption during the decomposition process limit the usage of the proposed CEEMDAN-based methodology in a real-time environment.In the future, patient-specific modification of the proposed methodology can eliminate these limitations and can be used in real-time feto-maternal heart rate assessment applications.

Figure 4 .
Figure 4.The process of finding the reference fetal R-peaks: (A) aECG signal of 2000 samples from record nif0404; (B) band passed signal (60-100-Hz band); (C) detected fetal R-peaks, which are considered as reference fetal QRS annotations for the "nifecgdb" dataset.

Figure 5 .
Figure 5.An example of: (A) a noisy aECG signal obtained from PhysioNet (adfecgdb) database; (B) a noisy synthetic aECG signal simulated in LabVIEW using the Bio-Physiological Signal Generator (BPSG) virtual instrument.

Figure 6 .
Figure 6.Flowchart of the proposed methodology.HRV, Heart Rate Variability; LPF, Low Pass Filter.

Figure 8 .
Figure 8. PhysioNet aECG processing: (A) sum of the selected IMFs for probable MHR and FHR information; (B) downsampled and bandpass signal; (C) peak detection of the envelope generated using squaring and low-pass filtering (S+L); (D) peak detection of the envelope generated using Hilbert transform and low-pass filtering (H+L).

Figure 9 .
Figure 9. Synthetic aECG processing: (A) sum of the selected IMFs for probable MHR and FHR information; (B) downsampled and bandpass signal; (C) peak detection of the envelope generated using squaring and low-pass filtering (S+L); (D) peak detection of the envelope generated using Hilbert transform and low-pass filtering (H+L).
2 ms RR = Consecutive R-peak interval in milliseconds (ms); mReff.= Reference maternal RR interval; mExtr.= Extracted maternal RR interval; fReff.= Reference fetal RR interval; fExtr.= Extracted fetal RR interval; NN = normal R-peak to normal R-peak intervals; RMSSD = square root of the mean of the sum of the squares of differences between adjacent NN intervals; SDNN = standard deviation of all NN intervals; NN50 = number of pairs of adjacent NN intervals differing by more than 50 ms; pNN50 = the ratio of the number of NN50 intervals and the number of all NN intervals; TINN = triangular interpolation of NN interval.

Figure 11 .
Figure 11.HRV analysis of the RR tachogram obtained from the reference signal and the signal extracted using CEEMDAN decomposition: (A) analysis of the maternal ECG components; (B) analysis of the fetal ECG components.

Table 3 .
Combination of maternal and fetal heart rates used to generate synthetic (syn) abdominal signals.MHR, Maternal Heart Rate; FHR, Fetal Heart Rate.

Table 4 .
Results of the maternal R-peak detection and its various performance measures.Se, Sensitivity; PDV, Positive Diagnostic Value; PI, Performance Index; r, recording.The nomenclature of the record name indicates the database to which it belongs; the first 2 digits indicate the record number and the last 2 digits indicate the channel number (for example, nif0403 belongs to the "nifecgdb" database with Record Number 4 and selected abdominal Channel 3).

Table 5 .
Results of the fetal R-peak detection and its various performance measures.The nomenclature of the record name indicates the database to which it belongs; the first 2 digits indicate the record number, and the last 2 digits indicate the channel number (for example, r1004 belongs to the "adfecgdb" database with Record Number 10 and selected abdominal Channel 4).

Table 6 .
Time domain statistics.