PCA-Based Preprocessing for Clustering-Based Fetal Heart Rate Extraction in Non-Invasive Fetal Electrocardiograms

Citation: Oyarzún, L.; Castillo, E.; Parrilla, L.; Meyer-Baese, U.; García, A. PCA-Based Preprocessing for Clustering-Based Fetal Heart Rate Extraction in Non-Invasive Fetal Electrocardiograms. Electronics 2024, 13, 1264. https://doi.org/10.3390/ electronics13071264 Academic Editors: Constantin Paleologu and Doru Florin Chiper Received: 28 February 2024 Revised: 24 March 2024 Accepted: 26 March 2024 Published: 28 March 2024 Copyright: © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Article PCA-Based Preprocessing for Clustering-Based Fetal Heart Rate Extraction in Non-Invasive Fetal Electrocardiograms Luis Oyarzún 1,2 , Encarnación Castillo 2, Luis Parrilla 2, Uwe Meyer-Baese 3and Antonio García 2,* 1Departamento de Sistemas Computacionales, Universidad Técnica de Manabí, Portoviejo 130105, Ecuador; [email protected] 2 Departamento de Electrónica y Tecnología de Computadores, Universidad de Granada, 18071 Granada, Spain; [email protected] (E.C.); [email protected].es (L.P.) 3Department of Electrical and Computer Engineering, FAMU-FSU College of Engineering, Tallahassee, FL 32310-6046, USA; [email protected] *Correspondence: [email protected]; Tel.: +34-958240474 Abstract: Non-invasive fetal electrocardiography (NI-ECG) is based on the acquisition of signals from electrodes on the mother’s abdominal surface. This abdominal ECG (aECG) signal consists of the maternal ECG (mECG) along with the fetal ECG (fECG) and other noises and artifacts. These records allow the acquisition of valuable and reliable information that helps ensure fetal well-being during pregnancy. This paper proposes a procedure based on principal component analysis (PCA) to obtain a single-channel master abdominal ECG record that can be used as input to fetal heart rate extraction techniques. The new procedure requires three main processing stages: PCA-based analysis for fECG-component extraction, polarity test, and curve fitting. To show the advantages of the proposal, this PCA-based method has been used as the feeding stage to a previously developed clustering-based method for single-channel aECG fetal heart rate monitoring. The results obtained for a set of real abdominal ECG recordings from annotated public aECG databases, the Abdominal and Direct Fetal ECG Database and the Challenge 2013 Training Set A, show improved efficiency in fetal heart rate extraction and illustrate the benefits derived from the use of such a master abdominal ECG channel. This allows us to achieve proper fetal heart rate monitoring without the need for manual inspection and selection of channels to be processed, while also allowing us to analyze records that would have been discarded otherwise. Keywords: fECG; PCA; clustering; fetal monitoring; ADFECG database; Challenge 2013 Training Set A 1. Introduction An essential part of modern obstetric services is electronic fetal monitoring (EFM), which allows diagnosis of fetal distress (FD) during pregnancy. Early detection of any cardiac condition is key to the prevention of after-birth complications or, even, permanent or fatal damage to the newborn’s heart [ 1 , 2 ]. For example, fetal hypoxia is an FD associated with severe perinatal morbidity and mortality, and it is due to an alteration in the placental function that reduces oxygen delivery to the fetus. Thus, the need to prevent heart diseases such as fetal hypoxia has led to the development of different techniques for heart monitoring. However, while this is a well-defined field for adult subjects, fetal monitoring during pregnancy remains a challenge. Auscultation, cardiotocography (CTG), echocardiography, and phonocardiography (PCG) are the most used methods in clinical practice for continuous fetal heart rate (fHR) monitoring [ 3 – 5 ], but those methods have some disadvantages such as their high sensitivity to different types of noise generated by maternal movements and artifacts [ 6 ]. However, even if those techniques, and others based on them, are the most used, they do not suffice to give a conclusive diagnosis due to their limited accuracy (even leading to unnecessary Electronics 2024,13, 1264. https://doi.org/10.3390/electronics13071264 https://www.mdpi.com/journal/electronics Electronics 2024,13, 1264 2 of 31 cesarean sections [ 7 ]). Even if CTG has been considered for years a standard of care for intrapartum fetal surveillance, its data interpretation will always depend on the expert’s judgment and medical diagnosis abilities, since it is based on a visual recognition of the fHR patterns (bradycardia, accelerations, decelerations, etc.), which limits its medical use [ 8 ]. This has caused a growing interest in fetal electrocardiograms (fECG), which have become a promising EFM method. The main reason for this is that the fECG signal carries valuable information, such as changes in the morphology of the fECG waveform associated with dysfunctions induced by FD [ 7 ]. A way to acquire a fECG is the use of a scalp electrode placed on the head of the newborn during labor, but this is a very invasive method. Thus, non-invasive options for fECG acquisition (NI-fECG) during pregnancy and/or labor have been developed [ 1 ]. They make use of some electrodes placed on the maternal abdomen, so each one of these electrodes provides a channel to analyze. While EFM techniques have been shown to help decrease perinatal mortality and morbidity [9] and cesarean sections [10] , NI-fECG can offer superior accuracy and reliability than monitoring based on ultrasound technology, and is less likely to display the maternal heart rate in place of the fetal heart rate [11] . Moreover, NI-fECG does not appear to be influenced by increased body mass index or fetal movement, and its high accuracy of heart-rate indices is promising for NI-fECG usage in remote settings [12] . As a matter of fact, results in [12] suggest that, regardless of signal loss, a trace where at least 10 min of fetal R waves are captured that demonstrate physiological FHR indices by displaying normal key FHR parameters such as short-term variation (STV), should be sufficient to demonstrate fetal wellbeing according to the Dawes–Redman criteria [13] . Additionally, NI-fECG devices deliver no energy and women can potentially wear monitors for longer durations without safety concerns, thus effectively enabling long-term monitoring. However, due to the noise produced by the mother’s and fetus’ movements and other artifacts, the resulting abdominal ECG (aECG) record may be seriously affected, resulting in signals either partially distorted (in some of the channels or in fragments of them) or even completely unusable. The extraction of fHR or other characteristics of the fECG signal thus requires sophisticated signal processing in order to overcome this [ 14 – 16 ]. As an example, in [ 14 ], a clustering-based method for singlechannel fetal heart rate monitoring from aECG was presented. However, in that proposal, not all the channels in aECG recordings could be successfully processed due to artifacts, excessive noise, saturation, etc. In some cases, even complete recordings were useless for fHR extraction as none of their channels allowed the identification of ECG characteristics. The aim of this research is focused on procedures for fHR extraction from NI-fECG recordings that enhance the performance of, in general, previous methods for fHR extraction and, more specifically, those based on clustering [ 14 ]. This improvement is oriented to avoid discarding some fragments of channels in aECG recordings, or even complete channels, due to noise or severe artifact contamination. The final objective is to achieve high accuracy in fHR extraction from any aECG recording, with no need to inspect and select channels or fragments. For this purpose, a representative signal of all the channels in an aECG recording will be composed, in which the noises and artifacts that may affect some channels are minimized or even removed. Principal component analysis (PCA) can be used for this task, since PCA is a technique that can reduce the dimensionality of large data sets by transforming a large set of variables into a smaller one that still contains most of the original information. Thus, in this paper, we propose the use of a PCA-based algorithm to extract a single master signal from the aECG recording, which is then fed to the clustering-based algorithm proposed in [ 14 ]. fHR results from annotated public aECG databases show the benefits associated with this new framework. The rest of the manuscript is organized as follows: Section 2describes the background of the different techniques and algorithms. Section 3is devoted to the description of the proposed algorithm. Results and discussion are presented in Section 4, while the final Section 5summarizes the main conclusions of this work. Electronics 2024,13, 1264 3 of 31 2. Background The fetal heart rate undergoes dynamic variations, reflecting both physiological and pathological influences. fHR changes not related to physiological events may be a sign of pathologies such as fetal hypoxia, so it is a priority to research how to monitor fHR, fetal movements, and uterine contractions simultaneously. In the following, some concepts in this field are summarized, along with the fundamentals of fECG signal processing. 2.1. NI-fECG Fundamentals The ECG is an important technique whose acquisition on the human skin involves the reading of very low voltages generated by the heart activity. These signals are contaminated by different sources of noise, and the difficulty is increased when the acquisition target is the fECG signal. The fHR can be obtained from fECG signals [ 14 , 17 – 19 ], and its monitoring is one of the most important methods to detect and identify diseases in the early stages of pregnancy [ 20 ]. As was commented above, the NI-fECG is a non-invasive method that consists of the use of a set of surface electrodes placed over the mother’s abdomen in order to acquire the aECG, which is composed of the fECG along with the maternal ECG (mECG). Examples of electrode positioning can be found in [ 1 , 2 ]. One of the main advantages of this non-invasive method is the fact that the fetus is not stressed during the acquisition of the signal, and it allows long-term monitoring during pregnancy [ 2 ]. Another advantage of NI-fECG is the absence of any radiated energies, which are mostly present in other methods such as X-rays and ultrasound; although the latter has been shown to be safe for fetal monitoring, prolonged exposure to ultrasound can heat fetal tissues [ 21 ]. Moreover, to obtain fECG signals from aECG signals it is important to consider that they are affected by different types of electrical interferences, such as other electric biological signals (mother’s respiration and contractions, fetus’ movement) and external electric artifacts (power line interferences, etc.). In addition, movements at the electrodes and a strong maternal component [ 22 ] compared to the weak fetal component may complicate fECG detection. Figure 1shows an example of a channel of an aECG recording, where these interferences can be observed. Figure 1. Example of an aECG channel, composed of mECG, fECG, and different noises. NI-fECG acquisition is usually performed according to two aECG acquisition procedures [ 23 ]: the first one is called abdominal electrode sources (AES), which only makes use of the leads placed on the maternal abdomen; meanwhile, the second procedure is called combined source (CS) because, apart from the leads placed on the maternal abdomen, at least an extra lead is placed on the maternal chest to capture the actual mECG signal. This last type of system usually includes an adaptative algorithm that, using those extra leads on the chest, tries to cancel the maternal component, which is the main and most challenging artifact in the aECG mixture [ 2 ] for fHR extraction. The relative weakness of the fetal signal and the noisy environment due to the different artifacts mentioned above may make the use of NI-fECG difficult for detecting structural defects of the fetal heart. However, it is possible to extract the fHR, and it is feasible to evaluate the blood conduction velocity, or the condition of tissues within the heart as well as various abnormalities [ 2 ]. In any case, the identification of QRS complexes and R-R intervals allows for monitoring heart rate variability. Thus, NI-fECG is an ideal tool for long-term FHR monitoring [22]. Electronics 2024,13, 1264 4 of 31 2.2. Processing of NI-fECG Different techniques for processing and extracting characteristics from aECG recordings are found in the literature [ 7 , 24 – 27 ]. As it was mentioned above, the fECG is mixed with other biopotentials and noises, such as the mECG, respiratory artifacts and contractions [ 1 ], and even the movement of the fetus or the mother, and it is very important to consider them as one of the main problems in the analysis of aECG records. Sudden movements from the fetus or the mother can even lead to saturation in the recordings. Additionally, some of these artifacts result in baseline wandering (BW). Wavelet denoising has been shown to be a useful tool for the removal of these artifacts, even suitable for implementation on embedded hardware [ 16 , 28 ]. Other methods of adaptive noise cancellation are least mean square (LMS) and Widrow’s multireference [ 29 – 31 ]. On the other hand, techniques based on wavelet transforms (WT) [ 22 , 32 – 34 ], empirical mode decomposition (EMD) [18,19,35], Kalman filtering [ 25 , 36 ], artificial neural networks [ 37 – 39 ], adaptative filtering [ 40 – 42 ], blind source separation [ 7 , 26 , 30 , 43 ], or even clustering-based algorithms [ 14 ], have been shown to be useful in the extraction of components more suitable to obtain relevant characteristics of the fECG, such as R-peak locations. The techniques to denoise aECG signals and to extract reliable fetus information used in this work are summarized in the following. 2.2.1. Wavelet Denoising The aim of aECG denoising is to recover a clean signal from the original aECG data, allowing further processing such as separation of the mECG and fECG signals, QRS complex detection, or parameter estimation (fHR, for example). The wavelet transform is one of the most used tools in signal processing and is an effective way to digitally remove noises within specific sub-bands for aECG signals [ 16 ]. This wavelet-based approach is used to remove the low frequency trend of a signal, but it introduces no latency and less distortion than digital filter-based approaches. In this way, an adequate wavelet decomposition can be used as a preprocessing method for BW and noise suppression in aECG signals. Moreover, due to the similar wavelet structure for the application of BW and noise suppression, it can be carried out with only one-step wavelet [ 14 , 16 , 22 ], which enhances efficiency without compromising accuracy, thus saving resources and facilitating hardware implementations [22]. In the following, the application of the wavelet preprocessing approach consists of a wavelet decomposition down to L levels, with the approximation coefficients at level L replaced by an all-zero vector. Additionally, for each level from i= 1 to M (with M<L ), the appropriate threshold limit and rule are applied to the detail coefficients. The wavelet reconstruction based on the zeroing approximations of level L , the modified details of levels 1 to M , and the original details of levels M+ 1 to L are computed to obtain the BW-corrected and denoised aECG signal [14,16,22]. 2.2.2. Clustering-Based Procedure for FHR Extraction Cluster analysis is one of the most used statistical methods to separate objects with the same characteristics in diverse groups (called clusters, defined as groups of similar objects that differ from objects in other groups), and is mainly used in data mining and as an unsupervised learning technique for the classification of objects into clusters [14,44]. Thus, clustering can be used to classify distinctive features of denoised aECG signals in order to separate fetal from maternal QRS complexes. After that, fetal QRS complexes can be used for fHR monitoring. The clustering-based method proposed in [ 14 ] uses for this classification the so-called candidate RS-peaks, which are local maxima followed by a local minimum. Concretely, the clustering algorithm comprises two stages: the first stage carries out the detection of maternal and fetal QRS combining the use of thresholds and clustering itself, while the second stage corrects false negatives (FN) and false positives (FP). In the first stage, the maternal and fetal QRS (mQRS and fQRS, respectively) detection is based on the Electronics 2024,13, 1264 5 of 31 location of candidate RS-peaks. These are local maxima followed by a local minimum in the preprocessed aECG signal. These RS-peaks correspond to two possible scenarios [14] , depending on whether or not the amplitudes of the fetal RS waves are similar to the amplitudes of the maternal RS waves, which should be normally larger. When these amplitudes are comparable, an additional feature is required to differentiate maternal and fetal candidates, which is the product of the amplitude and time distance of the candidate peak, rather than the amplitude. Thus, a smoothed version of the normalized distribution of the amplitudes of all candidate peaks in the classification window is classified into three possible cases, depending on the number of local maxima in this distribution [14]: • Case 1: detection of two local maxima from the first local minimum, which represents candidates that are actually noise. Each of these maxima provides information about the amplitude zones corresponding to fetal RS-peaks and maternal RS-peaks. If the distance between these maxima is 35% greater than the maximum candidate amplitude, maternal candidates are larger than fetal ones and amplitudes are the selected data to be classified. • Case 2: detection of two local maxima from the first local minimum, the distance between which is 35% less than the maximum candidate amplitude. This situation can be related to similar amplitudes for fetal and maternal RS-peaks and, thus, amplitudes multiplied by the number of samples of the candidates are the data to be classified. • Case 3: detection of a single local maximum from the first local minimum. This situation generally corresponds again to similar amplitudes for fetal and maternal RS-peaks, and amplitudes multiplied by the number of samples of the candidates are the selected data to be classified. Once the case and, thus, the data to be classified are determined, the candidates are classified using the k-medoids++ clustering algorithm into three groups: maternal R peaks, fetal R peaks, and noise. After this classification, the second stage of the FHR monitoring includes a method for the detection of FNs, non-detected fQRS complexes, and FPs, falsedetected FQRS complexes [22] . For this, several RR time distances are defined for each fQRS candidate to preceding and succeeding candidates, so FPs and FNs are detected imposing limits to the beat-to-beat variability of the detected heart rate [22]. Compared to threshold-based methods [ 22 , 45 ], clustering allows better detection of fetal QRS complexes since it is less affected by the amplitude variation of the R peak. In some cases in [ 14 ] all the channels in a given recording can be analyzed, but in most recordings, one or several channels were discarded due to the quality of the aECG. Thus, while the algorithm in [ 14 ] is able to successfully extract the fHR, it still presents some disadvantages, especially the need to preselect those channels in the recording that are apt for analysis, which usually requires a specialist and is not a trivial task to be automatically carried out. Thus, the objective of this work is to enhance the performance of the cluster algorithm in [ 14 ] for fHR extraction, which will be used to evaluate a single master channel obtained through a PCA-based algorithm of the complete aECG recording. 2.2.3. Blind Source Separation Fundamentals Blind source separation (BSP) methods have been widely used in the past, and they have been successfully used in fECG extraction [ 7 , 27 , 30 , 43 , 46 , 47 ]. The principle of these methods is the decomposition of the signal mixture (aECG) into the original source components (fECG, mECG, and noise), assuming that these components are statistically independent. BSP techniques involve statistical methods such as PCA, independent component analysis (ICA), and singular value decomposition (SVD) [ 23 , 48 ]. These methods have been proven to be suitable for fECG extraction from a multichannel aECG but, nevertheless, a larger number of input signals is associated with greater computational complexity and lower comfort for the pregnant subject during monitoring [ 7 ]. In any case, recordings with more channels will provide more accurate data. Our proposal is based on PCA, whose fundamentals are presented in the following. Electronics 2024,13, 1264 6 of 31 2.2.4. Principal Component Analysis (PCA) PCA is a powerful tool for data analysis, as it finds patterns in data of high dimensionality. It allows us to identify patterns in data and to highlight their similarities and differences, since in other methods these patterns may be hard to find in data of many dimensions [ 40 ]. The other main advantage of the PCA method is that, once data patterns are found, data can be compressed, i.e., the number of dimensions can be reduced without significant loss of information [ 23 ]. Considering the application of PCA methods to aECG signals, they allow us to combine all channels [ 40 ] into a reliable single signal. To perform the PCA analysis, the cross-correlation matrix between input signals is built, eigenvalues and eigenvectors are computed and, finally, input signals are transformed according to the eigenvectors and eigenvalues [40]. The covariance matrix, mtcov, can be calculated as: mtcov[xi,xj] = E(xi−E[xi])(xj−E[xj])T(1) where x is the array with the values to be analyzed, E denotes the expected value (mean), and T is the symbol of the transposed matrix. Let e vectork be the k -th eigenvector and evaluesk the corresponding eigenvalue of mtcov . Then, the k -th principal component is obtained using the following expression: yk=evectorT kx(2) The spread of the eigenvalues is a good indicator of how many principal components are of real interest in our input data, so the signals associated with small eigenvalues are usually dropped. As a result, we will have one or more components that represent the input signals in a minimum mean square sense. The principal components capture the power of most of the signal and, at the same time, all components are orthogonal to each other. That is why PCA sometimes is also called a whitening operation [ 15 , 40 , 49 ]. Thus, the PCA method fits perfectly with one of the main objectives of this work, which is the generation of a single aECG master signal without discarding any of the original channels from the aECG recording, in order to subsequently apply techniques for the extraction of fetal information of interest such as the fHR. 3. Proposed Framework As discussed extensively above, NI-fECG consists of a set of multiple signals obtained with electrodes placed on the maternal abdomen. Usually, most proposals to extract the fetal characteristics [ 7 , 14 , 15 , 48 ] require the search for the channel or channels from the aECG with less noise contamination, discarding even fragments of some channels when they are affected by severe artifacts [ 40 ], and even discarding complete aECG recordings. Thus, removing or discarding any channel could remove vital information about the fetus. The proposed PCA-based algorithm seeks to solve this problem and to prevent a complete record, signal, or fragment from being discarded. With that aim, a PCA-based [ 15 , 23 , 27 , 40 ] algorithm is proposed, which will extract a single aECG master signal from the aECG record without discarding any of its channels. This master signal will be fed to the clustering algorithm in [ 14 ] to extract the fetal QRS complexes and, thus, monitor fHR. The proof of concept has been implemented using GNU Octave 8.4 and Python 3.12, due to library compatibility, and then trained and tested using aECG recordings in two PhysioNet databases [ 50 ]. These databases are described below, while the new procedure will be detailed in Section 3.2. 3.1. Fetal ECG Datasets PhysioNet [ 51 ] offers free and public web access to large collections of recorded physiologic signals (PhysioBank), and the included databases are made available under the ODC Public Domain Dedication and License v1.0 [ 52 ]. Concretely, two of these databases provided by PhysioNet have been used for this work: Electronics 2024,13, 1264 7 of 31 • Abdominal and Direct Fetal Electrocardiogram Database (ADFECGDB) [ 50 ]: this database contains multichannel fECG recordings obtained from five different women in labor. Each recording comprises four 5 min differential signals acquired from the maternal abdomen, and the reference direct fECG registered from the fetus head. The recordings are sampled at 1 ksps with 16-bit resolution, and the signal bandwidth is 1–150 Hz. Moreover, the database includes a set of reference annotations indicating the fetal R-wave locations. The ADFECGDB will be used to train the proposed algorithm. • Challenge 2013 Training Set A (Challenge) [ 53 ]: these data consist of one-minute fECG recordings, sampled at 1 ksps, each one including four noninvasive abdominal signals as well as the reference annotations marking fetal R-wave locations [ 14 ]. The Challenge database is used to test and validate the proposed algorithm. 3.2. Proposed PCA-Based Framework As it was commented above, the objective of this work is to generate, through a PCAbased algorithm, an aECG master signal from the different channels in an aECG recording. As it is common practice in many digital signal processing algorithms, the signals will be divided into windows. Thus, the size of these windows must be adjusted to provide a correct analysis as well as to minimize the delay between acquisition and processing. It is also important to take into account that window-based processing can produce noncontinuous signals, while it is also possible to have opposing polarities in adjacent windows as a result of differing signs in the eigenvalues. Considering all these aspects, this work proposes to track the signal obtained after the PCA analysis in order to perform a polarity check followed by a curve fitting. A schematic of the proposed PCA-based framework is shown in Figure 2: PCA-based analysis, polarity check, and curve fitting. These stages and their steps are described in detail below, keeping first in mind that the aECG signals to be processed will have been preprocessed using the one-step wavelet-based algorithm described in [ 16 ]. Figure 3illustrates a fragment of the r01 record from the ADFECGDB training set before and after this wavelet-based preprocessing. The parameters [ 16 ] for the wavelet preprocessing stage of the databases used in this work were: wavelet function Daubechies 6, L = 7 (due to the maximum frequency component of the signals), M = 3, universal threshold, multiple rescaling, and soft thresholding. 3.3. Stage 1: PCA-Based Analysis for Fetal ECG Extraction Due to the presence of artifacts in aECG signals, many works in the literature [ 14 , 15 , 23 , 27 , 35 , 40 ] usually discard a fragment of a channel from aECG recordings, a complete channel, or even all of them, when extracting fetal characteristics such as fHR. Therefore, the application of the proposed PCA-based algorithm before extracting any information allows us to obtain a master aECG signal as the first principal component using the entire aECG record, including those fragments and channels that might have been discarded in previous models due to noise or artifacts. The application of PCA to the input signal avoids the need for an a priori analysis of the aECG records, searching for fragments in the channels that are not affected by artifacts or other noise. This limits the loss of information and provides a feasible signal for further analysis, even for studies outside of a controlled environment. A complete scheme of this stage with the three necessary steps is shown in Figure 4, while an algorithmic summary is shown in Algorithm 1 . This stage is detailed in the following. Electronics 2024,13, 1264 8 of 31 One-step wavelet-based preprocessing Stage 1: PCA-based analysis for fetal ECG extraction Step 1: Define window size Step 2: Remove offset Step 3: Perform PCA Stage 2: Polarization Check Step 1: Search for max-min-max-min points Step 2: Distance comparison Step 3: Determine correct QRS polarization Stage 3: Curve Fitting Step 1: Continuity correction Clustering for fHR monitoring Figure 2. Proposed processing framework. (a) (b) Figure 3. (a) Fragment of channel #1 in record r01 from ADFECGDB; (b) resulting signal after one-step wavelet-based preprocessing. Electronics 2024,13, 1264 9 of 31 Stage 1: PCA-based analysis for fetal ECG extraction Step 1: Define window size • Window size to be analyzed (ws, in blue) • Overlap window size (os, in red) • Total window size: os +ws +os Step 2: Remove offset Step 3: Perform PCA • Find covariance matrix (mtcov) of aECG records • Find eigenvalues and eigenvectors (evector,evalues) • Apply PCA evectorT×aECG √evalues Figure 4. Stage 1: PCA-based analysis for fetal ECG extraction. Algorithm 1 PCA-based analysis for fetal ECG extraction. 1: Step 1: Define window size 2: if signal is shorter than ws then ▷One-window case 3: left window extension ←0 4: right window extension ←0 5: else if first window then ▷First window 6: left window extension ←0 7: right window extension ←os 8: else if last window then ▷Last window 9: left window extension ←os 10: right window extension ←0 11: else ▷Any other case 12: left window extension ←os 13: right window extension ←os 14: end if 15: aux ←aECG[extended window(ws,os)] 16: Step 2: Remove offset 17: for j=1to num_channels do 18: aux[j, :]←aux[j, :]−mean(aux[j, :]) 19: end for 20: Step 3: Perform PCA 21: mt_cov ←cov(aux′)▷Find covariance matrix 22: [evectors,evalues]←eig(mt_cov)▷Find eigenvalues and eigenvectors 23: pca ←(evalues−0.5)×(evectors′)×(aux)▷PCA 24: Master ←pca[first_principal_component, :]▷aECG master signal 3.3.1. Step 1: Define Window Size In order to reduce the complexity of the proposed signal processing, especially if its future implementation on wearable devices is sought, it is important to select an appropriate window size ( ws ) that provides a balance between performance (in terms of data rate, resources, and execution time) and the quality of the resulting signal. Since the training database records, PhysioNet’s ADFECGDB [ 50 ], have a total of 300,000 samples, window Electronics 2024,13, 1264 16 of 31 It is important to finally note that it is not necessary to include additional scenarios for a possible T wave, since they will result in redundant evaluations. Thus, these four scenarios will suffice to determine the type of wave, as discussed below, that the pair of maxima and minima may match. The pseudocode in Algorithm 2 summarizes this first step of the second stage of the proposed PCA-based framework. Algorithm 2 Polarization check, step 1. 1: Step 1: Search for max-min-max-min points 2: Consecutive pairs of maximum followed by minimum 3: procedure MAXMIN(Master) 4: flag ←1, j←1 5: for mx ←2to length[Master]do 6: if (Master[mx]≥Master[mx −1]) and (Master[mx]>Master[mx +1]) and (mx ≤length[Master]) and (flag =1)then ▷Find the first maximum 7: max1[j]←mx 8: flag ←2 9: end if 10: if (Master[mx]≤Master[mx −1]) and (Master[mx]<Master[mx +1]) and (mx ≤length[Master]) and (flag =2)then ▷Find the first minimum 11: min1[j]←mx 12: flag ←3 13: end if 14: if (Master[mx]≥Master[mx −1]) and (Master[mx]>Master[mx +1]) and (mx ≤length[Master]) and (flag =3)then ▷Find the second maximum 15: max2[j]←mx 16: flag ←4 17: end if 18: if (Master[mx]≤Master[mx −1]) and (Master[mx]<Master[mx +1]) and (mx ≤length[Master]) and (flag =4)then ▷Find the second minimum 19: min2[j]←mx 20: flag ←1 21: j←j+1 22: end if 23: end for 24: end procedure 25: 26: Based on amplitude distances dPQ,dQR,dRS and ST, 27: compute four distance sets, related to four scenarios 28: procedure QRS(Master,max1,max2,min1,min2) 29: for j←2 to length[max1]do 30: ▷P-wave consideration 31: dPQ1[j−1]←Master[max1[j]] −Master[min1[j]] 32: dQR1[j−1]←Master[max2[j]] −Master[min1[j]] 33: dRS1[j−1]←Master[max2[j]] −Master[min2[j]] 34: dST1[j−1]←Master[max1[j+1]] −Master[min2[j]] 35: ▷R-wave consideration 36: dPQ2[j−1]←Master[max2[j−1]] −Master[min2[j−1]] 37: dQR2[j−1]←Master[max1[j]] −Master[min2[j−1]] 38: dRS2[j−1]←Master[max1[j]] −Master[min1[j]] 39: dST2[j−1]←Master[max2[j]] −Master[min1[j]] 40: ▷Q-wave consideration 41: dPQi1[j−1]←Master[max1[j]] −Master[min2[j−1]] 42: dQRi1[j−1]←Master[max1[j]] −Master[min1[j]] 43: dRSi1[j−1]←Master[max2[j]] −Master[min1[j]] 44: dSTi1[j−1]←Master[max2[j]] −Master[min2[j]] 45: ▷S-wave consideration 46: dPQi2[j−1]←Master[max2[j−1]] −Master[min1[j−1]] 47: dQRi2[j−1]←Master[max2[j−1]] −Master[min2[j−1]] 48: dRSi2[j−1]←Master[max1[j]] −Master[min2[j−1]] 49: dSTi2[j−1]←Master[max1[j]] −Master[min1[j]] 50: end for 51: end procedure 3.4.2. Step 2: Distance Comparison Following the considerations above on the four possible scenarios maxima and minima may match, all four sets of amplitude distances described in the previous section are evaluated for each pair of maximum–minimum points. Thus, the max1(j) , min1(j) , max2(j) , Electronics 2024,13, 1264 17 of 31 and min2(j) points in the j -th pair of maxima and minima, along with any required points from the (j− 1 ) -th and (j+ 1 ) -th set of pairs as illustrated in Figures 10–13, are used to evaluate the four different sets of amplitude distances for this j -th pair: dPQ1 , dQR1 , dRS1 , and dST1 (Figure 10), dPQ2 , dQR2 , dRS2 , and dST2 (Figure 11), dPQi1 , dQRi1 , dRSi1 , and dSTi1(Figure 12), and dPQi2 , dQRi2 , dRSi2 , and dSTi2 (Figure 13). It must be noted that when the pair of maxima and minima does not correspond to any cardiac component but is caused by noise and artifacts, as will usually be the case for most pairs of maxima and minima, none of the amplitude distance sets will match any of the scenarios. For this analysis, a set of conditions and thresholds are thus imposed to all the sets of amplitude distances in order to find, if any, the matching scenario. It must be kept in mind that it has been considered that the largest distances will correspond to the QRS complex, with dQR <dRS. Additionally, thresholds are imposed to avoid any distortion effects derived from the PCA processing. With this all, the following conditions must hold for the matching scenario: dPQ <dQR (4) dRS >dST (5) dQR >0.5 and dRS >1.0 (6) where the thresholds in (6) have been manually set after the training phase of the algorithm. The four sets of distances, for each pair of maxima and minima, are compared to the conditions in Equations (4)–(6), so if one of the distance sets is found to satisfy those conditions, the pair of maxima and minima is classified in the corresponding scenario. Thus, it is possible to compute over two different counters the number of detected QRS complexes that are correctly and incorrectly polarized, respectively. Algorithm 3 summarizes the procedure for this distance comparison. Algorithm 3 Polarization check, step 2. 1: Step 2: Distance comparison 2: procedure DISTANCE(distances) 3: k←1, l←1 4: for j←2to length[dQR1]do 5: if (dPQ1<0.5)and (dQR1>1)and (dRS1>1)and (dPQ1<dQR1) and (dRS1>dST1)then 6: q_peakn[k]←min1[j] 7: r_peakn[k]←max2[j] 8: s_peakn[k]←min2[j] 9: k←k+1 10: else if (dPQ2<0.5)and (dQR2>1)and (dRS2>1)and (dPQ2<dQR2) and (dRS2>dPQ2)then 11: q_peakn[k]←min2[j−1] 12: r_peakn[k]←max1[j] 13: s_peakn[k]←min1[j] 14: k←k+1 15: else if (dPQi1<0.5)and (dQRi1>1)and (dRSi1>1)and (dPQi1<dQRi1) and (dRSi1>dSTi1)then 16: q_peaki[l]←max1[j] 17: r_peaki[l]←min1[j] 18: s_peaki[l]←max2[j] 19: l←l+1 20: else if (dPQi2<0.5)and (dQRi2>1)and (dRSi2>1)and (dPQi2<dQRi2) and (dRSi2>dSTi2)then 21: q_peaki[l]←max2[j−1] 22: r_peaki[l]←min2[j−1] 23: s_peaki[l]←max1[j] 24: l←l+1 25: end if 26: end for 27: end procedure Electronics 2024,13, 1264 18 of 31 3.4.3. Step 3: QRS Polarization Correction It is now possible to adjust the polarity of the signal in case most of the QRS-complexes are inverted according to the analysis in the previous section, as summarized in Algorithm 4 . In this case, the whole window will be corrected, as illustrated in Figure 14. Algorithm 4 Polarization check, step 3. 1: Step 3: QRS Polarization Correction 2: procedure QRSCORRECTION(Master,r_peakn,r_peaki) 3: if length[r_peakn]<length[r_peaki]then 4: Invert aECG Master Signal 5: end if 6: end procedure (a) (b) Figure 14. QRS polarization correction: (a) fragment of the PCA-derived aECG master signal from r01 record with wrong polarity (red triangles mark R-peaks); (b) resulting fragment after polarity correction (blue triangles mark R-peaks). 3.5. Stage 3: Curve Fitting Due to the nature of windowed signal processing, continuity between adjacent windows is not guaranteed and there may be small offsets between them. Thus, curve fitting is the final stage of the proposed algorithm, composing a continuous signal by adjusting each of the windows that have been processed to the adjacent ones. This will result in the final record to be analyzed through clustering. A schematic of this stage is shown in Figure 15. As an example, Figure 16a shows how there is a discontinuity between the resulting signals from windows 3 (in blue) and 4 (in red) of record r01 of the ADFECGDB database after the PCA analysis and polarity correction described above. To solve this, the use of the parameter os as defined in Stage 1 creates an overlapping region between these two windows, as can be seen in Figure 16b. This overlapping region between adjacent windows allows to compute a weighted average (in green in Figure 16c) from the two different signals, where continuity can be enforced by making the weight of each component at each sample to be proportional to the sample position within the overlapping region. In this way, at the beginning of the overlapping region, the signal from the preceding window ( Pw ) is considered with unity weight, which decreases to zero at the end of the region; on the other hand, the weight from the component from the following window ( Cw ) moves from zero, at the beginning of the overlapping region, to one at the end of that region. Electronics 2024,13, 1264 19 of 31 Stage 3: Curve fitting. Step 1: Continuity correction • Define initial point from previous window Pw :x=nfp−os • Define final point from current window Cw :y=nip+os • Overlap range between previous and current windows: [x,y] • Size of the overlapping region: y−x • for iin [x,y] Avgi=[Pwi×(y−i)−Cwi×(x−i)] y−x Figure 15. Stage 3: Curve fitting. Concretely, and using the notation introduced in Figure 6, the overlapping region is defined by its initial, x, and final, y, sample positions: x=nfp−os (7) y=nic+os (8) where nic is the initial sample of the current window ( Cw ) and nfp is the final sample of the previous window ( Pw ). With this, the curve fitting procedure is computed as a weighted average in the [ x , y ] range. Thus, the i -th sample of this weighted average, Avgi , is computed as: Avgi=[Pwi×(y−i)−Cwi×(x−i)] y−x(9) where Pwi is the value of the component from the preceding window (in blue in Figure 16b) and Cwi is the value of the component from the current window (in red in Figure 16b). This all is summarized in Algorithm 5. Algorithm 5 Curve fitting. 1: Step 1: Continuity correction 2: procedure CURVEFIT(Pw,Cw,os) 3: x←Pw :nfpsubtracting os 4: y←Cw :nipadding os 5: for i←xto ydo 6: Master[i]←Pw[i]×(y−i)−Cw[i]×(x−i) y−x 7: end for 8: end procedure Electronics 2024,13, 1264 20 of 31 (a) (b) (c) Figure 16. (a) Discontinuity between preceding window ( Pw , in blue) and current window ( Cw , in red) without overlap; (b) overlapping region between Pw (in blue) and Cw (in red) defined by os ; and (c) curve fitting (in green) between Pw (in blue) and Cw (in red). 4. Validation and Results As previously described, the objective of the proposed framework is to produce an aECG master signal from the PCA-based analysis above that will be fed to the clustering algorithm proposed in [ 14 ]. It is important to note that in [ 14 ] a study to establish the optimal value of the number of samples in the data window for the clustering classification was carried out. This study indicated that, in general, there was no noticeable difference in the effectiveness of the classification results for windows from 10,000 to 60,000 samples. Thus, 50,000and 60,000-sample windows were selected in [ 14 ] for PhysioNet’s ADFECGDB [ 50 ] records and for PhysioNet’s Challenge 2013 Training Set A [ 53 ] records, respectively. In this work, a new study has been carried out in order to determine the optimal value of the number of samples in the data window for the PCA-based analysis described above and for the clustering-based fHR extraction, considering different ws values for each method. For this study, PhysioNet’s ADFECGDB records have been used to train the algorithm and to derive the optimum parameter values, while the Challenge 2013 Training Set A has been used as a test and validation database with these parameters determined during the training phase. This will be described in the following subsections after defining the metrics used to evaluate the performance of the proposed framework. 4.1. Performance Metrics In [ 14 ] the different records were evaluated through the analysis of individual channels, discarding those presenting artifacts and/or noisy fragments, since these provided poor Electronics 2024,13, 1264 21 of 31 results. For the evaluation of results [ 14 ], values such as sensitivity ( Se ), positive diagnostic value (PPV), accuracy (Acc), and F1[24] were used: Se =TD TD +FN (10) PPV =TD TD +FP (11) F1=2×PPV ×Se PPV +Se (12) where: F1=2×TD 2×TD +FN +FP. (13) In these equations, TD represents true-detected fetal QRS complexes, and FP and FN are false negatives and false positives, respectively [ 14 ]. In the following, F1 will be used as the control and comparison parameter. 4.2. Training The algorithm in [ 14 ] was trained using PhysioNet’s ADFECGDB records, which include fetal R-wave markers that allow the calculation of different accuracy and performance parameters. For the training of this algorithm [14] , the different recordings were used as a training set to establish the thresholds separating the different min-max scenarios for the signal, the optimal value of the data window size for the clustering classification, the number of times to repeat clustering using the new initial cluster medoid positions, and some other parameters. The results were classified into two groups, where group 1 included the records that a specialist doctor had previously selected and that were visually suitable for analysis. Group 2 included records with noisy fragments or artifacts, which led to a less efficient clustering-based identification of fetal heartbeats. As described above, the presented algorithm has been also trained using the five ADFECGBD database records, obtaining satisfactory results without the need to select suitable channels since a single aECG master channel resulting from the PCA analysis was fed to the clustering-based fHR monitor. In order to evaluate this approach in comparison to the analysis in [ 14 ], three control values have been considered for the F1 results: the first one is the maximum value of F1 ( Max_F1 ) obtained in any of the channels of the register evaluated in [ 14 ], the second one is the average value of the F1 data obtained from all channels in group 1 ( Mean_Group1 ), and the third control value is the average F1 from all channels in both group 1 and group 2 ( Mean_F1 ). All of these values were calculated using an analysis window of 50,000 samples in [ 14 ]. The results provided by the proposed algorithm and their respective comparison to the results in [ 14 ] are presented in Tables 1–5and in Figures 17–21 for different clustering windows sizes. Concretely, Table 1and Figure 17 show data with 15,000-sample windows for clustering, Table 2and Figure 18 correspond to 20,000-sample windows, Table 3and Figure 19 show data for 30,000-sample windows, Table 4and Figure 20 correspond to 50,000-sample windows and, finally, Table 5and Figure 21 show data for 60,000-sample windows. To differentiate values and understand them, italics are used in these tables to highlight when the ws values used in the PCA analysis are multiples of the ws used in the clustering classification, the grey box is used to highlight when the ws of the PCA and the clustering classification are the same, and the bold typeface is used to highlight the best F1 values. On the other hand, the last two rows of these tables show the discarded channels and noisy channels in [ 14 ]. Furthermore, the values of Max_F1 , Mean_Group1 , and Mean_F1 used in the Figures correspond to the maximum value of all the Max_F1 results, the average value of all the Mean_Group1 values, and the minimum value of all the Mean_F1 results that were obtained in [ 14 ], respectively. Finally, the green dots in the Figures indicate when the PCA ws is a multiple of the clustering window size in the analysis and the red squares mark when they are not multiples. Electronics 2024,13, 1264 22 of 31 Table 1. Results of clustering for fHR monitoring with clustering windows of 15,000 samples (results in italics indicate when the PCA ws is a multiple of the clustering window size, results in bold correspond to the best results for each record, grey box highlights coincident window sizes for both PCA and clustering). ADFECGDB Record r01 r04 r07 r08 r10 PCA ws F1(%) F1(%) F1(%) F1(%) F1(%) Training results derived from the proposed framework 15,000 93.10 87.4 91.94 89.35 94.98 20,000 97.42 97.77 98.56 91.45 94.29 30,000 98.60 97.06 99.60 84.52 95.29 50,000 96.8 96.72 99.60 92.16 94.98 60,000 98.68 97.13 99.60 82.59 95.14 100,000 99.3 96.15 99.52 89.39 94.83 150,000 99.46 97.05 99.92 82.76 94.98 300,000 99.38 97.61 99.76 90.52 94.98 Result summary from [14] Max_F199.38 97.62 98.96 99.08 98.09 Mean_Group198.91 97.41 97.99 97.86 98.02 Mean_F193.39 95.19 97.99 89.72 95.95 Discarded channels – 1 1 – 3 Noisy channels 2, 3 3 – 2, 3 4 Figure 17. Plot of F1 training results for ADFECGDB with clustering windows of 15,000 samples (green dots correspond to PCA ws that are multiples of the clustering window size). Electronics 2024,13, 1264 23 of 31 Table 2. Results of clustering for fHR monitoring with clustering windows of 20,000 samples (results in italics indicate when the PCA ws is a multiple of the clustering window size, results in bold correspond to the best results for each record, grey box highlights coincident window sizes for both PCA and clustering). ADFECGDB Record r01 r04 r07 r08 r10 PCA ws F1(%) F1(%) F1(%) F1(%) F1(%) Training results derived from the proposed framework 15,000 91.21 88.49 99.12 88.45 95.29 20,000 99.15 96.90 97.16 82.39 94.44 30,000 98.91 96.75 98.31 78.14 93.47 50,000 96.76 95.53 98.31 85.16 93.47 60,000 98.76 97.54 98.72 84.99 94.44 100,000 99.46 96.26 99.44 88.68 95.29 150,000 99.53 94.44 99.28 86.15 94.29 300,000 99.46 95.65 99.68 92 95.29 Result summary from [14] Max_F199.38 97.62 98.96 99.08 98.09 Mean_Group198.91 97.41 97.99 97.86 98.02 Mean_F193.39 95.19 97.99 89.72 95.95 Discarded channels – 1 1 – 3 Noisy channels 2, 3 3 – 2, 3 4 Figure 18. Plot of F1 training results for ADFECGDB with clustering window of 20,000 samples (green dots correspond to PCA ws that are multiples of the clustering window size). Electronics 2024,13, 1264 24 of 31 Table 3. Results of clustering for fHR monitoring with clustering windows of 30,000 samples (results in italics indicate when the PCA ws is a multiple of the clustering window size, results in bold correspond to the best results for each record, grey box highlights coincident window sizes for both PCA and clustering). ADFECGDB Record r01 r04 r07 r08 r10 PCA ws F1(%) F1(%) F1(%) F1(%) F1(%) Training results derived from the proposed framework 15,000 93.10 87.40 91.94 89.35 94.98 20,000 97.42 97.77 98.56 91.45 94.29 30,000 98.6 97.06 99.60 84.52 95.29 50,000 96.80 96.72 99.6 92.16 94.98 60,000 98.68 97.13 99.6 82.59 95.14 100,000 99.30 96.15 99.52 89.39 94.83 150,000 99.46 97.05 99.92 82.76 94.98 300,000 99.38 97.61 99.76 90.52 94.98 Result summary from [14] Max_F199.38 97.62 98.96 99.08 98.09 Mean_Group198.91 97.41 97.99 97.86 98.02 Mean_F193.39 95.19 97.99 89.72 95.95 Discarded channels – 1 1 – 3 Noisy channels 2, 3 3 – 2, 3 4 Figure 19. Plot of F1 training results for ADFECGDB with clustering window of 30,000 samples (green dots correspond to PCA ws that are multiples of the clustering window size). Electronics 2024,13, 1264 25 of 31 Table 4. Results of clustering for fHR monitoring with clustering windows of 50,000 samples (results in italics indicate when the PCA ws is a multiple of the clustering window size, results in bold correspond to the best results for each record, grey box highlights coincident window sizes for both PCA and clustering). ADFECGDB Record r01 r04 r07 r08 r10 PCA ws F1(%) F1(%) F1(%) F1(%) F1(%) Training results derived from the proposed framework 15,000 89.06 88.56 98.49 90.98 94.76 20,000 91.71 91.98 86.90 93.08 95.22 30,000 88.25 97.77 99.28 84.54 95.06 50,000 99.46 97.37 99.68 89.79 95.22 60,000 87.81 96.5 99.68 86.75 94.76 100,000 99.22 97.13 99.84 88.17 95.22 150,000 99.22 97.37 99.84 87.01 95.22 300,000 99.07 96.48 99.84 89.73 95.22 Result summary from [14] Max_F199.38 97.62 98.96 99.08 98.09 Mean_Group198.91 97.41 97.99 97.86 98.02 Mean_F193.39 95.19 97.99 89.72 95.95 Discarded channels – 1 1 – 3 Noisy channels 2, 3 3 – 2, 3 4 Figure 20. Plot of F1 training results for ADFECGDB with clustering window of 50,000 samples (green dots correspond to PCA ws that are multiples of the clustering window size).