New Estimators and Guidelines for Better Use of Fetal Heart Rate Estimators with Doppler Ultrasound Devices

Characterizing fetal wellbeing with a Doppler ultrasound device requires computation of a score based on fetal parameters. In order to analyze the parameters derived from the fetal heart rate correctly, an accuracy of 0.25 beats per minute is needed. Simultaneously with the lowest false negative rate and the highest sensitivity, we investigated whether various Doppler techniques ensure this accuracy. We found that the accuracy was ensured if directional Doppler signals and autocorrelation estimation were used. Our best estimator provided sensitivity of 95.5%, corresponding to an improvement of 14% compared to the standard estimator.


Introduction
Continuous monitoring of fetal parameters has shown their advantages in estimating fetal wellbeing [1]. According to the report of the Society for Maternal-Fetal Medicine [2], continuous fetal heart rate monitoring has reduced infant mortality. The development or the improvement of noninvasive methods dedicated to continuous fetal monitoring is therefore of major interest.
One important parameter in assessing fetal wellbeing is the variability of the fetal heart rate. This parameter, which corresponds to the variation between intervals of two consecutive heart beats, is an indicator of central nervous system development [3][4][5]. It characterizes fetal behavior states [6][7][8] and can be an indicator of further neurological evolution [9]. Variability analysis provides a good indication of fetal distress [10] and identifies fetuses with intrauterine growth retardation [11]. According to Dawes criteria [12], variability of 4 ms is a predictor of the lack of acidosis, while a value of 2.6 ms is critical for the fetus. In the normal fetal heart rate range (110-160 bpm), a time variability of 4 ms corresponds to a cardiac frequency variability of 0.81 bpm, while a time variability of 2.6 ms corresponds to a cardiac frequency variability of 0.53 bpm (the values of 0.53 bpm and 0.81 bpm are obtained as follows: 60/(60/110−2.6 ms)−110 = 0.56 bpm and 60/(60/110 − 4.0 ms) − 110 = 0.81 bpm, resp.). Other authors [13] suggest that it is necessary to estimate the heart rate with an accuracy of 0.25 bpm in order to analyze fetal heart rate variability correctly. Reliable estimation of fetal heart rate and hence of heart rate variability is therefore essential.
Several methods are available to assess the fetal heart rate. These methods differ both at the Doppler signal level (directional or nondirectional) and at the level of the algorithm that estimates the heart rate. For example, existing devices on the market such as Sonicaid Oxford (Oxford Sonicaid Instruments, Abington, UK) [14], Hewlett-Packard 8030A (Hewlett-Packard, Palo Alto, CA, USA) [15], and Philips Avalon F40 (Philips, Amsterdam, Netherlands) use the envelope of the nondirectional Doppler signal. Other authors [16,17] have used the envelope of the directional Doppler signal. Several algorithms based on autocorrelation are commonly used to estimate the heart rate (Oxford Sonicaid, Hewlett-Packard 8030A, and Philips Avalon F40). These algorithms have been applied either directly to the Doppler signal envelope (directional or nondirectional) or to 2 Computational and Mathematical Methods in Medicine the signals resulting from discrete wavelet decomposition of the envelope. In the latter case, the final estimated heart rate involves the combination of different estimates [16].
In this study, we first verified whether the pulsed Doppler techniques currently used in commercial devices ensure such an accuracy of 0.25 bpm, and we propose here some recommendations regarding parameter settings. We also compared the techniques used in commercial devices with other pulsed Doppler techniques that use directional Doppler signals. We evaluated the efficacy of all these techniques empirically (error of estimation of the fetal heart rate, sensitivity, and false negative rate).
The originality of this study lies in the recommendations on the parameter settings of the system and in the description of the individual limitations of each technique. Finally, to improve detection probability, a new method based on the combination of heart rates obtained from directional signals is proposed.

Materials and Methods
In this section, we describe the Doppler system we developed, the patients, and the synthetic signals used for the comparison of each estimator. The synthetic signals were inferred from real signals. Finally, we present the various existing techniques for estimation of fetal heart rate and we present a new technique based on a combined procedure.

The Doppler
System. In order to evaluate fetal wellbeing objectively and to classify the fetus, we codeveloped the pulsed, multitransducer, multichannel Doppler Actifoetus unit with Althaïs Technologies (Tours, France).
Our system comprised a personal computer (PC) and our Actifoetus unit. The Actifoetus unit contained three groups of four transducers and a Doppler acquisition board. The detailed operating functions of the acquisition board were presented in [18].
The transducers exploring the fetal heart were nonfocused and monoelement. They were circular in shape, with a diameter of 13.5 mm and an acoustic power of 1 mW/cm 2 . Geometrically, the transducers were located at the center of gravity and at the top of an equilateral triangle with sides measuring 40.7 mm.
The transducers were placed on the mother's abdomen. They transmitted a sinusoidal pulse at 2.25 MHz with a pulse repetition frequency (PRF) of 1 kHz. Note that a theoretical accuracy of 60/2/1000 = 0.03 bpm can be achieved with this value of 1 kHz and accuracy can be still further improved by performing interpolation of the correlation function. The wave was propagated through the mother's abdomen towards the fetal heart. The backscattered signal was recorded from five different depths, annotated 1 , . . . , 5 . Note that only one channel was considered in the present study.
The ultrasound signal received was converted into an electrical signal and amplified to compensate for the attenuation of 1 dB/cm/MHz. The signal was then demodulated in phase ( ) and quadrature ( ) [19]. After demodulation, the signals were digitized. The digital outputs of the converters represented the digital Doppler signal.

Patients.
The Doppler signals were acquired at the CHRU "Bretonneau" Tours, France. The consent of each patient was obtained and the study was approved by the Ethics Committee of the Clinical Investigation Centre for Innovative Technology of Tours (CIC-IT 806 CHRU of Tours). Patients were older than eighteen years and all pregnancies were single. The recordings were made during the twenty-fifth and fortieth gestational weeks. Evolution during pregnancy was normal for all fetuses.

Simulation.
Because it was difficult to quantify the effectiveness of the estimation techniques directly on real signals and because there was no suitable model, we generated synthetic signals. These synthetic signals were used as a groundtruth to evaluate the effectiveness of each estimator. To make these signals as realistic as possible, we proceeded in two stages: an analyzing stage deducing the characteristics of the real Doppler signal envelope and a synthetic phase providing realistic simulated signals. Figures 1(a)   and the two envelopes of corresponding directional signals obtained from the and signals [19,20]. The synthetic envelope signal was calculated as follows: where ( ) and ( ) are the envelopes of directional Doppler signals produced by the scatters that approach and move away from the transducer, respectively.
We verified (Figures 1(a) and 1(b)) that the envelope of the real nondirectional Doppler signal had the signatures of both envelopes of the directional signals. For example, at around 500 ms, the envelope of the nondirectional signal was mainly influenced by scatters that approached the transducer, while at around 400 ms, we observed the influence of movements away from the transducer. The alternating influence of these two movements determined the envelope of the nondirectional signal. Figure 2 shows 2000 ms of the envelope of a real directional Doppler signal. In order to find the important parameters required for the synthesis of this signal, we extracted its intrinsic features (the number and amplitudes of the peaks, the lags between the peaks, and the differences in amplitude between the peaks). The values of these parameters were evaluated by considering a quasiconstant fetal heart rate. Figure 2 represents a sequence of several patterns. These patterns were made up of peaks that corresponded to cardiac wall and valve movements of the fetal heart. As suggested by Shakespeare et al. [21], although six peaks (atrial contraction, ventricular contraction, opening and closing of the mitral valves, and opening and closing of the aortic valves) could be detected theoretically, only a few peaks were in practice detected in the nondirectional Doppler signal. From our analysis, it appeared that the most likely pattern was that with four peaks. Note that this signature composed of four peaks could vary considerably from one beat to another, and it was similar to that identified by Jezewski et al. [13]. Among all these patterns, the most likely was the pattern with peaks in the order 2143; that is, the highest peak 1 was in second position, the second highest peak 2 was in first position, and so forth. For the 2143-pattern, we evaluated the amplitudes of each peak ( 1 , 2 , 3 , and 4 ), the peak durations ( 1 , 2 , 3 , and 4 ), the lags between two consecutive peaks ( 1 2 , 2 3 , and 3 4 ), and the differences in amplitude between two consecutive peaks ( 2 1 , 1 4 , and 4 3 ). The results of this statistical analysis are reported in Table 1.

Analysis of Real Directional Signals.
As the patterns observed in Figure 2 were noisy, we decided to assess the noise level in order to simulate noisy synthetic Doppler signals. We assessed the signal to noise ratio (SNR) as follows: where and are the powers of the active and passive regions, respectively. We considered the active region as the area containing the pattern peaks, whereas there were none in the passive region. Using (2), we found that the SNR calculated on our real signals corresponded to a Gaussian law: SNR ∼N (11, ( √ 2.5) 2 (dB)). x F

Synthesis of Synthetic Directional
x B x e I Q Figure 3: General diagram of the Doppler data processing acquired using one transducer and one channel. and represent the envelopes of the directional signals, and represents the envelope of the nondirectional signal. 1 , 2 are the two autocorrelation estimators.  3 4 represent the differences between the positions of two adjacent maxima over time; 2 1 , 1 4 , 4 3 represent statistical differences between two adjacent maxima over time, and represents the statistics of the peak durations. We found that the statistics of the four peaks of , = 1, . . . , 4, were identical. these characteristics. Equation (3) shows the two possible components of such a signal: where ( ) is the noise, is the peak amplitude, = 1/(2 ) is the peak frequency, Rect ( ) is the unit rectangular function centered on with width and = − , is the pattern duration, and is the synthetic signal period. We set a constant interval between the highest peaks of two consecutive patterns of the synthetic signal, as illustrated in Figure 2. We also chose the pattern period as 50% of the synthetic cardiac cycle period , since this period can vary between 40 and 60% [22].
Using (3), we generated two synthetic noisy envelopes corresponding to the envelopes of the directional signals. The envelope of the nondirectional synthetic signal was modeled using (1), being the sum of the envelopes of the both directional synthetic signals. In order to simplify our study, only ( ) was calculated, ( ) being a delayed and amplified version of ( ). To simulate realistic signals, we introduced a lag between the directional components: where ( ) and ( ) are the envelopes of directional signals and is the lag between the two envelopes. is a factor that represents the amplitude ratio between the two types of envelope. From Figure 1, ≈ 40 ms and ≈ 2.

Estimators.
In this section, we describe the different estimators used in our study. Each estimator that was based on the autocorrelation function was denoted by , = 1, . . . , 8, as illustrated in Figure 3. Each estimator operated on different signals: ( ), ( ), and ( ). Devices existing on the market currently use the ( ) envelope and autocorrelation. The estimators that used these configurations were 7 and 8 . The mathematical expression Computational and Mathematical Methods in Medicine 5 of the two autocorrelation estimators is given in [23] and is represented thereafter by where is the size of the analyzing window, is the time for which the estimator is computed, and is the lag. ( ) represents one of the signals analyzed ( ( ), ( ), or ( )).
2.4.1. Algorithm. The algorithm to estimate the fetal heart rate was the same for all three signals ( ( ), ( ), and ( )). The steps of the algorithm were as follows.
(1) Extract from each signal under consideration ( ( ), ( ), or ( )) a limited number of samples, being the window size.
(4) From the position of peaks of 1 ( ) and 2 ( ), determine the durations between consecutive peaks with = 1, 2, . . . , − 1. This conditional test also permits removal of cardiac frequency estimates that are half the expected value, as are sometimes observed (Shakespeare et al. [21]).
(6) Calculate the average cardiac frequency (FHR) from CF not exceeding 35 bpm [13]: As an illustration, consider a window of 4.096 s. Whenever the cardiac frequency was 240 bpm, 16 peaks were observed in the autocorrelation function. Using an empirically set threshold, the duration between each peak was measured Note that such an algorithm is not perfect since it is hypothesised that the FHR is constant during the process. Sometimes the second peak of the autocorrelation can be lower than the third and the FHR estimate is incorrect. A process must be performed to remove outliers.

Elimination of Outliers.
In order to eliminate outlier estimates associated with estimator dysfunction, we introduced a postprocessing step. This postprocessing step was applied only in the case of real signals. An estimate was considered to be an outlier if it laid outside the statistic computed from 40 previous estimates, or if it differed between two consecutive analysis windows by 35 bpm.

Combination.
In order to improve the effectiveness of the FHR estimation, we combined estimations. For 2 and 5 , the two values of the fetal heart rate estimated on signals ( ) and ( ) were combined. The estimate on the two signals was achieved using 1 or 2 . The combination rule we used was as follows: (i) if the heart rate was detected on a single signal, the combined value took this value; (ii) if the heart rate was detected on both signals, the combined value was the average of the two values.
Note that, in contrast to Kret's study [16], we combined the fetal heart rates estimated on both directional Doppler signals, while Kret's technique was based on combination of fetal heart rate estimations computed after discrete wavelet decomposition of the envelope of the directional Doppler signal. Since Kret's technique was applied only for continuous Doppler signals, it was not taken into consideration in our study.

Results
To find the best estimators and the conditions in which they could be used, we performed a series of simulations and experiments. Using simulations, we sought configurations that ensured an error of estimation, that is, the expected accuracy below 0.25 bpm, the highest sensitivity, and the lowest false negative rate. Experimentally, we sought the best configuration for optimal use of the estimator.

Simulated Signals.
We present the results from two types of simulation. In the first series of simulations, we sought parameter settings that ensured the desired accuracy of 0.25 bpm. In the second series of simulations, we evaluated the effectiveness of each estimator in terms of true positive rate and false negative rate. . We varied the signal periodicity between 1000 and 250 ms, as these values corresponded to the standard range of exploration (60-240 bpm) of different fetal monitors. The SNR range varied between 0 and 14 dB, in order to include our measured SNR values on the real signals and in order to take into account the worst cases. The range of size varied between 512 and 4096 ms. The highest fetal heart rate could be obtained with a window size of 512 ms, although we limited the maximum window to 4096 ms to reduce computation time. Figures 4 and 5 show the smallest window size analyzed ( = 4096 ms) of all estimators tested that ensured the expected accuracy of 0.25 bpm in the range of 60-240 bpm and that ensured a SNR at least greater than 0.6 dB. Note that for estimator 2 reported in Figure 5, we showed that there was no size which ensured the desired accuracy, whatever the SNR or the frequency. To test estimator robustness in relation to the increasing complexity of the simulated signals, the lag varied between 0 and 40 ms, this value of 40 ms being taken from Figure 1.
The results derived from Figure 4 showed that accuracy for the envelope signal (estimator 7 ) was no longer achieved   Figure 5: Duration of the analyzing window W required to reach an error of at least 0.25 bpm with the autocorrelation estimator ( 2 ) and with a SNR > 0.6 dB.
for certain frequencies, but it still was for directional signals. Finally, Figure 5, shows the best estimators ( 1 , 2 , 4 ) and their respective best parameter settings ( = 4096) that ensured an accuracy of 0.25 bpm with a SNR > 0.6 dB in the 60-240 bpm range.
To summarize, these first results showed the superiority of 1 compared to 2 and the superiority of the envelope of directional signals compared to that of nondirectional signals. We therefore recommend the use of 1 and the estimators ( 1 , 2 , and 4 ) based on the envelope of directional signals.

Performance Levels of Estimators.
In this study, the performance levels of estimators we wanted to compute were sensitivity and the false negative rate. Fetal heart rates were evaluated every 250 ms from noisy signals with = 4096 ms. Sensitivity was computed with the equation: = TP/(TP + FN), where TP was the true positive rate and FN was the false negative rate. Estimates of simulated heart rate were considered to be false negative if they did not ensure the expected accuracy; otherwise, they were true positive. Sensitivity and the false positive rate were evaluated as the True positive rate (%) for estimator I 1 True positive rate (%) average of 30 values. Each value was determined after analysis of a noisy signal of 30 s where sensitivity and the false positive rate had converged to the highest value and to the lowest value, respectively. Convergence was reached for a minimum SNR of 6 dB. The results of this second series of simulations are presented in Figures 6 and 7 for = 4096 ms and SNR ≥ 6 dB. Note that the estimation of error of 25 bpm was obtained only for 1 whatever the frequency, whereas accuracy for 2 was obtained only for 100, 150, 200, 220, and 240 bpm.
The results set out in Figure 7 show that estimators based on 1 ( 1 , 2 , 4 , and 7 ) had an average (average obtained from the cardiac frequency) false negative rate of 1.5%, while those based on 2 ( 3 , 5 , 6 , and 8 ) presented a higher average false negative rate of approximately 14.8%. The 97.5% average true positive rate of 1 was slightly lower than that of estimators based on 2 , which was 100%. Finally, when the accuracy of 0.25 bpm was reached, we observed that the estimators based on 1 were generally more accurate than those based on 2 , although the average false negative rate was not zero. Figure 8 shows the error of estimation corresponding to different values of SNR when a zero false negative rate was imposed. This zero false negative rate was obtained by modifying the detection threshold, and a direct consequence was an increase in the estimation. The results derived from Figure 8 showed that the zero false negative rate for 1 was ensured for a SNR ≥ 2dB (below the SNR measured on real signals) and for an error of estimation of 0.8 bpm. In the case of 2 and directional signals, we obtained an error of estimation of 4 bpm, whereas for a nondirectional signal it was 6 bpm. Table 2 summarizes the effectiveness of each estimator in terms of FHR error of estimation, SNR, and average false negative rate. To reach an error of estimation, that is, the expected accuracy of 0.25 bpm, we recommend 1 , that is, autocorrelation-based estimators ( 1 , 2 , 4 , and 7 ), the price to be paid being an average false negative rate of 1.5%. To reach an average false negative rate of 0%, we recommend 1 autocorrelation-based estimators ( 1 , 2 , 4 , and 7 ), where the price to be paid is an error of estimation of 0.8 bpm far from the expected accuracy of 0.25 bpm.

Results
Obtained on Real Signals. We recorded 580 minutes for the analysis of real Doppler signals. We selected areas where signals had the cardiac activity signature. The performance levels on these signals were evaluated on the envelopes of both the nondirectional and the directional signals. The FHR estimation obtained with a commercial device (Oxford Sonicaid) was used as a reference to evaluate sensitivity which was evaluated for each estimator. All estimators were evaluated using a size of = 4096 ms.
The results obtained for all the signals are presented in Table 3. The estimators based on directional signals ( 1 , 4 , 3 , and 6 ) provided a higher level of sensitivity compared to those which used nondirectional signals ( 7 , 8 ). In the case of directional signals, the results of 1 and 2 were close but slightly better for 1 . This result confirmed the results obtained by simulations. We therefore recommend the use of estimators based on 1 calculated on the directional signals ( 1 , 4 ).

Discussion and Conclusion
In this study, we focused on different settings of the estimators (window size, lag) that ensured a fetal heart rate estimation with a maximum authorized error of 0.25 bpm. We found in simulation that only estimators based on 1 and directional signals could ensure such an accuracy of 0.25 bpm. The size necessary in this case was = 4096 ms. Note that, although we proposed synthetic signals that were as realistic as possible, we are aware that the plotted performance levels are representative only of our simulations and not of all cases encountered in practice. It is likely that the performance levels of the algorithms tested can be reduced in the presence of artefacts. However, the 95% sensitivity obtained from real signals suggests that our proposed estimators may be trusted.
In the case of real signals, the sensitivity was quantified. Since in our study the estimated SNR on the real signals was greater than the threshold of 6 dB (deduced on simulated signals) required to reach the desired accuracy of 0.25 bpm, a denoising filter was not necessary. However, in cases of a SNR lower than 6 dB, a denoising process (Wiener, wavelet) could be introduced.
Sensitivity was quantified using a size of 4096 ms. We found that the estimators 1 , 2 , and 4 based on 1 had slightly greater sensitivity than those based on 2 . We therefore recommend the use of 1 .
Various cases were considered on the basis of this study, that is, those that do not require a precise estimate of the fetal heart rate and those for which accuracy is critical. The accuracy of fetal heart rate estimation in the first case is not important but the false negative rate should be as low as possible. For example, if the goal of a monitoring system is simply to verify that the fetal heart rate is in the normal range (110-160 bpm), very high accuracy is not needed. In this case, an error of estimation of 0.8 bpm is sufficient. Our computations showed that for an error of estimation of 0.8 bpm, a zero false negative rate in the zones when the rhythm is quasi-constant could be ensured. In the second case, an error of estimation of 0.25 bpm is required for a system in which the goal is not only to estimate the heart rate, but also to evaluate fetal wellbeing. Our study showed that for this type of system, the false negative rate may be slightly higher than zero. It is important to note that this error of estimation was guaranteed for a quasi-constant heart rate. This is not a constraint for such a system, since the variability of fetal heart rate must be evaluated in these ranges to predict fetal distress.
Applied to real signals, the estimators based on 1 provided sensitivity close to those of 2 , and the most efficient of these estimators were those that used directional signals ( 1 , 2 , and 4 ).
A 7% increase in sensitivity compared to estimators based on individual directional signals was possible when we combined the two heart rates calculated on the directional signals. A 14% increase in sensitivity compared to estimators based on individual nondirectional signals was possible when we combined the two heart rates calculated on the directional signals. When a combination was used, both signals were processed in parallel, thus doubling the number of operations.
The good levels of performance of our estimator based on this combination suggest first that it can be adapted to multitransducer, multichannel configurations and second that such an estimator will improve fetal diagnosis.