A novel multi-step adaptive Kalman filtering method based on dynamic noise estimation for FECG extraction

Abstract

Objective

With the increase in women's childbearing age, the risk of fetal developmental abnormalities and fetal abortion is rising. The existing fetus monitoring methods based on Doppler ultrasound are inconvenient to use, require precise positioning of the fetal heart, and are particularly difficult to maintain for continuous long-term monitoring. Fetal electrocardiogram (FECG), as a very important physiological signal, can intuitively reflect the health status of the fetus. Affected by the strong noise and complex external components, the commonly used algorithm for extracting FECG from the abdominal electrocardiogram (AECG) of pregnant women cannot perform well enough.

Methods

Herein, we present a novel FECG extraction method based on three adaptive Kalman filters (KF). Regarding the first KF, to whiten the colored noise of the original AECG, the Expectation Maximization algorithm is used to iteratively solve the optimal parameters, and a new pseudo measurement variable named the “measurement time difference” is constructed. In addition, based on the residual vector $e_k$ and innovation $d_k$ , combined with the forgetting factor ${\alpha }_k$ , the measurement noise covariance matrix R and the process noise covariance matrix Q are adaptively updated, which can help the second and third KFs extract a more pure FECG.

Results

The proposed method demonstrates superior performance in the extraction accuracy and quality across various datasets, including 25 groups of AECG from the FECGSYN toolbox, 10 groups of clinical AECG from the ADFECGDB database, 10 groups of clinical AECG from the FECGDARHA database, and 3 groups of real AECG collected by the independently designed, portable, and low-cost AECG hardware system, when compared to 3 other commonly used algorithms.

Conclusion

Our method enables obstetricians to more accurately assess fetal physiological state, leading to more informed preventive measures and treatment plans, ultimately improving the health outcomes for both mothers and fetuses.

Introduction

With the changes of people's lifestyles, the childbearing age of women has been increasing and the proportion of pregnant women in older age has also been rising in recent years. ¹ Due to the gradual decline in physiological functions, the risk of accidental diseases that may occur in these women with advanced maternal age significantly increases, such as gestational hypertension, diabetes, preterm delivery, and fetal developmental abnormalities. Therefore, long-term portable monitoring of fetal health has become an important actual demand. Currently, fetal heart monitoring is a crucial method for obstetrics departments in hospitals to assess the intrauterine condition of fetuses. Fetal heart rate (FHR) curves can be real-time recorded through the ultrasound probes, by which fetal abnormalities such as tachycardia, bradycardia and intrauterine hypoxia can be promptly detected. Nevertheless, the above fetal heart monitoring method necessitates the identification of the fetal heart position prior to utilization. The complexity of operating these ultrasound devices poses a significant challenge for inexperienced pregnant women, which restricts their use in domestic environments. ² Furthermore, these devices tend to be relatively cumbersome and substantial, rendering them impractical for comfortable wear.

The abdominal electrocardiogram (AECG), comprising both the maternal electrocardiogram (MECG) and fetal electrocardiogram (FECG), is a critical physiological signal in perinatal medicine. It is essential for physicians to accurately evaluate the health status of the mother and fetus, thereby facilitating informed clinical diagnostic decisions and enhancing the safety of both mother and fetus. ³ Presently, the acquisition of AECG predominantly relies on large-scale hardware devices equipped with high-performance amplifiers and filters. However, these devices are associated with high costs, complex operation, and limited feasibility for widespread adoption among pregnant women. Consequently, low-cost and portable FECG devices hold substantial significance in the field of maternal and fetal health monitoring, particularly for long-term continuous home monitoring.^4,5 Nevertheless, these devices frequently encounter challenges such as the strong noise and considerable external interference, ⁶ which can substantially compromise signal quality and the accuracy of monitoring results. In conclusion, the development of methods capable of extracting high-quality FECG from low signal-to-noise ratio (SNR) AECG holds considerable medical meaning and clinical importance.

Extracting FECG from AECG is commonly achieved through three kinds of methodologies: blind source separation (BSS), template subtraction (TS), and filtering-based techniques. Within the realm of BSS, techniques such as the parallel linear predictor (PLP), principal component analysis (PCA), and independent component analysis (ICA) are employed. PLP ⁷ operates on the principle of linear prediction to estimate the target signal. However, it is prone to introducing errors during the reconstruction of non-stationary signals, which can compromise the quality of the target signal. By reducing the dimensionality of the original signal, PCA^8–10 effectively removes noise and redundant information from the original signal, facilitating efficient reconstruction of the target signal. Nonetheless, when the MECG and FECG are significantly overlapping, the reconstruction capabilities of PCA are often constrained. ICA^11,12 has the ability to decompose the mixed signals into several independent components. However, it is sensitive to initial conditions, which can yield unstable extraction outcomes. To mitigate these limitations, hybrid approaches have been proposed, such as combining ICA with Empirical Mode Decomposition ¹³ or adaptive filtering ¹⁴ to improve FECG extraction robustness. Nevertheless, it is also a challenging task to accurately identify which independent component corresponds to the target signal.

TS constructs the simple cycle template using the average of multiple cycles of the QRS complex, which template is then subtracted from the corresponding position of AECG at every occurrence of the QRS complex to extract the FECG. Due to the ease of use and straightforwardness, TS becomes popular in FECG extraction. Kanjilal ¹⁵ employs singular value decomposition (SVD) based on the spectrum of singular value ratio (SVR). The MECG and the FECG components are identified in terms of the appropriately configured matrices of the different decomposed modes. Vullings ¹⁶ presented a new technique for MECG and FECG extraction, which technique operated by the linear prediction of each mode from AECG. Hasan et al. ¹⁷ enhanced the accuracy of FECG extraction by optimizing the scaling factor to minimize the difference between MECG and its template. Martens et al. ¹⁸ adjusted the template based on the time sequence characteristics of MECG and FECG, aiming to improve the extraction effect of FECG. However, the performance of the aforementioned template-based algorithms heavily relies on the quality of the constructed template. Any slight mismatch between the template and the target signal can potentially be amplified during processing, thereby affecting the extraction effect of FECG.

Filtering techniques mainly include adaptive filtering, wavelet transform, and KF. Among them, adaptive filtering can dynamically adjust filter parameters to optimize the extraction effect of FECG, but it usually relies on additional reference signals and significantly slowdowns in convergence when dealing with rapidly time-varying or low SNR signals.^19,20Wavelet transform decomposes various frequency components from the original signal by expanding and compressing the mother wavelet, but its effectiveness highly depends on the chosen wavelet basis and the set thresholds.^21–23 Due to its streamlined architecture, low computational requirement, strong versatility, and flexibility, KF is widely used in dynamic signal filtering across various industries.^24–26 In recent years, many scholars have also explored the application of KF to the field of FECG extraction. Niknazar ²⁷ employed an extended Kalman filter (EKF) with fixed initial parameters for extracting FECG from single-channel AECG. However, this approach struggles to adapt to the time-varying characteristics of the target signal. Panigrahy ²⁸ enhanced the EKF by incorporating the differential evolution (DE) algorithm, yet DE's instability and propensity for local optimization can compromise prediction accuracy. Furthermore, Panigrahy ²⁹ also optimized the EKF using particle swarm optimization (PSO). However, the EKF based on PSO may suffer from premature convergence during the iterative process. Suganthy ³⁰ proposed an improvement to EKF through weighted least squares regression (WLSR), but the sensitivity of WLSR to outliers can affect the reliability of the filtering outcomes.

Although prior methods have improved FECG extraction, they often struggle with colored noise, strong nonstationarity, and fixed filter parameters. Bryson ³¹ introduced pseudo-measurement variables to model temporally colored noise, while Hesar ³² and Akhlaghi ³³ showed that EM-based parameter learning and forgetting-factor mechanisms can improve adaptivity in state-space filtering. However, these ideas have largely remained isolated and have not been systematically tailored to AECG or FECG extraction. Specifically, a unified framework tailored to AECG that jointly addresses temporally correlated noise, time-varying signal characteristics, and data-driven adaptation of filter parameters has yet to be established.

To address the limitations of existing approaches in terms of initial-parameter sensitivity, external noise interference, and rigid filter tuning, a new AECG-to-FECG extraction method is proposed (Figure 1). Herein, three main contributions are made: (i) an EM-integrated initial KF with a newly defined pseudo-measurement is applied to FECG extraction for the first time; (ii) a unified forgetting factor ${\alpha }_k$ enabling the simultaneous adaptive updating of both the measurement noise covariance matrix ( $R$ ) and the process noise covariance matrix ( $Q$ ) is proposed for the first time in the FECG extraction field; (iii) a coherent three-stage adaptive KF pipeline for AECG preprocessing, MECG extraction, and FECG extraction, is firstly organized, aiming to improve both morphological integrity and SNR of the extracted FECG.

Firstly, the original AECG is subjected to a second-order Butterworth bandpass filter with the frequency band of 3–40 Hz to remove the motion artifacts and baseline drift components. Then, to accurately track the dynamic changes of measurement noise and convert the colored noise to Gaussian noise, its temporal correlation is fully considered, and a new pseudo-measurement variable named “measurement time difference” is constructed. Besides, the EM algorithm is used to adaptively iterate and optimize the parameters of the first KF. Hence, the various external noise from the original AECG, excluding MECG and FECG, is more thoroughly removed.

To enhance the modeling accuracy of MECG and simultaneously reduce the residual signal in FECG, the second KF is designed using an MECG template fitted by nine Gaussian kernel functions in this article. Then, based on the residual vector $e_k$ and the innovation $d_k$ , the forgetting factor ${\alpha }_k$ is introduced to adjust the influence of the historical signal on the current signal. Besides, a low-pass filtering strategy is adopted to adaptively update R and Q, improving the filter's robustness in AECG with different SNR and achieving high-precision extraction effect of MECG. After the process mentioned in (2), the remaining signal mainly consists of FECG, a little residual MECG, and some residual noise.

To extract the target FECG from the residual signal of AECG processed by Steps (1) and (2) in order, the third KF is designed based on the FECG template fitted by seven Gaussian kernel functions. The adaptive updating strategy similar to that of Step (2) is also applied, where R and Q are also adaptively updated. Ultimately, the target FECG with high SNR and complete morphological characteristics is successfully extracted.

Figure 1.

The flow chart of the proposed method for FECG extraction in this article.

In summary, the proposed method for extracting FECG from AECG constructs the first adaptive KF based on the EM to suppress various noise components in the original AECG, remaining the mixed signal of MECG and FECG. Subsequently, combining the residual vector $e_k$ , the innovation $d_k$ with the forgetting factor ${\alpha }_k$ , the two adaptable KFs are designed, the R and Q matrix of which can be dynamically updated with the process noise and measurement noise's change. The second and third filters are successively applied to the remaining signal processed by the EM-based KF, yielding the pure and high-quality FECG. Our method enables obstetricians to more accurately assess fetal physiological state, leading to more informed preventive measures and treatment plans, ultimately improving the health outcomes for both mothers and fetuses.

Method

Dataset

Data acquisition

The extraction accuracy, SNR and morphological integrity of the target FECG are evaluated on 25 groups of AECG sourced from the FECGSYN toolbox,^34–36 10 groups of AECG clinically collected from the ADFECGDB, and three groups of AECG acquired by the hardware acquisition system designed in this article. As shown in Table 1, the above 25 groups of AECG consist of five independent sets generated by the FECGSYN toolbox, with each set being configured at five different SNR levels (3 dB, 5 dB, 6 dB, 7 dB, and 10 dB), all of which were sampled at 250 Hz, with a record duration of 5 minutes per record. The ADFECGDB database comprises the AECG from singleton pregnant women spanning gestational ages from 38 to 41 weeks, the original AECG signals of which are both sampled at the frequency of 1 kHz. Besides, record durations are variable across subjects and expert annotations of maternal and fetal R-peak positions are provided. In comparison, the FECGDARHA database ³⁷ primarily targets late-term pregnancies ranging from 32 to 42 weeks, which is sampled at a frequency of 500 Hz with a duration of approximately 20 minutes per record, and equally includes standardized expert annotations. Moreover, to preliminarily validate the accuracy of the proposed method for FHR calculation in fetal monitoring applications, this study was designed as a pilot experimental study. Three pregnant volunteers aged 33, 34, and 28 years, with gestational ages of 37, 38, and 40 weeks, respectively, were recruited at Liyuan Hospital, Tongji Medical College, Huazhong University of Science and Technology in March 2025. A common Doppler-based fetal monitoring device is used as the standard reference, and a custom-designed hardware system for recording AECG is employed synchronously to collect the original signal at 500 Hz on the abdominal surface of pregnant women with a duration of approximately 20 minutes, as depicted in Figure 2(a). The above system features an electrocardiogram module, a microprocessor module, a Bluetooth communication module, a power management module, and a lithium-ion battery module on the back of the board. Among them, ADS1299 is the core of the electrocardiogram module, whose electrodes need to be laid out on the surface of the pregnant woman's abdomen as shown in Figure 2(b). By the above hardware system, AECG can be collected with a common mode rejection ratio of 110 dB, which can ensure the quality of the original signal to a certain extent.

Figure 2.

The custom-designed system with the mainboard for recording AECG (a); the schematic diagram of electrode layout for ADS1299 (b).

Table 1.

The dataset characteristics used to validate the effectiveness of the proposed algorithm for FECG extraction.

Dataset	Sample size (Selected)	Sampling rate	Record duration	Source
FECGSYN	25	250 Hz	5 min	Synthetic
ADFECGDB	10	1 kHz	variable	Clinical
FECGDARHA	10	500 Hz	20 min	Clinical
Ours	3	500 Hz	20 min	Clinical

The evaluation index of the extraction effect

The American National Standards Institute (ANSI) and the Association for the Advancement of Medical Instrumentation (AAMI) have proposed a series of precision indicators to evaluate the extraction effect of MECG and FECG from AECG. As shown in Formulas (1) to (4).

$$\mathrm{SE}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\times 100\text{\%}}$$

(1)

$$\mathrm{PPV}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}\times 100\text{\%}}$$

(2)

$$\mathrm{ACC}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}\times 100\text{\%}}$$

(3)

$$\mathrm{F}1=2{\displaystyle \frac{\mathrm{PPV}\times \mathrm{SE}}{\mathrm{PPV}+\mathrm{SE}}}$$

(4)

In the above formula, TP is true positive, which refers to the correctly detected R peak. FN is false negative, which means the missed R peak. FP is false positive, which refers to the noise peak similar to R peak. SE represents sensitivity, PPV represents positive predictive value, and ACC represents accuracy. F1 is a statistic that is the harmonic average of SE and PPV. Combining the above four indexes, the accuracy and stability of the FECG extraction algorithm can be estimated comprehensively. In addition, SNR and morphological integrity of the extracted FECG are also two key indicators for quantitative evaluation on the algorithms’ performance, the specific computational methods of which are shown in Formulas (5) to (6).

$$\mathbf{S}\mathbf{N}\mathbf{R}=10{\log }_{10}{\displaystyle \frac{{\sum }_{\mathit{n}=1}^{\mathit{N}}\mathit{x}{(\mathit{n})}^2}{{\sum }_{\mathit{n}=1}^{\mathit{N}}{(\hat{\mathit{x}}(\mathit{n})-\mathit{x}(\mathit{n}))}^2}}$$

(5)

$$\mathbf{R}\mathbf{M}\mathbf{S}\mathbf{E}=\sqrt{{\displaystyle \frac{1}{\mathit{N}}\sum _{\mathit{n}=1}^{\mathit{N}}{(\hat{\mathit{x}}(\mathit{n})-\mathit{x}(\mathit{n}))}^2}}$$

(6)

where N is the length of the original signal, $\hat{x}(n)$ is the pure FECG, $x(n)$ is the extracted FECG, SNR is the signal-to-noise ratio, and RMSE is the root mean square error. The higher SNR represents the better quality of the extracted FECG. The lower RMSE, the better the morphological integrity of the extracted FECG.

AECG preprocessing

The principle of the classic KF

The classic KF can estimate unknown variables from the observation signal containing uncertain noise. ³⁸ In the discrete-time domain, the general state equation and measurement equation of KF are expressed as shown in Formulas (7) to (8).

$${\mathit{x}}_{\mathit{k}}={\mathit{A}}_{\mathit{k}}{\mathit{x}}_{\mathit{k}-1}+{\mathit{w}}_{\mathit{k}-1}$$

(7)

$${\mathit{y}}_{\mathit{k}}={\mathit{H}}_{\mathit{k}}{\mathit{x}}_{\mathit{k}}+{\mathit{v}}_{\mathit{k}}$$

(8)

In the above formula, $x_k$ , $y_k$ , $A_k$ , $H_k$ , $v_k$ , and $w_k$ are the state vector, measurement vector, state transition matrix, measurement matrix, measurement noise, and process noise in the $k$ th iteration, respectively. Among them, the measurement noise $v_k$ and process noise $w_k$ are modeled using normal white Gaussian noise. Besides, the process noise covariance matrix Q and the measurement noise covariance matrix R are used to describe the uncertainty during measurement.

The KF mainly includes the following two steps.

Prediction step:

$${\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}={\mathit{A}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}-1|\mathit{k}-1}$$

(9)

$${\mathit{P}}_{\mathit{k}|\mathit{k}-1}={\mathit{A}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}-1|\mathit{k}-1}{\mathit{A}}_{\mathit{k}}^{\mathit{T}}+{\mathit{Q}}_{\mathit{k}-1}$$

(10)

Update step:

$${\mathit{d}}_{\mathit{k}}={\mathit{y}}_{-}\mathit{H}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}$$

(11)

$${\mathit{S}}_{\mathit{k}}={\mathit{H}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}|\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k}}$$

(12)

$${\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}|\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}{\mathit{S}}_{\mathit{k}}^{-1}$$

(13)

$${\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}}={\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}{\mathit{d}}_{\mathit{k}}$$

(14)

$${\mathit{P}}_{\mathit{k}|\mathit{k}}={\mathit{P}}_{\mathit{k}|\mathit{k}-1}-{\mathit{K}}_{\mathit{k}}{\mathit{S}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$

(15)

where ${\hat{x}}_{k|k-1}$ and ${\hat{x}}_{k|k}$ represent the predicted and updated values of the state vector in the $k$ th iteration. $P_{k|k-1}$ and $P_{k|k}$ represent the predicted and updated values of the covariance matrix in the $k$ th iteration. Besides, $d_k$ , $K_k$ , and $S_k$ are each the innovation, the Kalman gain, and the covariance matrix of measurement prediction error.

The adaptive KF based on the EM algorithm

In order to first reduce the interference of various external noises and ensure the quality of the residual mixed signal of MECG and FECG, it is necessary to preprocess the original AECG. The backpropagation method of Rauch–Tung–Striebel (RTS) smoother is conducted to further smooth the output of the KF ³⁹ in this article. By working forward from the end of the sequence, the RTS smoother utilizes the output of the KFs’ prediction and update steps to optimize the past state vector and covariance matrix. As shown in Formulas (16) to (18).

$${\mathit{L}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}|\mathit{k}}{\mathit{A}}_{\mathit{k}+1}^{\mathit{T}}{\mathit{P}}_{\mathit{k}+1|\mathit{k}}^{-1}$$

(16)

$${\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}={\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}}+{\mathit{L}}_{\mathit{k}}({\hat{\mathit{x}}}_{\mathit{k}+1|\mathit{N}}-{\hat{\mathit{x}}}_{\mathit{k}+1|\mathit{k}})$$

(17)

$${\mathit{P}}_{\mathit{k}|\mathit{N}}={\mathit{P}}_{\mathit{k}|\mathit{k}}+{\mathit{L}}_{\mathit{k}}({\mathit{P}}_{\mathit{k}+1|\mathit{N}}-{\mathit{P}}_{\mathit{k}+1|\mathit{k}}){\mathit{L}}_{\mathit{k}}^{\mathit{T}}$$

(18)

In the above formula, ${\hat{x}}_{k|N}$ and $P_{k|N}$ are respectively the state vector and covariance matrix after the k th step of inverse smoothing, where N represents the length of the original signal.

According to Formulas (9) to (15), KF's performance is directly affected by the parameter set $\theta$ including the state transition matrix A, the measurement matrix $H$ , the process noise covariance matrix Q, the measurement noise covariance matrix R, and the covariance matrix P. Therefore, the EM algorithm is utilized in this article to iteratively refine the parameter set $\theta$ multiple times, leading to further optimization of the KF's filtering effect on the original AECG.

Specifically, the EM algorithm is a parameter estimation method for probability models containing some latent variables, which iteratively optimizes the models’ parameters by alternately performing E-step named expectation and M-step named maximization. By respectively calculating the latent variable's expectation and maximizing the likelihood function, the optimal parameters that maximize the likelihood of the observed signal can be figured out. ³² The process of EM is as follows.

For the observation signal $Y_N=\{y_1,y_2,\dots ,y_N\}$ , the maximum likelihood estimation of the probability parameter set $\theta$ is obtained by maximizing the marginal likelihood of $y_{1:N}$ with respect to $\theta$ .

$${\hat{\mathit{\theta }}}_{\mathit{M}\mathit{L}}=\arg \underset{\mathit{\theta }}{max}(\log {\mathit{p}}_{\mathit{\theta }}({\mathit{y}}_{1:\mathit{N}}))$$

(19)

Step E: Given the estimated parameters of $\theta$ after the i th iteration, the expectation of the logarithmic likelihood function $O(\theta ,{\theta }_i)=E[\log P_{{\theta }_i}(y_{1:N},x_{1:N}|y_{1:N})]$ is calculated as follows.

$$\mathit{E}[{\mathit{x}}_{\mathit{k}}|{\mathit{y}}_{1:\mathit{N}}]={\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}$$

(20)

$${\mathit{S}}_{1,1,\mathit{k}}=\mathit{E}[{\mathit{x}}_{\mathit{k}}{\mathit{x}}_{\mathit{k}}^{\mathit{T}}|{\mathit{y}}_{1:\mathit{N}}]={\mathit{P}}_{\mathit{k}|\mathit{N}}+{\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}^{\mathit{T}}$$

(21)

$${\mathit{S}}_{1,0,\mathit{k}}=\mathit{E}[{\mathit{x}}_{\mathit{k}}{\mathit{x}}_{\mathit{k}-1}^{\mathit{T}}|{\mathit{y}}_{1:\mathit{N}}]={\mathit{P}}_{\mathit{k},\mathit{k}-1|\mathit{N}}+{\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}{\hat{\mathit{x}}}_{\mathit{k}-1|\mathit{N}}^{\mathit{T}}$$

(22)

$${\mathit{S}}_{0,0,\mathit{k}}=\mathit{E}[{\mathit{x}}_{\mathit{k}-1}{\mathit{x}}_{\mathit{k}-1}^{\mathit{T}}|{\mathit{y}}_{1:\mathit{N}}]={\mathit{P}}_{\mathit{k}-1|\mathit{N}}+{\hat{\mathit{x}}}_{\mathit{k}-1|\mathit{N}}{\hat{\mathit{x}}}_{\mathit{k}-1|\mathit{N}}^{\mathit{T}}$$

(23)

$${\mathit{P}}_{\mathit{k},\mathit{k}-1|\mathit{N}}={\mathit{P}}_{\mathit{k}|\mathit{k}}{\mathit{L}}_{\mathit{k}-1}^{\mathit{T}}+{\mathit{L}}_{\mathit{k}}({\mathit{P}}_{\mathit{k}+1,\mathit{k}|\mathit{N}}-{\mathit{P}}_{\mathit{k}|\mathit{k}}){\mathit{L}}_{\mathit{k}-1}^{\mathit{T}}$$

(24)

The above formula is initialized by Formula (25).

$${\mathit{P}}_{\mathit{N},\mathit{N}-1|\mathit{N}}=(1-{\mathit{K}}_{\mathit{N}}){\mathit{P}}_{\mathit{N}-1|\mathit{N}-1}$$

(25)

Step M: Figure out the parameter set $\theta$ that can maximize the probability of the given observation signal $Y_N=\{y_1,y_2,\dots ,y_N\}$ occurring by maximizing the likelihood function O. The updated parameters in the i th iteration can be acquired, as shown in Formulas (26) to (30).

$${\mathit{A}}^{(\mathit{i})}=\left(\sum _{\mathit{k}=2}^{\mathit{N}}{\mathit{S}}_{1,0,\mathit{k}}\right){\left(\sum _{\mathit{k}=2}^{\mathit{N}}{\mathit{S}}_{0,0,\mathit{k}}\right)}^{-1}$$

(26)

$${\mathit{Q}}^{(\mathit{i})}={\displaystyle \frac{1}{\mathit{N}-1}\left(\sum _{\mathit{k}=1}^{\mathit{N}}{\mathit{S}}_{1,1,\mathit{k}}-{\mathit{A}}^{(\mathit{i})}{\left(\sum _{\mathit{k}=1}^{\mathit{N}}{\mathit{S}}_{1,0,\mathit{k}}\right)}^{\mathit{T}}\right)}$$

(27)

$${\mathit{H}}^{(\mathit{i})}=\left(\sum _{\mathit{k}=1}^{\mathit{N}}{\mathit{y}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}^{\mathit{T}}\right){\left(\sum _{\mathit{k}=1}^{\mathit{N}}{\mathit{S}}_{1,1,\mathit{k}}\right)}^{-1}$$

(28)

$${\mathit{R}}^{(\mathit{i})}={\displaystyle \frac{1}{\mathit{N}}\sum _{\mathit{k}=1}^{\mathrm{\infty }}({\mathit{y}}_{\mathit{k}}-{\mathit{H}}^{(\mathit{i})}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}+{\mathit{y}}_{\mathit{k}}-{\mathit{H}}^{(\mathit{i})}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{N}}{)}^{\mathit{T}}-{\mathit{H}}^{(\mathit{i})}{\mathit{P}}_{\mathit{k}|\mathit{N}}({\mathit{H}}^{(\mathit{i})}{)}^{\mathit{T}}}$$

(29)

$${\mathit{P}}_0^{(\mathit{i})}={\mathit{P}}_{1|\mathit{N}}{\mathit{m}}_0^{(\mathit{i})}={\hat{\mathit{x}}}_{1|\mathit{N}}$$

(30)

Suppressing the colored measurement noise.

Due to the colored nature of measurement noise, the filtering effect of the standard KF is always difficult to achieve the best. Based on the “measurement time difference” between the measurement value of adjacent moments, a pseudo measurement variable $z_k$ is constructed based on the standard KF, as shown in Formula (31), cleverly reducing the influence of the colored measurement noise on the filtering effect. The core of this strategy lies in introducing a pseudo-measurement variable, which fully considers the temporal correlation of measurement noise, thereby enabling more accurate tracking of its dynamic changes.^31,40 As shown in Formulas (32) to (36).

$${\mathit{v}}_{\mathit{k}}={\mathit{\psi }}_{\mathit{k}-1}{\mathit{v}}_{\mathit{k}-1}+{\mathit{\epsilon }}_{\mathit{k}-1}{\mathit{z}}_{\mathit{k}}={\mathit{y}}_{\mathit{k}+1}-{\mathit{\psi }}_{\mathit{k}}{\mathit{y}}_{\mathit{k}}={\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\mathit{x}}_{\mathit{k}}+{\mathit{H}}_{\mathit{k}}{\mathit{w}}_{\mathit{k}}+{\mathit{\epsilon }}_{\mathit{k}}$$

(31)

$${\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}={\mathit{H}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}-{\mathit{\psi }}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}$$

(32)

$${\mathit{S}}_{\mathit{k}}={\mathit{Q}}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}$$

(33)

$${\mathit{R}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}={\mathit{H}}_{\mathit{k}}{\mathit{Q}}_{\mathit{k}}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k}}$$

(34)

$${\mathit{Q}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}={\mathit{Q}}_{\mathit{k}}+{\mathit{S}}_{\mathit{k}}{\mathit{R}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}^{-1}{\mathit{S}}_{\mathit{k}}^{\mathit{T}}$$

(35)

$${\mathit{A}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}={\mathit{A}}_{\mathit{k}}-{\mathit{J}}_{\mathit{k}}{\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\mathit{J}}_{\mathit{k}}={\mathit{S}}_{\mathit{k}}{\mathit{R}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}^{-1}$$

(36)

where ${\psi }_k$ as the measurement noise transfer matrix can describe the noise type existing during the whole measurement. When ${\psi }_k$ is small, it indicates that the measurement noise is approximately stationary white Gaussian noise. Conversely, when ${\psi }_k$ is large, it indicates that the stability of measurement noise is poor. ${\epsilon }_k$ represents white Gaussian noise. Moreover, $H_{k,new}$ , $R_{k,new}$ , $Q_{k,new}$ , and $A_{k,new}$ represent the new measurement matrix, the new measurement noise covariance, the new process noise covariance, and the new state transition matrix, respectively.

Accordingly, the prediction and update formula of the standard KF also need to be modified, as shown in Formulas (37) to (43).

(i) Prediction step:

$${\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}={\mathit{A}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\hat{\mathit{x}}}_{\mathit{k}-1|\mathit{k}-1}+{\mathit{J}}_{\mathit{k}}({\mathit{z}}_{\mathit{k}-1}-{\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\hat{\mathit{x}}}_{\mathit{k}-1|\mathit{k}-2})$$

(37)

$${\mathit{P}}_{\mathit{k}|\mathit{k}-1}={\mathit{A}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\mathit{P}}_{\mathit{k}-1|\mathit{k}-1}{\mathit{A}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}^{\mathit{T}}+{\mathit{Q}}_{\mathit{k}-1,\mathit{n}\mathit{e}\mathit{w}}$$

(38)

(ii) Update step:

$${\mathit{d}}_{\mathit{k}}={\mathit{z}}_{\mathit{k}}-{\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}$$

(39)

$${\mathit{S}}_{\mathit{k}}={\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}{\mathit{P}}_{\mathit{k}|\mathit{k}-1}{\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}$$

(40)

$${\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}|\mathit{k}-1}{\mathit{H}}_{\mathit{k},\mathit{n}\mathit{e}\mathit{w}}^{\mathit{T}}{\mathit{S}}_{\mathit{k}}^{-1}$$

(41)

$${\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}}={\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}{\mathit{d}}_{\mathit{k}}$$

(42)

$${\mathit{P}}_{\mathit{k}|\mathit{k}}={\mathit{P}}_{\mathit{k}|\mathit{k}-1}-{\mathit{K}}_{\mathit{k}}{\mathit{S}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$

(43)

Iteration stop criterion of EM

The essence of the first KF algorithm is to optimize the unknown core parameter set $\theta$ of KF using the EM framework. By iterating E-step and M-step, the true noise statistical characteristics of the original signal can be gradually approached, thus making KF maintain good estimation performance in unknown noise environments. In order to balance the filtering effect and running efficiency of the algorithm, the cosine similarity between the original signal and the filtered signal is calculated as the termination condition for iteration. The closer its value is to 1, the higher the morphological similarity between the filtered signal and the original signal. The calculation method for cosine similarity is shown in Formula (44), where $AEC{G}_S$ is the original AECG and $AEC{G}_f$ is the filtered AECG by the KF based on the EM.

$$\mathrm{Cosine}\mathrm{Similarity}={\displaystyle \frac{AEC{G}_S\cdot AEC{G}_f}{AEC{G}_{S}AEC{G}_f}}$$

(44)

As shown in Figure 3(a), the black curve shows $AEC{G}_f$ . In contrast, the blue curve represents $AEC{G}_S$ . Obviously, $AEC{G}_f$ has a lower noise level compare to $AEC{G}_S$ . As shown in Figure 3(b), during the iteration process of the EM algorithm, the cosine similarity between $AEC{G}_f$ and $AEC{G}_S$ gradually increases with the growth of iteration steps at the beginning. As the number of iteration steps continues to increase, its increasing trend slows down (Figure 3(c)). Then, it gradually approaches a value close to 1, which also indicates the stability and convergence of the EM algorithm. In order to improve the running efficiency of the algorithm, the stopping condition for the iteration of the EM algorithm is set to be that the increment of the cosine similarity is less than 0.05.

Figure 3.

The original AECG (blue) processed only by band-pass filtering and the AECG (black) processed by KF based on the em algorithm (a); the change curve of cosine similarity between $\mathit{A}\mathit{E}\mathit{C}{\mathit{G}}_{\mathit{f}}$ and $\mathit{A}\mathit{E}\mathit{C}{\mathit{G}}_{\mathit{S}}$ as the number of iterations increases (b); the increment of the cosine similarity between $\mathit{A}\mathit{E}\mathit{C}{\mathit{G}}_{\mathit{f}}$ and $\mathit{A}\mathit{E}\mathit{C}{\mathit{G}}_{\mathit{S}}$ as the number of iterations increases (c).

MECG extraction

After the original AECG is processed by the first KF based on the EM algorithm, most of the various external noise in the original signal is significantly reduced. The mixed residual signal mainly includes FECG, MECG, and a small amount of residual noise. To ensure the morphological integrity of the extracted MECG, the MECG template is fitted by the Gaussian kernel function and is employed in the construction of the state equation of KF. Additionally, so as to improve the estimation accuracy and robustness of KF, the forgetting factor ${\alpha }_k$ is introduced to adaptively update the process noise covariance matrix Q and the measurement noise covariance matrix R based on the innovation $d_k$ and residual vector $e_k$ , separately. The state vector of MECG is mainly composed of phase and amplitude. For MECG of each heartbeat, the linear interpolation in the [- π, π] interval is used to allocate its phase and the amplitude is described by the template equation composed of Gaussian kernel functions.

Locate R peak of MECG

In order to ensure the accuracy of the R peak of MECG, the QRS detector based on the Pan-Tompkins algorithm ⁴¹ is used to identify the R peak of MECG from the mixed residual signal. Firstly, the mixed signal is differentiated to enhance the peak characteristics of the QRS complex. Subsequently, the moving integral operation with a certain window width is performed on the differentiated signal, which smoothes the waveform and highlights the QRS complex. Next, the fixed amplitude threshold and absolute refractory period parameter are set to locate the R peak of MECG, as shown in Figure 4.

Figure 4.

The R peak of MECG in the mixed residual signal.

MECG template construction

In order to reduce the impact of the residual MECG on subsequent FECG extraction, nine Gaussian kernel functions are selected to construct the MECG template,^36,42 the two state equations of which can be described by Formulas (45) to (46). Besides, the reason and details can be found at Figure S1 in Supporting Information (SP).

$${\mathit{z}}_{\mathit{k}+1}={\mathit{z}}_{\mathit{k}}-\sum _{\mathit{i}=1}^{\mathit{N}}{\displaystyle \frac{{\mathit{\alpha }}_{\mathit{i}}{\mathit{\omega }}_{\mathit{t}}\mathit{d}\mathit{t}}{{({\mathit{b}}_{\mathit{i}})}^2}({\mathit{\phi }}_{\mathit{k}}-{\mathit{\theta }}_{\mathit{i}}){\mathit{e}}^{-\frac{{({\mathit{\phi }}_{\mathit{k}}-{\mathit{\theta }}_{\mathit{i}})}^2}{2{({\mathit{b}}_{\mathit{i}})}^2}}}$$

(45)

$${\mathit{\phi }}_{\mathit{k}+1}={\mathit{\phi }}_{\mathit{k}}+{\mathit{\omega }}_{\mathit{k}}\mathit{d}\mathit{k}$$

(46)

In the above formula, $z_k$ and ${\phi }_k$ individually represent the amplitude and the phase of the original signal. ${\alpha }_i$ , $b_i$ , and ${\theta }_i$ correspond to the amplitude, width, and center parameter of the $i$ th Gaussian kernel function, respectively. ${\omega }_k$ is the angular velocity, where the coefficient $w=[{\mu }^T,{\omega }_k{]}^T=[{\alpha }_1,b_1,{\theta }_1,\dots ,{\alpha }_9,b_9,{\theta }_9,{\omega }_k{]}^T$ . $dk$ represents the minimum sampling time interval.

T and N represent the duration of a single heartbeat and the number of discrete points in a single heartbeat, respectively. ${\omega }_k$ can be calculated by Formula (47).

$${\mathit{\omega }}_{\mathit{k}}={\displaystyle \frac{2\mathit{\pi }}{\mathit{T}}}$$

(47)

According to the definitions of the two state equations, the initial value, upper and lower thresholds of all the elements in the Gaussian kernel function coefficient $\mu$ are determined, as shown in Formulas (48) to (53).

$${\mathit{\alpha }}_{\mathit{i},0}={\mathit{A}}_{\mathit{i}}$$

(48)

$${\mathit{b}}_{\mathit{i},0}={\displaystyle \frac{2\mathit{\pi }{\mathit{W}}_{\mathit{i}}}{5\mathit{N}}}$$

(49)

$${\mathit{\theta }}_{\mathit{i},0}={\displaystyle \frac{2\mathit{\pi }{\mathit{L}}_{\mathit{i}}}{\mathit{N}}}$$

(50)

$${\mathit{\alpha }}_{\mathit{i}}\in [{\mathit{A}}_{\mathit{m}\mathit{i}\mathit{n}}{\mathit{\alpha }}_{\mathit{i},0},{\mathit{A}}_{\mathit{m}\mathit{a}\mathit{x}}{\mathit{\alpha }}_{\mathit{i},0}]$$

(51)

$${\mathit{b}}_{\mathit{i}}\in [0.01,{\mathit{b}}_{\mathit{i},0}+1]$$

(52)

$${\mathit{\theta }}_{\mathit{i}}\in \left[{\mathit{\theta }}_{\mathit{i},0}-{\displaystyle \frac{\mathit{\pi }}{5},{\mathit{\theta }}_{\mathit{i},0}+{\displaystyle \frac{\mathit{\pi }}{5}}}\right]$$

(53)

In the above formula, $A_i$ , $W_i$ , and $L_i$ represent the amplitude, width, and center position of the i th Gaussian kernel function, respectively.

The optimal value of each Gaussian kernel function coefficient $\mu$ is determined based on the nonlinear least square method, where the upper and lower thresholds of $\mu$ are as shown in Formulas (54) to (57).

$$\underset{\mathit{\mu }}{min}\mathit{z}(\mathit{\mu },{\mathit{\omega }}_{\mathit{k}}{)}^2$$

(54)

$$\mathbf{s}.\mathbf{t}.{\mathit{\mu }}_{\mathit{L}}\le \mathit{\mu }\le {\mathit{\mu }}_{\mathit{H}}$$

(55)

$$\begin{array}{rl}{\mathit{\mu }}_{\mathit{L}}=& [{\mathit{A}}_{\mathit{m}\mathit{i}\mathit{n}}{\mathit{\alpha }}_{1,0},0.01,{\mathit{\theta }}_{1,0}-{\displaystyle \frac{\mathit{\pi }}{5},}\\ & \dots ,{\mathit{A}}_{\mathit{m}\mathit{i}\mathit{n}}{\mathit{\alpha }}_{9,0},0.01,{\mathit{\theta }}_{9,0}-{\displaystyle \frac{\mathit{\pi }}{5}}]\end{array}$$

(56)

$$\begin{array}{rl}{\mathit{\mu }}_{\mathit{H}}=& [{\mathit{A}}_{\mathit{m}\mathit{a}\mathit{x}}{\mathit{\alpha }}_{1,0},{\mathit{b}}_{1,0}+1,{\mathit{\theta }}_{1,0}+{\displaystyle \frac{\mathit{\pi }}{5},}\\ & \dots ,{\mathit{A}}_{\mathit{m}\mathit{a}\mathit{x}}{\mathit{\alpha }}_{9,0},{\mathit{b}}_{9,0}+1,{\mathit{\theta }}_{9,0}+{\displaystyle \frac{\mathit{\pi }}{5}}]\end{array}$$

(57)

After the above steps, the MECG template is constructed based on nine Gaussian kernel functions using the first 30 cycles of the above mixed residual signal, as shown in Figure 5.

Figure 5.

The MECG template fitted by nine Gaussian kernel functions.

The dynamic estimation of MECG based on the second KF

As shown in Formulas (58) and (59), the state equation and measurement equation comprise the foundation of the KF. The state vector $x_k=[{\phi }_k,z_k{]}^T$ consists of the phase ${\phi }_k$ and amplitude $z_k$ of the MECG model. while the phase ${\phi }_k$ and amplitude $s_k$ of the collected MECG form the measurement vector $y_k=[{\phi }_k,s_k{]}^T$ .

$$x_k=F(x_{k-1},u_k)$$

(58)

$$y_k=H{x}_k+v_k=x_k+v_k$$

(59)

where H is the measurement matrix, representing the relationship between the state and measurement vectors. $u_k\sim N(0,Q)$ and $v_k\sim N(0,R)$ are the process noise and measurement noise of the KF. The subscript k represents the sampling time of the MECG.

Different from the strategy of the fixed parameters of the first KF based on the EM algorithm, R and Q is adaptively updated in each iterative step of the second KF by incorporating the forgetting factor ${\alpha }_k$ based on the statistical properties of the innovation $d_k$ and the residual vector $e_k$ . ³³ The initial value of the parameters of the second KF is shown in Formulas (60) to (64).

$${\mathit{x}}_0=\left[\begin{array}{c}-\mathit{\pi }\\ 0\end{array}\right]$$

(60)

$${\mathit{P}}_0=\left[\begin{array}{c}{(2\mathit{\pi })}^2\\ 10{|{\mathit{s}}_{\mathit{m}\mathit{a}\mathit{x}}|}^2\end{array}\right]$$

(61)

$${\mathit{R}}_0=\left[\begin{array}{c}{10}^{-5}{\left({\displaystyle \frac{\mathit{\omega }}{{\mathit{f}}_{\mathit{s}}}}\right)}^2\\ 10{\mathit{\sigma }}_{\mathit{s}}^2\end{array}\right]$$

(62)

$${\mathit{Q}}_0=\mathit{d}\mathit{i}\mathit{a}\mathit{g}\left[\begin{array}{ccc}{\mathit{p}}_{\mathit{\alpha }}^2& {\mathit{p}}_{\mathit{b}}^2& {\mathit{p}}_{\mathit{\theta }}^2\end{array}\begin{array}{cc}{\mathit{\sigma }}_{\mathit{\phi }}^2& {\mathit{\sigma }}_{\mathit{s}}^2\end{array}\right]$$

(63)

$$\mathbf{H}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$

(64)

where $s_{max}$ is the maximum value of MECG amplitude, $\omega$ is the angular frequency, and $f_s$ is the sampling frequency. ${\sigma }_s$ and ${\sigma }_{\phi }$ are the standard deviations of amplitude and angular frequency, respectively. In addition, $p_{\alpha }$ , $p_b$ , and $p_{\theta }$ represent the row vectors of all Gaussian kernel function coefficients ${\alpha }_i$ , $b_i$ , ${\theta }_i$ .

Formula (58) reveals that the state equation possesses nonlinear characteristics, necessitating linearization prior to its application in subsequent processing stages, as shown in Formula (65). ${\hat{x}}_{k|k-1}$ and ${\hat{x}}_{k-1|k-1}$ represent the prior and posterior estimation of x at timestamp k and $k-1$ , respectively. It can be seen that $A_k$ and $B_k$ can be derived by computing their respective Jacobian matrices (Formulas (66) to (67)), where N is the number of Gaussian kernel functions.

$$x_k=F(x_{k-1},u_k)\approx {\hat{x}}_{k|k-1}+A_k(x_{k-1}-{\hat{x}}_{k-1|k-1})+B_k{u}_k$$

(65)

$$A_k=\left[\begin{array}{cc}1& 0\\ {\displaystyle \frac{\mathrm{\partial }{z}_k}{\mathrm{\partial }{\phi }_{k-1}}}& 1\end{array}\right]$$

(66)

$$B_k=\left[\begin{array}{cccc}{0}_{1\times N}& 0_{1\times N}& 0_{1\times N}& 0_{1\times N}\\ {\displaystyle \frac{\mathrm{\partial }{z}_k}{\mathrm{\partial }{\alpha }_i}{|}_{1\times N}}& {\displaystyle \frac{\mathrm{\partial }{z}_k}{\mathrm{\partial }{b}_i}{|}_{1\times N}}& {\displaystyle \frac{\mathrm{\partial }{z}_k}{\mathrm{\partial }{\theta }_i}{|}_{1\times N}}& {\displaystyle \frac{\mathrm{\partial }{z}_k}{\mathrm{\partial }{\omega }_{k-1}}}\end{array}\right]$$

(67)

The amplitude of MECG Gaussian model exhibits a high degree of nonlinearity. To facilitate the construction of the second KF, the expression of MECG template is linearized in this paper through the application of the Taylor expansion formula. According to Formula (45), the state vector and its error covariance matrix are as depicted in Formulas (68) to (69). $d{\theta }_i={\phi }_{k-1}-{\theta }_i$ represents the phase difference between the phase of the state vector and the center position of the Gaussian kernel function and

$${\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}=\left[\begin{array}{c}{\hat{\mathit{\phi }}}_{\mathit{k}-1}+\mathit{\omega }\mathit{d}\mathit{k}\\ {\hat{\mathit{z}}}_{\mathit{k}-1|\mathit{k}-1}-\mathit{d}\mathit{t}\sum _{\mathit{i}=1}^{\mathit{N}}\left({\displaystyle \frac{\mathit{\omega }{\mathit{\alpha }}_{\mathit{i}}}{{\mathit{b}}_{\mathit{i}}^2}}\right)\mathit{d}{\mathit{\theta }}_{\mathit{i}}{\mathit{e}}^{-\frac{{(\mathit{d}{\mathit{\theta }}_{\mathit{i}})}^2}{2{({\mathit{b}}_{\mathit{i}})}^2}}\end{array}\right]$$

(68)

$${\mathit{P}}_{\mathit{k}|\mathit{k}-1}={\mathit{A}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}-1|\mathit{k}-1}{\mathit{A}}_{\mathit{k}}^{\mathit{T}}+{\mathit{Q}}_{\mathit{k}-1}$$

(69)

$R$ and Q respectively describe the uncertainty of the measurement and the system itself, which jointly ensure the accuracy and robustness of the second KF. ⁴³ In this article, based on the residual vector $e_k$ and the innovation $d_k$ , the forgetting factor ${\alpha }_k$ is introduced to update R and Q iteratively.

The adaptive updating of the measurement noise covariance matrix R.

In the framework of KF, the positive definiteness of R is the key to ensure the stability and convergence of KF. Because the performance parameters of the sensor itself are usually affected by temperature, humidity, radiation and other factors, the measurement noise characteristics of which often fluctuate with the change of time and environment. Therefore, the fixed R is difficult to adapt to the above dynamic change. On the other hand, the adaptive method can adjust R according to the real-time signal to ensure that the current noise level is actually reflected. The residual vector $e_k$ is defined as the difference between the actual measurement value and the posterior estimate value of the state vector, which is an important indicator to describe the adaptability of the KF model. By monitoring the statistical characteristics of $e_k$ , the elements of R are adaptively adjusted to always maintain positive definiteness of R during the update process, avoiding the degradation of filter performance. The specifics of the above adjustment are outlined in Formulas (70) to (71).

$${\mathit{e}}_{\mathit{k}}={\mathit{y}}_{\mathit{k}}-{\mathit{H}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}}$$

(70)

$${\mathit{R}}_{\mathit{k}}=\mathit{E}[{\mathit{e}}_{\mathit{k}}{\mathit{e}}_{\mathit{k}}^{\mathit{T}}]+{\mathit{H}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}|\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}$$

(71)

With the passage of time, the characteristic change of the measurement noise accumulated by the sensor for a long time may lead to the saturation of the estimated value of $R_k$ , thus losing the accurate description of the current noise characteristics. To eliminate the influence of $R_k$ transient fluctuation on its stability, the strategy of low-pass filtering is used to smooth and update $R_k$ , as shown in Formula (72). For non-stationary signals, newly collected signal by sensor often contains more useful information that can reflect the trend of signal changes compared to the old signal. The forgetting factor ${\alpha }_k$ can balance the weight of old and new signal in the update process. Moreover, the initial value of ${\alpha }_k$ is determined by empirical analysis. A large ${\alpha }_k$ prioritizes the new signal, sacrificing KF stability for faster convergence. Conversely, a smaller ${\alpha }_k$ emphasizes the old signal, enhancing stability while slowing convergence.

$${\mathit{R}}_{\mathit{k}}=(1-{\mathit{\alpha }}_{\mathit{k}}){\mathit{R}}_{\mathit{k}-1}+{\mathit{\alpha }}_{\mathit{k}}({\mathit{e}}_{\mathit{k}}{\mathit{e}}_{\mathit{k}}^{\mathit{T}}+{\mathit{H}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}|\mathit{k}-1}{\mathit{H}}_{\mathit{k}}^{\mathit{T}}){\mathit{\alpha }}_{\mathit{k}}=\mathit{\lambda }/(1-{\mathit{\lambda }}^{\mathit{k}+1})$$

(72)

Therefore, the coefficient $\lambda$ is introduced to balance the stability and convergence speed of the second KF by iteratively optimizing the forgetting factor ${\alpha }_k$ according to Formula (69). In this paper, the coefficient $\lambda$ is set to 0.95, as determined through empirical analysis in the SP (Figure S2–3 and Table S1). With the increasing of iteration rounds, the change process of ${\alpha }_k$ is shown in Figure 6. It can be seen that at the initial stage, ${\alpha }_k$ is a relatively large value, which accelerates the convergence speed of KF. Subsequently, ${\alpha }_k$ decreases exponentially and gradually reaches the rapid stability stage, so as to ensure the stability of KF.

The adaptive updating of the process noise covariance matrix Q.

Figure 6.

The update curve of ${\mathit{\alpha }}_{\mathit{k}}$ with the increase of iteration rounds.

For the classical KF, the Q matrix reflects the model uncertainty caused by process noise in the state estimation stage. Inaccurate Q estimation may lead to an error in filter gain calculation, which will affect the effect of the filter. Due to the change of external environment and the position of the fetus in the abdomen in the measurement process, the characteristics of the AECG process noise often change dynamically with time. The innovation $d_k$ is the difference between the actual measured value and the prior estimated value, which reflects the prediction error. If Q can be adjusted based on the changing characteristics of process noise, the filter can quickly adapt to the change of process noise and adjust the effect of state estimation in time.

The strategy of low-pass filtering is utilized in this paper, similar to the adaptive update strategy of R. To ensure the stability of the second KF, Q is adaptively updated based on the innovation $d_k$ , the residual vector $e_k$ and the forgetting factor ${\alpha }_k$ . The details of the adaptive update strategy of Q are provided in Formulas (73) to (74).

$${\mathit{d}}_{\mathit{k}}={\mathit{y}}_{\mathit{k}}-{\mathit{H}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}|\mathit{k}-1}$$

(73)

$${\mathit{Q}}_{\mathit{k}}=(1-{\mathit{\alpha }}_{\mathit{k}}){\mathit{Q}}_{\mathit{k}-1}+{\mathit{\alpha }}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}{\mathit{d}}_{\mathit{k}}{\mathit{d}}_{\mathit{k}}^{\mathit{T}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$

(74)

As shown in Figure 7, MECG in the preprocessed original AECG can be accurately extracted by the second adaptive KF with little MECG component left, which lays the foundation for the next step of high-quality extraction of FECG.

Figure 7.

The original AECG processed by EM (a); the extracted MECG from AECG by the first KF based on EM and the second KF serially (b); the residual signal with little MECG component left (c).

In KF, the diagonal elements of the error covariance matrix P correspond to the variance of the estimation error of each state variable, while the non-diagonal elements reflect the error correlation between state variables. As shown in Formula (75), the decrease of the Frobenius norm of P means that the estimation error variance of the state vector shrinks overall and the estimation accuracy synchronously improves. As P gradually decreases and stabilizes, the error compensation of KF prediction and update reaches dynamic balance, which is also the core indicator of filter convergence. The dynamic variation curve of the norm of the P of the second KF during the filtering process is shown in Figure 8.

$$P_F=\sqrt{\sum _{i=1}^m\sum _{j=1}^n{|p_{ij}|}^2}$$

(75)

Figure 8.

The dynamic variation curve of the $\mathit{P}$ norm of the second KF during the filtering process.

FECG extraction

After the processing of the first and second KF, most of the noise and MECG component in the original AECG have been basically filtered out except for FECG and a little MECG residual component. In order to further improve the SNR of the target FECG, seven Gaussian kernel functions are selected to construct the FECG template, which can be seen in Figure S1 of SP for details. Then, a similar strategy in MECG extraction is also employed to adaptively update R and Q matrix, so as to restore the time domain characteristics of Q, R, S, P, T waves of FECG to the greatest extent and ensure the morphological integrity of the extracted FECG.

Since the time interval between two adjacent R peaks of FECG is generally shorter than that of MECG, and the amplitude of its R peak is also relatively low. Therefore, it is necessary to shorten the absolute refractory period and reduce the amplitude threshold of the QRS detector based on the Pan-Tompkins method, thereby ensuring the accuracy of R peak detection. The detection result is shown in Figure 9. In addition, based on the innovation $d_k$ , the residual vector $e_k$ and the forgetting factor ${\alpha }_k$ , the third KF is also performed on the residual signal which mainly contains FECG and little residual MECG. Figure 9(b) shows the template of FECG. As shown in Figure 9(c), the black curve is the pure FECG signal and the blue curve is the extracted FECG, indicating that the similarity between the extracted FECG and the pure FECG template is high and the morphological integrity is also rather strong.

Figure 9.

The R peak position of FECG (a); the FECG template fitted with seven Gaussian kernel functions (b); the extracted FECG (blue) and the pure FECG (black) (c).

Similar to Figure 8, Figure 10 represents the dynamic variation curve of the norm of P of the third KF during the filtering process. It can be seen that the above curve initially drops rapidly, then slows down and finally approaches a certain small value, indicating the third KF has reached a convergence state.

Figure 10.

The dynamic variation curve of the norm of P of the third KF during the filtering process.

Results

In order to quantitatively compare the extraction effect of the proposed algorithm and three widely used algorithms (PCA, EKF and TS) in extracting FECG from AECG, 25 groups of AECG generated by the FECGSYN toolbox, 10 groups of clinical AECG collected from ADFECGDB database, 10 groups of clinical AECG sourced from the FECGDARHA database and three groups of real AECG collected from the AECG recording system designed are used to comprehensively evaluate the extraction accuracy, SNR and morphological integrity of the extracted FECG. The above four algorithms are implemented in MATLAB 2024B environment.

Firstly, the 25 groups of AECG generated by FECGSYN toolbox are preprocessed by the first KF based on the EM algorithm, which aims to suppress various noise components to the greatest extent and retain the complete, high SNR, residual signal mixed by FECG and MECG. Here, a group of AECG with the SNR of 6 dB is utilized as an example, its preprocessing result is shown in Figure 11. The blue curve represents the original AECG, while the black curve represents the preprocessed signal by EM. It can be seen that the noise components, including power frequency interference, motion artifacts, EMG interference, baseline drift and so on, have been effectively filtered out.

Figure 11.

The comparison of a group of 6 dB AECG before (blue) and after (black) preprocessing.

Then, the four algorithms mentioned above are respectively applied to 25 groups of AECG after preprocessing to extract the target MECG. Similarly, taking the 6 dB AECG processed by EM algorithm as an example. The extracted MECG are shown in Figure 12. Figure 12(a) shows the pure MECG, Figure 12(b) shows the MECG extracted by the algorithm proposed and Figure 12(c) to (e) shows the MECG extracted by PCA, EKF and TS algorithms, respectively. As depicted in Figure 12, when the QRS complexes of FECG and MECG overlap, the proposed algorithm effectively extracts R peak morphologies from the MECG that closely match the pure MECG template. In contrast, the R peak morphologies extracted by the other three algorithms show considerable fluctuations. Due to the superfluous residual MECG or excessive extraction to FECG, the time domain features of P, T, and U waves are poor in periodicity, which makes it difficult to ensure the noise level and morphological integrity of the remaining FECG.

Figure 12.

The original MECG template (a); MECG signals obtained through four different algorithms (the proposed algorithm, PCA, EKF and TS) from 6 dB AECG processed by EM algorithm (b–e).

Finally, the four different algorithms are again applied to 25 groups of residual signals after noise removal and MECG extraction. These signals mainly include FECG, the residual MECG component and a small amount of noise. Equally, taking the above AECG with an SNR of 6 dB as an example, the extraction results are shown in Figure 13. Figure 13(a) shows the original pure FECG as a reference signal. Figure 13(b) shows the FECG extracted from the residual signal by the algorithm proposed in this article. Figure 13(c) to (e) shows the FECG extracted by PCA, EKF and TS algorithms, respectively. Through detailed comparison, it can be found that the FECG extracted by the other three algorithms have the phenomenon of irregular waveform characteristics and abnormal amplitude. Conversely, the FECG extracted by the algorithm proposed in this article is very close to the pure FECG for the amplitude and morphology of P, Q, R, S, T, U waves, which proves the accuracy and stability of the algorithm proposed in FECG extraction.

Figure 13.

The original FECG template (a); FECG extraction from 6 dB AECG using em algorithm preprocessing followed by two consecutive applications of our proposed algorithm, PCA, EKF and TS, respectively (b–e).

In order to comprehensively evaluate the performance of the proposed algorithm for extracting FECG from AECG, 25 groups of AECG with different SNR generated by FECGSYN toolbox are quantitatively analyzed and the ACC, F1, SE and PPV are calculated. Table 2 and Figure 14 show the average of the extraction results of the 25 groups of AECG. The detailed statistics of 25 groups of AECG can be seen in Table S2–6 in SP. Figure 14 shows that all four algorithms have a positive correlation with the SNR of the original AECG on the ACC, F1, SE and PPV of the R peak. In this paper, the average of ACC, F1, SE, and PPV of the proposed algorithm is more than 76%. Especially, in the case of the original AECG with high SNR, the four indicators are more than 90%, which is significantly higher than PCA, EKF, and TS. In addition, SNR is used to describe the noise level of the final extracted FECG. As shown in Table 2, the average SNR of FECG extracted by the proposed algorithm remains above 3.6 dB, while the average SNR of FECG extracted by PCA, EKF, and TS are all significantly lower than that of the algorithm proposed in this paper. The morphological integrity of the extracted FECG is evaluated using the RMSE, which is the root mean square error between the original and extracted FECG. As shown in Table 2, the RMSE of FECG obtained by the proposed algorithm is kept within 370 mV, while the RMSE of FECG obtained by PCA, EKF, and TS is more than 440 mV. These above experimental results further show that the proposed algorithm has significant advantages in the morphological integrity of FECG. In conclusion, the algorithm proposed in this paper has good extraction effect for extracting FECG from AECG, which is hope to provide an important reference for doctors to accurately analyze the clinical characteristics of FECG and formulate more reasonable clinical diagnosis scheme in the future.

Figure 14.

The line chart of the extraction effect (SNR, RMSE, ACC, F1, PPV and SE) of 25 groups of AECG with the proposed algorithm and other three algorithms.

Table 2.

The average of the extraction effect of the 25 groups of AECG with the proposed algorithm and the other three algorithms.

Method SNR		ACC	F1	SE	PPV	SNR(dB)	RMSE(mV)
Proposed	3 dB	76.47%	85.07%	81.91%	88.65%	3.6985	363.8333
	5 dB	84.56%	91.34%	88.98%	93.89%	3.8372	367.5977
	6 dB	86.84%	92.79%	90.73%	94.98%	4.1484	359.8935
	7 dB	89.24%	94.21%	92.36%	96.18%	4.1605	364.0441
	10 dB	90.90%	95.10%	93.43%	96.87%	4.1992	350.1645
PCA	3 dB	67.12%	77.36%	73.22%	82.11%	2.3208	555.5079
	5 dB	75.23%	84.07%	80.50%	88.05%	2.8487	515.1599
	6 dB	75.93%	84.73%	81.04%	88.89%	2.9315	502.9774
	7 dB	81.87%	89.69%	86.44%	93.23%	3.2578	483.5033
	10 dB	83.15%	90.21%	86.87%	93.86%	3.1340	476.1205
EKF	3 dB	39.21%	52.90%	47.73%	60.97%	−4.8124	530.6829
	5 dB	49.10%	63.78%	57.14%	73.36%	−3.0123	603.6893
	6 dB	48.07%	61.32%	55.84%	69.61%	−4.4195	559.0838
	7 dB	60.79%	73.06%	67.26%	80.98%	−3.4043	584.5207
	10 dB	55.24%	68.85%	63.12%	76.74%	−4.0478	548.6258
TS	3 dB	65.55%	75.59%	71.61%	80.15%	2.5074	516.7678
	5 dB	72.98%	82.55%	78.76%	86.83%	3.2759	469.8702
	6 dB	73.80%	83.31%	79.43%	87.67%	3.3846	456.9305
	7 dB	80.35%	88.80%	84.83%	93.21%	3.5569	445.6205
	10 dB	79.65%	88.01%	83.78%	92.75%	3.4018	442.6383

Moreover, in order to further verify the effectiveness of the algorithm proposed in this paper for the actual collection AECG from pregnant women, the performance of the above four algorithms are compared in extracting FECG from AECG based on 10 groups of clinical collected AECG of pregnant women in ADFECGDB database. Since the database does not include the noiseless FECG as the standard reference signal, the previously mentioned SNR and RMSE cannot be used to accurately evaluate the SNR and morphological integrity of the extracted FECG. Therefore, the average performance (ACC, F1, SE and PPV) of the four algorithms on the mentioned 10 groups of AECG are mainly evaluated. As shown in Table 3, the average performance of the algorithm proposed in the above four indicators is higher than 90%, which is superior to the other three algorithms. PCA and EKF have the approximate performance, slightly lower than the algorithm proposed. TS has the worst performance, which is 10% lower than the algorithm proposed. To sum up, the above results show that the proposed algorithm is more suitable for extracting FECG from the actual collection AECG compared with the other three algorithms.

Table 3.

The average of the extraction effect of 10 groups of AECG from ADFECGDB database with the proposed algorithm and other three algorithms.

Method	ACC(%)	F1(%)	SE(%)	PPV(%)
Proposed	89.842	94.033	94.429	94.597
PCA	83.931	90.506	88.141	93.081
EKF	87.219	92.382	90.978	93.903
TS	73.146	82.94	78.921	88.031

Besides, to further assess the robustness of the proposed method on late-term clinical recordings, the proposed algorithm with the other three methods (PCA, EKF, and TS) were compared using 10 groups of AECG from the FECGDARHA database. Since standardized expert annotations are available, ACC, F1, SE, and PPV were adopted for evaluation. As reported in Table 4, the proposed method achieves the best overall performance. PCA and TS yield comparable but consistently lower results, whereas EKF shows a clear degradation across all metrics. It is likely because the FECGDARHA exhibits stronger nonstationary interference and larger morphological variability, under which the EKF's fixed modeling assumptions and parameter tuning are less able to track rapid changes, leading to degraded extraction performance. Overall, these results indicate that the proposed algorithm provides more reliable FECG extraction on FECGDARHA than the other three algorithms.

Table 4.

The average of the extraction effect of 10 groups of AECG from FECGDARHA database with the proposed algorithm and other three algorithms.

Method	ACC(%)	F1(%)	SE(%)	PPV(%)
Proposed	85.18	91.46	90.68	92.55
PCA	82.23	89.77	92.55	87.31
EKF	63.03	74.58	73.69	76.98
TS	82.28	89.58	90.98	88.23

To assess whether the proposed method significantly outperforms the baselines, we performed statistical tests on the paired results for ACC, F1, SNR and RMSE. Normality of the paired differences was first checked using the Shapiro–Wilk test. Depending on the outcome, a paired t-test was applied when the differences were normally distributed, whereas the Wilcoxon signed-rank test was used otherwise.

As summarized in Table 5, the proposed method shows highly significant improvements on FECGSYN across all metrics (P < .001). On ADFECGDB, it remains significantly better in ACC and SE (P < .01), while differences in F1 and PPV versus EKF and PCA are not significant (P > .05); the supplementary box plots (Figure S4–6) further indicate a tighter spread, suggesting better robustness. On FECGDARHA, the proposed method is significantly better than EKF and TS in ACC and PPV (P < .05); compared with PCA, most metrics are comparable (ns), while PPV remains significantly higher.

Table 5.

The statistical significance test results of performance metrics on FECGSYN, ADFECGDB and FECGDARHA databases with the proposed algorithm and other three algorithms.

Metric	Comparison	FECGSYN	ADFECGDB	FECGDARHA
ACC	Proposed vs PCA	<0.001 (***)	<0.001 (***)	0.145 (ns)
	Proposed vs EKF	<0.001 (***)	<0.001 (***)	0.001 (**)
	Proposed vs TS	<0.001 (***)	0.002 (**)	0.042 (*)
F1	Proposed vs PCA	<0.001 (***)	0.002 (**)	0.163 (ns)
	Proposed vs EKF	<0.001 (***)	0.084 (ns)	0.002 (**)
	Proposed vs TS	<0.001 (***)	0.002 (**)	0.053 (ns)
SEN	Proposed vs PCA	<0.001 (***)	<0.001 (***)	0.203 (ns)
	Proposed vs EKF	<0.001 (***)	<0.001 (***)	0.002 (**)
	Proposed vs TS	<0.001 (***)	<0.001 (***)	0.871 (ns)
PPV	Proposed vs PCA	<0.001 (***)	0.051 (ns)	0.008 (**)
	Proposed vs EKF	<0.001 (***)	0.261 (ns)	0.003 (**)
	Proposed vs TS	<0.001 (***)	<0.001 (***)	0.002 (**)
SNR and RMSE	Proposed vs PCA	<0.001 (***)	–	–
	Proposed vs EKF	<0.001 (***)	–	–
	Proposed vs TS	<0.001 (***)	–	–

In addition, three pregnant volunteers are recruited, and their AECG are collected using the hardware system designed in this paper. The FHR curve is recorded synchronously using the most common and most popular FHR monitoring device based on the principle of Doppler ultrasonography. Then, the FHR curves are acquired by the same QRS detector based on the FECG extracted by the four algorithms. The results show that the correlation (R²) between the FHR curve calculated by the proposed algorithm and the FHR curve output by Doppler equipment is the strongest, indicating that the proposed algorithm is more suitable for real scenarios. As shown in Table 6, the R² between the FHR curve calculated by FECG based on the algorithm proposed and the FHR curve output by Doppler equipment is significantly better than the other three algorithms with an average of about 92%. The results of the other three algorithms are at a relatively poor level with the average of about 73%, 80%, and 85%, respectively. Thus, it can be concluded that the portable hardware signal acquisition device integrated with the proposed algorithm proposed in this paper has great application potential in the field of portable and low-cost FHR monitoring.

Table 6.

The correlation (R²) between the FHR curve calculated by the above four algorithms and the FHR curve recorded by Doppler equipment from all three groups of AECG of three pregnant volunteers.

R2	Proposed	PCA	EKF	TS
1	0.9278	0.7467	0.7784	0.8684
2	0.9403	0.7620	0.8143	0.8671
3	0.8992	0.6962	0.8303	0.7467

In order to better demonstrate the stability of the proposed algorithm in practical application, Bland–Altman analysis was performed based on the actual validation results of three patients. As shown in Figure 15, the mean bias is approximately 0, and the 95% limits of agreement range from −1.85 to 1.92 bpm. Moreover, the calculated mean absolute error (MAE) is 0.74 bpm. Most data points fall within the limits of agreement, indicating that the system has high consistency with the clinical standard.

Figure 15.

Bland–Altman plot comparing FHR estimated by the proposed method and the Doppler reference across three subjects.

Discussion

At present, PCA, EKF and TS are commonly used to extract FECG from AECG. However, when FECG and MECG in the original AECG coincide in time domain, PCA and TS cannot clearly separate them, which will seriously affect the quality of the extracted FECG. The performance of EKF is closely related to R and Q. If the two matrices are fixed, EKF will not be able to adapt to the time-varying measurement uncertainty from the external environment and model uncertainty from the change of fetal position, thereby affecting the robustness and stability of the filtering result. Existing research on parameter adaptation still faces significant technical bottlenecks; the method proposed by Niknazar ²⁷ is limited by rigid parameter settings, lacking dynamic adaptability to non-stationary environments. Meanwhile, although Panigrahy^28,29 and Suganthy ³⁰ introduced strategies such as DE, PSO, and WLSR, these stochastic optimization algorithms inherently entail high computational burdens and output uncertainty, making it difficult to guarantee the real-time stability of the filtering process. Due to the above problems, the EM algorithm is employed to optimize the initial filtering parameters of the first KF and a new pseudo measurement variable named “measurement time difference” is constructed, whitening the colored noise to Gaussian noise and further suppressing a variety of noise components in the original AECG. Then, based on the innovation $d_k$ and the residual vector $e_k$ , the forgetting factor ${\alpha }_k$ is also introduced to adaptively update the R and Q of the second and third KFs with the strategy of low-pass filtering. After the above two KFs are executed continuously, the relatively pure FECG is extracted from the residual mixed signal processed by EM algorithm.

Furthermore, in order to quantitatively evaluate the performance of the proposed algorithm and the other three algorithms (PCA, EKF and TS) in the extraction accuracy, SNR and morphological integrity of the extracted FECG, the comparative experiments are conducted based on 25 groups of AECG generated by the FECGSYN toolbox, 10 groups of clinical collected AECG of pregnant women in ADFECGDB database and 10 groups of clinical AECG sourced from the FECGDARHA database. Experimental results show that the proposed algorithm outperforms the other three algorithms in terms of extraction accuracy, SNR and morphological integrity. Further, in order to verify the feasibility of the algorithm proposed in the real FHR monitoring scene, a portable, low-cost, home-based AECG recording system for pregnant women is designed in this paper, which has the advantages of small size, low power consumption and convenient use. Based on this system, AECG of pregnant women can be collected continuously and dynamically. Combined with the algorithm proposed, the system can also extract relatively pure FECG from AECG of primordial pregnant women. In addition, three pregnant women volunteers are also recruited and the hardware system designed in this paper is used to collect the three pregnant women's AECG. At the same time of collecting AECG from pregnant women, the most popular FHR monitoring equipment based on the Doppler ultrasonic principle is used to record the FHR curve. Then, the accuracy of four algorithms for FHR is evaluated based on the collected AECG. The results show that the FHR curve calculated by the proposed algorithm has the strongest correlation with the FHR curve output by Doppler equipment. The above results show that the proposed algorithm and the designed AECG acquisition hardware for pregnant women are hopeful to a novel replacement of the FHR monitoring scheme based Doppler in the near future. However, the above real-world validation of the custom hardware system has a limitation. Due to the pilot nature of this study and practical constraints in subject recruitment and synchronous reference acquisition, only three subjects were included, which may not fully reflect variability across gestational ages, fetal positions, and maternal BMI. Future work will involve a larger clinical study with more diverse participants to further validate the reliability of the method and support its use in portable, low-cost FHR monitoring.

In order to evaluate the time efficiency of the proposed algorithm, the experiment was executed based on the FECGSYN, ADFECGDB, FECGDARHA datasets and the real AECG obtained from the self-designed signal acquisition system in MATLAB 2024b on the PC environment (Windows 10, R7–9750 CPU). The detailed performance is shown in Table 7. The proposed algorithm was settled to extract the MECG and FECG from the AECG, with an average time consumption about 150 ms per heartbeat. The time consumption is significantly shorter than the typical heartbeat duration, highlighting its enormous potential in FECG real-time extraction.

Table 7.

The time consumption of the proposed algorithm.

Datasets	Time consumption
FECGSYN	143.246 ms/heartbeat
ADFECGDB	149.477 ms/heartbeat
FECGDARHA	150.356 ms/heartbeat
The AECG obtained by the designed hardware	148.531 ms/heartbeat

Portable and low-cost FHR monitoring equipment can make it easy and convenient for pregnant women to monitor fetal health status, whether in hospitals, clinics or family environments, which plays an important role on improving the accessibility and convenience of medical services. In particular, regular monitoring at home can significantly reduce the time and economic burden of pregnant women who frequently go to hospitals due to fetal monitoring and help to detect potential health problems of pregnant women and fetuses in the early stage. Besides, compared with the existing expensive FHR monitoring equipment, the low-cost system designed in this paper can reduce the medical expenditure and make the limited medical resources more reasonable distribution, especially for families and regions with poor economic conditions. Therefore, it will be easier to widely promote so that more families can afford high-quality fetal health monitoring services, improving the overall maternal and fetal health level of the society.

Conclusion

In this paper, a novel FECG extraction method based on three adaptive KFs is proposed, which can not only optimize the extraction accuracy of FECG, but also improve the morphological integrity and the SNR of FECG. Compared with the traditional KF, the first KF proposed in this study integrates the EM algorithm for parameter optimization and introduces a novel pseudo-measurement variable based on “measurement time difference.” These innovations effectively whiten colored noise into Gaussian noise and reduce the overall noise level of the original AECG. In order to improve the performance of the second and third KFs in extracting FECG from AECG processed by the first KF, the residual vector $e_k$ and the innovation $d_k$ are combined with the forgetting factor ${\alpha }_k$ to adaptively update R and Q, by which the adaptability of KF to the abnormal amplitude and the model uncertainty can be adjusted in real time. In addition, the proposed algorithm and three widely used algorithms (PCA, EKF and TS) are applied to 25 groups of AECG generated by the FECGSYN toolbox, 10 groups of clinically collected AECG of pregnant women in the ADFECGDB database, 10 groups of clinical AECG sourced from the FECGDARHA database, and three groups of real AECG collected by the portable and low-cost AECG acquisition hardware system designed in this article. The performance of the above algorithms in the extraction accuracy, SNR and morphological integrity of the extracted FECG are comprehensively evaluated. The results show that the proposed algorithm has better performance in the above three aspects. In the future, the portable and low-cost FECG monitoring system can not only help pregnant women carry out fetal heart monitoring at any time and place more conveniently, but also help doctors more accurately evaluate pregnant women and fetuses’ physiological status and draw up more reasonably formulate health management and clinical treatment plans, which is of great significance to reduce social medical expenses and improve pregnant women and fetuses’ health level as far as possible.

Supplemental Material

sj-docx-1-dhj-10.1177_20552076261435069 - Supplemental material for A novel multi-step adaptive Kalman filtering method based on dynamic noise estimation for FECG extraction

Supplemental material, sj-docx-1-dhj-10.1177_20552076261435069 for A novel multi-step adaptive Kalman filtering method based on dynamic noise estimation for FECG extraction by Yingbin Liu, Longxi Li, Yanbin Guo, Zibo Feng and Guoping Wang in DIGITAL HEALTH