Fetal QRS detection and heart rate estimation: a wavelet-based approach

Electronic fetal heart rate (FHR) analysis, introduced into clinical practice about 40 years ago, is now the most widely used fetal monitoring technique in industrialized countries (FIGO 1995) and has provided extensive knowledge on intrauterine oxygenation. The most prominent scientific associations (Rooth et al 1987, FIGO 1995, RCOG 2001, ACOG 2005) recommend the surveillance of selected pregnancies from 24–26 weeks of gestation to term and fetal monitors are widespread in modern obstetric units.

There are systems for computerized and automatic analysis of FHR, which may provide real-time alerts for healthcare professionals when changes associated with fetal hypoxia are detected (Nunes et al 2013). Research has been performed on the application of alternative linear and nonlinear FHR indices, which have proven to be useful to detect cases of lower umbilical artery blood pH (Gonçalves et al 2006), particularly in intrauterine growth restricted fetuses when additional information such as the gender of the fetus is also considered (Gonçalves et al 2013). However, despite recent advances, there are still inherent technical limitations (Sameni and Clifford 2010, Clifford et al 2014), bounding its contribution to the improvement on perinatal indicators (Nunes et al 2013).

One of the main challenges in this field is to extract accurate and useful information from the external fetal electrocardiogram (FECG), which may provide a better non-invasive characterization of the fetal cardiovascular system during the third trimester of pregnancy. A more accurate detection of the fetal cardiac rhythms, from the maternal abdominal ECG, is a current research topic (Silva et al 2013, Clifford et al 2014). Such improvements in FHR extraction may allow better performance of currently used FHR indexes, as well as the application and development of alternative FHR indexes.

There are several methods for QRS detection in related literature, among which wavelet transform (WT) based strategies can be found (Elgendi 2014). The use of a derivative WT for QRS detection over standard ECG signals was proposed by Li et al (1995). An extended and enhanced algorithm was later developed and validated in Martínez et al (2004), with good results over standard databases. The WT provides a description of the signal in the time-scale domain, allowing the representation of its temporal features at different resolutions (scales) according to their frequency content. Thus, regarding the purpose of locating different waves with typical frequency characteristics, avoiding noise and artifacts, the WT seems a suitable tool for QRS location over FECG.

Since the maternal abdominal ECG is composed of the maternal and fetal ECGs, it becomes obvious to apply methods of source separation such as independent component analysis (ICA). One of the main limitations of ICA is the difficulty of automatically assigning each of the ICA components to the maternal or fetal ECG (Hyvárinen and Oja 2000, Sameni et al 2006).

An adapted version of the algorithm described in Martínez et al (2004), focusing on fetal QRS (FQRS) detection over abdominal FECG recordings, which allows for location of both maternal and fetal QRS complexes was developed as part of the Noninvasive Fetal ECG PhysioNet/Computing in Cardiology Challenge 2013 (CinC_Ch) (Almeida et al 2013, Vidaurre et al 2011). A wide variety of techniques have been used in an attempt to locate FQRS from abdominal FECG and WT have been used typically for denoising (Vidaurre et al 2011). In our approach WT is explicitly employed for detection by a set of filters especially suited to singularities location, using the same strategy as for an adult ECG. The objective of this work is to present a enhanced version of that described in Almeida et al (2013) and fully validate it. The algorithm was applied to ECG leads, as well as to derived ICA components, from two different databases.

Two datasets of four abdominal FECG leads with reference FQRS annotations available in PhysioBank of PhysioNet (Goldberger et al 2000) were used, denoted here by Silesia dataset and Challenge dataset. The first dataset also includes a channel with the simultaneous direct FECG. The proposed approach is shown as a scheme in figure 1. The general approach for abdominal (indirect) FECG consisted of four steps:

• (a)  
  initial pre-processing of each lead j(  j = 1, 2, 3, 4);
• (b)  
  application of the original QRS detection strategy of Martínez et al (2004) to detect maternal QRS complexes, locating their boundaries and identifying associated WT's extrema (maximum modulus lines—MML);
• (c)  
  application of the adapted QRS detection strategy to identify the FQRS over abdominal FECG (SL_j, j = 1, 2, 3, 4), which excludes maternal QRS complexes using the information from step (b);
• (d)  
  combination of single leads marks, aiming to improve FQRS identification (SLR).

An alternative approach, represented in the gray box of figure 1, consisted in feeding the ICA components, obtained after step (a), to step (c), leading to another version of identified FQRS (ICA_k, k = 1, 2, 3, 4). For the Silesia dataset, the original QRS detection was also directly applied to preprocessed FECG, leading to FQRS (direct).

Zoom In Zoom Out Reset image size

2.1. Data and pre-processing

• Silesia dataset Data from the Abdominal and Direct Fetal Electrocardiogram Database consisting of five files (r01, r04, r07, r08 and r10) of 5 min, sampled at 1000 Hz recorded from women in labor between 38 and 41 weeks of gestation (Matonia et al 2006, Kotas et al 2011). Each file includes one channel with a direct FECG, acquired from a fetal scalp electrode, and 4-leads of noninvasive abdominal FECG signals. The recordings were acquired in the Department of Obstetrics at the Medical University of Silesia, by means of the KOMPOREL system for acquisition and analysis of FECG (ITAM Institute, Zabrze, Poland). The reference marks for FQRS locations are available for all five recordings. Initial R-wave locations automatically determined in the direct FECG signal using the KOMPOREL system were verified by a group of cardiologists. The reference annotations have been stored together with the signals in EDF/EDF + format and were accessed using BioSig for Octave and Matlab (biosig4octmat) (Vidaurre et al 2011).
• Challenge dataset Data from the sets A and B of the CinC_Ch consisting, respectively, of 75 and 100 files of 1 min 4-leads of noninvasive abdominal ECG signals, sampled at 1000 Hz (Vidaurre et al 2011, Clifford et al 2014). These files are part of a larger set, which includes Silesia dataset files and comprises a total of 447 files, divided by three sub datasets named A, B and C. ECG samples corresponding to invalid observations present the special value − 32 768. Reference marks for FQRS complex locations are available only for 74 of set A files, while set B is used for blind testing (Challenge organizers scoring), as described in Clifford et al (2014). Both signals and annotations are in MIT-BIH ECG data format and were accessed and preprocessed using the BioSigBrowser Matlab tool (Bolea et al 2009).

Samples with the special value − 32 768 in the Challenge dataset files were ignored in further processing. A reduced bandwidth notch filter around frequency ω₀ = 0, according to the system function:

$$\begin{eqnarray}&&H\left(z\right)={{b}_{0}}\frac{1-2 \cos {{\omega}_{0}}{{z}^{-1}}+{{z}^{-2}}}{1-2r \cos {{\omega}_{0}}{{z}^{-1}}+{{r}^{2}}{{z}^{-2}}}\end{eqnarray} \tag{ 1 }$$

with b₀ = 1 and r = 0.95, was applied to all signals for detrending.

Independent component analysis (ICA) was applied to abdominal FECG recordings using the fast fixed-point algorithm implemented in the FastICA package for MATLAB (Hyvárinen and Oja 2000) to obtain ICA derived leads.

2.2. Wavelet-based detection method

A single-lead based delineation system using the WT combined with a derivative prototype wavelet (Martínez et al 2004) is used here for QRS detection. For the selected prototype wavelet the WT is implemented using FIR low and high-pass filters with impulse responses

$$\begin{eqnarray}&&h\left[n\right]=1/8\cdot \left\{\delta \left[n+2\right]+3\delta \left[n+1\right]+3\delta \left[n\right]+\delta \left[n-1\right]\right\}g\left[n\right]=2\cdot \left\{\delta \left[n+1\right]-\delta \left[n\right]\right\}.\end{eqnarray} \tag{ 2 }$$

Using these wavelet filters, the WT at scale 2^m, w_{x, m}[n] is proportional to the derivative of the filtered version of the signal x[n] with a smoothing impulse response at scale 2^m. Thus, signal wave peaks correspond to zero crossings in the WT and signal maximum slopes correspond to the WT maxima and minima (maximum modulus lines—MML), as can be seen in figure 2 where direct and abdominal FECG signals and the respective WT signals are plotted.

Zoom In Zoom Out Reset image size

The detection of the fiducial points is carried out across the adequate WT scales and attending to the dominant frequency components. QRS waves are located across scales 2¹ to 2⁴, by searching candidates for MML as local extreme points over (root mean square based) scale dependent thresholds. Isolated and redundant candidates are eliminated and only the ones that appear as MML pairs of opposite polarity across adequate scales are considered. Any MML line is considered to be isolated and discarded if at scale 2¹ differs 0.15 s or more from the closest neighbour at the same scale. Positive MML differing less than timepair = 0.12 s from each other are considered to be redundant and only the closest to the MML with negative polarity is kept. Additionally, MML of the same polarity differing less than timepair from the same MML of opposite polarity are also considered to be redundant. From those, the one whose amplitude normalized by the time interval to the MML with opposite polarity is at least 1.2 times higher than the other is kept, or the one closest to the MML with opposite polarity. The QRS location is taken as the zero crossing between MMLs of opposite polarity differing less than timepair. A 275 ms refractory period is included and search back performed if too long an RR interval is found. QRS onset and end are located using slope based criteria over the WT at scale 2². The above referred parameters were empirically tuned and validated over ECG signals of standard databases during the development of an automatic delineation system described in Martínez et al (2004) and set as default.

When simultaneous leads are available, multilead based global marks for main peak location are taken as the median over SL based locations for QRS candidates found in at least K of the available leads, while boundaries are taken as onset[end] of the first[last] SL mark with at least one neighbour mark in the other lead (SLR—single lead plus rules approach). In this work K = 2 was considered.

The above described method, using the default parameters of Martínez et al (2004), is applied over the direct FECG fetal QRS (FQRS). However, this same strategy applied to abdominal FECG is expected to locate only the maternal QRS complexes, as illustrated in figure 2.

Adaptations, attending to the fetal physiology and signal contamination with maternal ECG, are therefore clearly required for FQRS location over abdominal FECG. In adults, QRS complex content can range from almost 0 to 40 Hz while the frequency content of FQRS is over 20 Hz (Matonia et al 2006). This means that the WT scale 2⁴ is not useful as its equivalent frequency band does not reach 20 Hz (Li et al 1995). This fact is also illustrated in figure 2 in which no FQRS related MML are found for W_{abdominalFECG,4} signal. Thus for the adapted strategy only scales 2¹ to 2³ are used.

The lower power of the fetal contribution requires a lower threshold for fetal MML detection and a 75% reduced threshold, calculated without considering the maternal QRS intervals in its computation, is used. MML lines previously associated with maternal QRS are excluded after isolated and redundant candidates elimination, but before a polarity check of the fetal MML to define pairs. The time interval timepair used for no redundancy between MML of the same polarity and for MML pair association was reduced to 25 ms for adaptation to the shorter duration of FQRS complexes: MML associated with secondary QRS waves (like Q and small S waves) should appear closer. All numeric changes in parameters from the default were guided by the physiological differences of FECG with respect to the adult recordings, but the specific values were obtained by trial and error. Nevertheless, in the case of absence of FECG components the algorithm will detect maternal T/P waves instead, as illustrated in the lower right panel of figure 2, producing erroneous detections. Each available abdominal lead and ICA component is processed separately to produce SL and ICA based sets of locations.

For multilead based locations, a FQRS is accepted if it is detected in at least two out of the four leads (SL sets only) within a 250 ms neighborhood, as illustrated at the upper right panel of figure 2. The final SLR location is taken as the median mark. This allows us to produce a unique annotator (SLR set of locations) from which RR and HRV fetal measures could be taken.

With respect to the preliminary version presented in (Almeida et al 2013) this algorithm has incorporated slight changes. In the previous version, maternal MML was excluded at a more initial stage, before the redundancy check. Also the SLR strategy was more conservative, as it only accepted FQRS detected in at least three out of four leads.

2.3. Performance evaluation

The FQRS marks obtained that differed less than 100 ms from the reference marks were considered to be true positives (TP). Sensitivity (Se) and positive predictivity (P + ) are computed from the number of TP, false positives (FP) and false negatives (FN) as:

$$\begin{eqnarray}&&\text{Se}=\text{TP}/\left(\text{TP}+\text{FN}\right)*100\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\text{P}+=\text{TP}/\left(\text{TP}+\text{FP}\right)*100.\end{eqnarray} \tag{ 4 }$$

The FQRS location errors were computed for each record as the difference in ms from obtained marks and reference marks considering only TP beats. Additionally FHR was estimated as the inverse of the mean (or median) RR interval, and computed for each minute for both reference and SLR marks, considering all QRS marks (TP and FP). The errors in FHR estimation are taken as estimated minus reference values. The scores relative to FHR and RR series defined for the CinC_Ch as described in Clifford et al (2014), were obtained for the Challenge dataset and are also presented.

All the above described metrics were obtained for data from the Silesia dataset and Challenge dataset A, for which reference annotations were available. For Challenge dataset B only the mean scores were provided by the Challenge organizers, and thus are the only performance metrics we can present for these data.

The modified detector was applied over each of the abdominal leads (SL_j, j = 1, 2, 3, 4) and ICA (ICA_k, k = 1, 2, 3, 4) derived leads. The SLR marks were then obtained by combining the four SL based sets obtained over the abdominal leads. The original SL approach was also applied over the direct FECG data in the Silesia Data, for comparison purposes. In figure 3(a) (left panel) are plotted the P + and Se values, including the values regarding direct FECG, for the 5-min Silesia recordings of the SLR approach and considering each of the four available abdominal FECG leads (SL) and four ICA derived signals separately. In the right panel of the same figure are plotted the distributions of the FQRS location errors, for both direct, SLR, SL₄ and ICA₃, the SL and ICA leads with a higher number of TP detections. In this and similar plots the central box limits correspond to the first and third quartiles, with the median marked as a horizontal line inside the box. Values out of the 1.5 of the inter-quartile range (IQR) are marked as ' + ', as is usual in box plots. A more local analysis is presented in figure 3(b), where the Se and P + values per minute of file were plotted, considering direct FECG, SLR, SL₄ and ICA₃. Table 1 summarizes the reference FHR and respective estimation errors, considering both mean and median.

Zoom In Zoom Out Reset image size

Regarding Challenge dataset A, the obtained distributions across files of P + and Se for SLR, SL_j and ICA_k derived signals are presented in the upper panel of figure 4, with the FQRS location error distributions for SLR in the lower panel. In figure 5, the FHR estimates based on SLR locations are plotted versus the FHR based on the reference marks for each file, along with the error distributions in FHR estimation (middle panel), considering both mean and median as estimators, and score distributions (right panel). Mean|median FHR estimation errors were 2.03|0.26 bpm, respectively. The correlation between reference and estimated FHR was significant (p < 10 ^− 2) and found to be 0.30|0.45, increasing to 0.48|0.67 (p < 10 ^− 4) when considering only records for which estimated FHR > 110 bpm. This corresponds to excluding 8 out of 74 files (circles on figure 5), all with negative errors, representing FHR underestimation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering 1 min based FHR estimation in the data from both databases with reference annotations, the mean|median errors were 1.77|0.13 bpm and the correlation between reference and estimated FHR was 0.35|0.48 (p < 10 ^− 3), increasing to 0.60|0.73 (p < 10 ^− 9)when considering only estimations with FHR > 110 bpm (excluding 8 out of the 99 1 min segments).

The mean scores relative to the FHR and FRR series defined for the CinC_Ch obtained for dataset A were 386.1 and 29.7, respectively, while for dataset B were 513.1 and 35.3, respectively, exhibiting some improvements with respect to the results reported in Almeida et al (2013).

The original detection algorithm (Martínez et al 2004) applied to direct FECG is able to achieve a high detection performance, with a global error lower than 1% and a maximum location error of 20 ms. This same strategy applied to abdominal FECG was expected to be able to locate maternal QRS. Visual inspection allowed us to conclude that it performed correctly, as illustrated in figure 2 where maternal QRS complexes are clearly visible and associated MML lines are marked. Nevertheless, no systematic validation of that fact could be done due to the lack of reference annotations for maternal QRS locations in the data currently available. This is a limitation of this work and we expect that it can be solved in the future.

The P + and Se values for the proposed SLR approach were found to be better or equivalent to the best SL or ICA approach for all records in the Silesia data, as well as for most of the 1 min segments. Similar results were achieved regarding FHR estimated values computed in the Challenge data, where SLR also outperforms SL and ICA.

Data from the Silesia database have been previously used in the validation of algorithms for the automatic FQRS detection from FECG, namely in the works of Kotas et al (2011) and Castillo et al (2013). However the direct comparison between our results and those previously published is not possible, as those authors did not report which files were considered or did not use the same validation criteria.

Kotas et al (2011) proposed spatio-temporal multichannel filtering to construct a new signal, aiming to enhance fetal cardiac activity. The FQRS were detected over each of the four original abdominal FECG, a manually selected ICA component and that constructed signal, using several detection approaches. Three 5 min files from the Silesia dataset were used in the published performance evaluation, with results highly dependent on the file: P + and Se ranged from 100% to close to 85% for the proposed method, or as low as 35% for a single lead based result. As the authors did not report which files were considered, the direct comparison of our results with those previously published is not possible.

Castillo et al (2013) used a wavelet based pre-processing strategy, followed by threshold based FQRS detection; the threshold values are file specific, requiring a training dataset, as they depend on both physiological and technical factors. The files and channel evaluated there were manually selected according to the quality of FECG related activity present and reference marks were manually checked by a medical specialist who had also validated the reported results. The authors included four out of the five Silesia files, for the selected leads, reporting average Se and P + values of 98.31% and 98.22%, respectively. The lower performance of our results with respect to that work could be related to the differences in the validation criteria. Our results reported here did not involve any kind of manual selection or verification of the signals nor annotations and use all files and leads. All the processing was strictly automatic. Also, the criteria to decide if an obtained mark corresponds to a true FQRS is constant, as we discard all candidates differing by more than 100 ms from the reference mark as a FP. As a matter of fact, a 100 ms tolerance can be considered to be small, as FQRS durations above 70 ms were reported in normal fetuses (Chia et al (2005)).

The lack of annotated abdominal FECG databases was one of the main difficulties in the validation of the automatic detection methods. Most of the published approaches use their own data or annotations, which invalidates a correct comparison. The CinC_Ch is, to our knowledge, the larger database of such data. Unfortunately, reference FQRS locations are only available for set A, as the main part of the data is reserved for blind validation, using a score defined for the Challenge. Therefore we did not have access to the number of TP, FP or FN to validate our method over the set B of the Challenge data. Other open source databases including non-invasive recordings lack reference FQRS marks.

The adaptation of the original QRS detector for processing abdominal FECG used in this work allowed the algorithm to detect the FQRS by using a lower threshold for fetal MML detection. The price of a 75% reduction is a higher number of candidates for MML. This did not represent a problem since the same protections against isolated and redundant local maxima were sufficient to eliminate them in most cases.

The FHR, as measured from FECG, presents high variability with many statistical outliers; that is, values out of [Q₁ − 1.5IQR, Q₃ + 1.5IQR], where Q_i stands for the ith quartile. The clinically useful measure regarding 1-beat based FHR is not the beat-by-beat instantaneous value, but rather a representative value of the minute, that is the central tendency evaluated on the time interval. Usually the mean is used, nevertheless it is not a robust measure of central tendency, as a single large outlier can throw it off. Considering the reference annotations for set A Challenge data, more than 26% of the 1 min long files has a percentage of fetal RR intervals that are statistical outliers higher than 5%, according to the above criteria. More than 13% have a percentage of outliers higher than 10% . In the presence of outliers, very frequent in this kind of data, as seen above, the median should be used instead, as it still has a breakdown point of 50%, by definition, being a more robust central tendency measure.

SLR allowed estimation of the median FHR, with an error of less than 5 bpm for all but two 1 min subsegments out of 25 in the Silesia data. For the Challenge data the error was lower than 5 bpm for more than 50% of the processed files and lower than 20 bpm for more than 80% of the files. The median error in 1 min based FHR estimation considering both databases was 0.13 bpm, and the correlation between reference and estimated median FHR was 0.5. This allowed us to conclude that the proposed method is able to provide a clinically useful estimation of the FHR baseline. Nevertheless, the proposed method relevantly underestimated FHR for 8 out of the 99 min in which reference fetal RR can be obtained, with both estimated mean and median FHR below 110 bpm, producing a bradycardia false positive. One of these cases is illustrated in the lower right panel of figure 2, in which no FECG components are visible and maternal P/T waves are detected instead. The gestational age and type of presentation (breech or cephalic) are important factors that may have been related to the records associated with lower quality. For instance, the influence of the vernix caseosa between 28 and 32 weeks of gestation leads to a lower amplitude of the FECG. As these factors were unknown for the main records in the considered data, it was not possible in this work to confirm whether they may have been associated with records presenting lower detection performance. This is an important issue for future evaluation. Excluding those eight files allowed us to increase the correlation between reference and estimated median FHR to 0.73.

This work was focused in the correct location of the FQRS locations. The proposed wavelet based methodology does not require a specific transformation/separation method regarding the FECG analysis. The use of ICA did not improve the performance, with SLR presenting the best results. No post-processing with regard to cardiac rhythm was considered. The proposed approach seems promising for assessing fetal cardiac rhythms from abdominal ECGs, despite the fact that a high number of errors were present for some files. The performance of the detector strongly depends on the quality of the data, and thus pre-processing methods for discarding very low quality signals should be considered. The results obtained allowed us to conclude that the proposed methodology is able to provide a clinically useful estimation of the FHR.

CMUP research funded by the European Regional Development Fund through the programme COMPETE and by the Portuguese Government through the FCT—Fundação para a Ciência e a Tecnologia under the project PEst–C/MAT/UI0144/2013. CIBER-BBN is an initiative funded by the VI National R&D&i Plan 2008-2011, Iniciativa Ingenio 2010, Consolider Program, CIBER Actions and financed by the Instituto de Salud Carlos III with assistance from the European Regional Development Fund. This study was partially supported by projects TEC2010-21703-C03-02 and TEC2013-42140-R from the Ministerio de Economia y Competitividad (MINECO) with European Regional Development Fund (FEDER), Spain, and by Grupo Consolidado BSICoS from DGA (Aragón, Spain) and European Social Fund (EU). Hernâni Gonçalves is funded by a post-doctoral grant from the Fundação para a Ciência e a Tecnologia, Portugal (FCT: SFRH/BPD/69671/2010).