New method for beat-to-beat fetal heart rate measurement using Doppler ultrasound signal

sensors Article New Method for Beat-to-Beat Fetal Heart Rate Measurement Using Doppler Ultrasound Signal Tomasz Kupka 1,*, Adam Matonia 1, Michal Jezewski 2, Janusz Jezewski 1, Krzysztof Horoba 1, Janusz Wrobel 1, Robert Czabanski 2and Radek Martinek 3 1 Łukasiewicz Research Network—Institute of Medical Technology and Equipment, PL41800 Zabrze, Poland; [email protected] (A.M.); [email protected] (J.J.); [email protected] (K.H.); [email protected] (J.W.) 2Department of Cybernetics, Nanotechnology and Data Processing, Silesian University of Technology, PL44100 Gliwice, Poland; [email protected] (M.J.); [email protected] (R.C.) 3Department of Cybernetics and Biomedical Engineering, VSB—Technical University of Ostrava, 70800 Ostrava-Poruba, Czech Republic; [email protected] *Correspondence: [email protected] Received: 27 May 2020; Accepted: 20 July 2020; Published: 22 July 2020   Abstract: The most commonly used method of fetal monitoring is based on heart activity analysis. Computer-aided fetal monitoring system enables extraction of clinically important information hidden for visual interpretation—the instantaneous fetal heart rate (FHR) variability. Today’s fetal monitors are based on monitoring of mechanical activity of the fetal heart by means of Doppler ultrasound technique. The FHR is determined using autocorrelation methods, and thus it has a form of evenly spaced—every 250 ms—instantaneous measurements, where some of which are incorrect or duplicate. The parameters describing a beat-to-beat FHR variability calculated from such a signal show significant errors. The aim of our research was to develop new analysis methods that will both improve an accuracy of the FHR determination and provide FHR representation as time series of events. The study was carried out on simultaneously recorded (during labor) Doppler ultrasound signal and the reference direct fetal electrocardiogram Two subranges of Doppler bandwidths were separated to describe heart wall movements and valve motions. After reduction of signal complexity by determining the Doppler ultrasound envelope, the signal was analyzed to determine the FHR. The autocorrelation method supported by a trapezoidal prediction function was used. In the final stage, two different methods were developed to provide signal representation as time series of events: the first using correction of duplicate measurements and the second based on segmentation of instantaneous periodicity measurements. Thus, it ensured the mean heart interval measurement error of only 1.35 ms. In a case of beat-to-beat variability assessment the errors ranged from − 1.9% to − 10.1%. Comparing the obtained values to other published results clearly confirms that the new methods provides a higher accuracy of an interval measurement and a better reliability of the FHR variability estimation. Keywords: fetal monitoring; Doppler ultrasound signal; fetal heart rate; beat-to-beat variability 1. Introduction Cardiotocography is a crucial part of modern perinatal medicine, which enables to monitor and assess a fetal condition during pregnancy and labor. It involves recording a fetal heart rate (FHR) signal against uterine contractile activity and fetal movements [ 1 – 5 ]. Normal fetal heart activity is an indirect evidence of adequate fetal oxygenation and preservation of the central nervous system functions [ 6 ]. The FHR signal values in the range 110 ÷ 150 bpm confirm a normal fetal heart activity. Accelerations Sensors 2020,20, 4079; doi:10.3390/s20154079 www.mdpi.com/journal/sensors Sensors 2020,20, 4079 2 of 25 or decelerations are defined as episodes of increasing or slowing of the heart rate. Another activity feature is an instantaneous variability, represented by two components: short-term variability (STV) describing FHR fluctuations occurring from beat-to-beat and long-term variability (LTV) defining periodic fluctuations in direction and extent of the STV. The highest accuracy of the FHR signal determination is ensured by recording primary electrical excitations of a fetal heart, i.e., the fetal electrocardiogram (FECG) [ 7 , 8 ], which ensures very precise locations of characteristic points—R waves, representing subsequent fetal heartbeats. Based on time T i between two consecutive R waves, the instantaneous FHR i value is determined in beats per minute [bpm]: FHRi=60, 000 Ti (1) where: Tiis the interval between two consecutive heartbeats expressed in milliseconds. Unfortunately, the use of fetal electrocardiography during pregnancy faces serious problems [ 9 – 12 ]. When acquiring an indirect FECG from maternal abdominal surface, its amplitude is much lower than the electrical activity from mother, and in addition, a high level of other interferences often prevents successful T i measurement [ 13 – 15 ]. A good quality signal can be obtained invasively from an electrode directly attached to a fetal head, which can only be carried out during an advanced labor. These limitations caused that already in early 70s a Doppler ultrasound (US) method [ 16 – 22 ] was applied for monitoring of mechanical activity of the fetal heart. The Doppler echo signal reflected from moving heart is characterized by a complex structure originating from opening and closing of valves, as well as from movements of heart walls [ 23 – 28 ]. In addition, it is characterized by high variability in time due to fetal and maternal movements [ 3 , 29 , 30 ] and changing of transducer position in relation to a moving signal source—the fetal heart [ 31 – 33 ]. Identification of the heartbeat occurrence, based on a simple peak detection in looking for characteristic component, is very inaccurate due to its variable shape [ 23 , 34 ]. In addition, any interfering impulses of sufficiently high amplitude can be identified as successive beats [ 23 , 35 – 37 ]. All of these makes practically impossible to indicate points equivalent to the R-waves in FECG [ 38 ]. Therefore, for FHR determination based on the US signal, correlation methods considering the full impulse shape are used [33,35,39–42]. The dominant method is the autocorrelation function (AF), which comparing to the cross correlation, is easier to perform and does not require selecting a heartbeat pattern that makes it resistant to signal interferences [ 43 – 45 ]. Unfortunately, only the periodicity is estimated, not the exact locations of subsequent heartbeats in the analyzed signal. Since measurement of the T i is repeated in windows covering at least two complete cardiac cycles, a single heart interval can be represented by more than one measurement. Such signal form makes it impossible to reliably assess a beat-to-beat variability. The physiological FHR range is from 60 to 240 bpm (which corresponds to T i intervals from 1000 to 250 ms). Therefore, not to miss any cardiac cycle, subsequent FHR values must be determined at least every 250 ms. Such repetition time is a widely accepted standard at which Doppler ultrasound monitor provides a new FHR value at its output. To date, the possibilities of parameterization and objectification of the FHR variability evaluation relate to accurate quantitative signal analysis carried out by dedicated computer-aided fetal monitoring systems [ 46 – 51 ]. Automated signal analysis ensures repeatability of interpretation [ 52 – 56 ], as well as the extraction of information hidden for visual assessment, such as instantaneous FHR variability at a beat-to-beat level [ 57 ]. At turn of the 60 s and 70 s, several different indices were defined for quantitative description of the instantaneous FHR variability [ 58 ]. Most of these definitions assume that the analyzed signal comes from an invasive direct FECG, i.e., the indices are calculated based on T i intervals determined between the R waves. In other words, the analyzed FHR signal has a form of time series of events, in which each cardiac cycle is represented by a single measurement only. Nevertheless, a direct application of the definitions originating from FHR acquired via FECG, to the FHR determined from Doppler US signal, does not give good results. Especially the short-term variability indices, calculated from measurements occurring every 250 ms, show significant errors [59–62]. Sensors 2020,20, 4079 3 of 25 Therefore, attempts have been made to process the FHR signal from a time series of evenly spaced measurements into a time series of events representing subsequent heartbeats [ 63 – 65 ]. Despite the improved method for determination of FHR variability indices [ 65 ], their accuracy was still limited due to a wide autocorrelation window resulting in averaging of several neighboring intervals. The window width selection is a compromise between noise immunity and measurement accuracy. In practice, the AF windows of a few seconds are used, which ensures the continuity of periodicity measurement even in the presence of interferences from fetal movements, however at a cost of the measurement accuracy. The determined interval value does not then represent a single cardiac cycle, but an averaged value of few last beats [ 36 , 62 , 64 , 66 ]. It does not affect the results of classical FHR analysis—the baseline estimation or detecting accelerations and decelerations, because such analysis is carried out on a signal averaged in 2.5 s periods [ 62 ]. However, averaging of measurements practically prevents a reliable STV determination regarding beat-to-beat changes [62,65]. It seems that analysis of the US signal inside a fetal monitor do not use its full potential. This results partly from a lack of requirements for this type instrumentation, due to relatively short history of fetal monitoring systems enabling automated FHR analysis with beat-to-beat accuracy. The question arises whether, by improving the built-in algorithm for determining the periodicity in a raw Doppler ultrasound signal, it is possible to obtain measurements close to reference values from the FECG. Peters [ 66 ] developed an algorithm for determining the FHR from US signal, providing it in a form of time sequence of events, which enabled correct determination of indices describing the FHR instantaneous variability [ 66 ]. The sequence of events was obtained by using markers that approximately corresponded to subsequent heartbeats. They were used to select the signal fragments for which the AF was determined. Peters showed a high correlation with the FHR signal determined from direct FECG. However, he examined only one short signal and did not evaluate an interval measurement accuracy. Such evaluation was made in [ 67 ], where an autocorrelation window comprising two cardiac cycles was positioned in US signal based on the simultaneously recorded FECG. The research concerned different values of window shift in relation to its initial position, determined by an occurrence of the R wave. Unfortunately, the work is only theoretical, because the method, as requiring an additional signal source, is not suitable for any practical use in fetal monitor. Al-Angari has also investigated the problem of determining the FHR signal from the Doppler ultrasound, as time series of events [ 68 ]. However, the obtained FHR signals were characterized by unusually high variability. Moreover, the source of reference FHR was a FECG signal, but recorded indirectly from a maternal abdomen. Jezewski also dealt with the analysis of a raw US signal [ 36 ] recorded together with a reference direct FECG signal. The width and step of the AF window were adapted to the previously determined interval value. The efficiency of a new method was estimated by the interval measurement and de Haan’s variability index determination errors [ 58 ]. Thus, the obtained results confirmed the high efficiency, but only in relation to a single FHR variability index [64]. The cited works do not allow us to answer the previously asked question about the possibility of determining the FHR from a raw US signal with the accuracy similar to the reference direct FECG. Therefore, the authors of this work conducted comprehensive research on the Doppler ultrasound signal describing the mechanical activity of a fetal heart. The aim of this research was to develop new methods that will both improve the accuracy of FHR determination and provide the resulting signal in form of time series of events. Achieving of these two goals should result in a more accurate description of instantaneous FHR variability, and thus a more effective early detection of fetal distress. The goal of our research was to develop a complete channel for processing the Doppler ultrasound signal enabling precise periodicity measurement. First, the information content of the US signal was analyzed-two subranges of Doppler bandwidths were separated using bandpass filtering to describe slower heart wall movements and valve motions of a higher speed. In a next stage, the US envelope was determined, and a classic low-pass filtering and Hilbert transform were considered. To determine the FHR signal, the autocorrelation function supported by the trapezoidal prediction function was Sensors 2020,20, 4079 4 of 25 used. In final stage, we proposed two different methods providing the FHR signal in a form of time series of events: first using correction of duplicate measurements and second based on segmentation of instantaneous periodicity measurements. In addition, very important was the original methodology for evaluating efficiency of the proposed algorithms, which was based on automated synchronization of examined FHR signals. All these aspects are important because there are no other publications comprehensively evaluating the FHR signal determination accuracy, both by estimating cardiac interval Ti accuracy and calculating clinically crucial indices describing the instantaneous FHR variability and additionally using direct fetal electrocardiography as reference data. 2. Methodology Using an ultrasound method, determination of a fetal heart rate is based on analyzing the Doppler effect, i.e., a frequency change of US beam reflected from moving parts of a fetal heart. The received returning echo is a signal with a frequency of the US wave modulated by the speed of a given moving structure. After amplification it is demodulated for separation of particular Doppler components. Then the signal is converted into a digital form, and a further processing is carried out—leading to the FHR determination. When all phases of a cardiac cycle, originating both from valve and wall movements are present (Figure 1) much more complex signal structure (US1, US2) in relation to a simultaneously recorded FECG can be noticed. The amplitude and shape of individual components change during long-term monitoring, as the object (fetus) can change its position relative to a transducer [ 31 ]. The FECG additionally presented in Figure 1, allows for detailed identification of cardiac phases observed in signals of mechanical heart activity. Sensors 2020, 20, x FOR PEER REVIEW 4 of 25 original methodology for evaluating efficiency of the proposed algorithms, which was based on automated synchronization of examined FHR signals. All these aspects are important because there are no other publications comprehensively evaluating the FHR signal determination accuracy, both by estimating cardiac interval Ti accuracy and calculating clinically crucial indices describing the instantaneous FHR variability and additionally using direct fetal electrocardiography as reference data. 2. Methodology Using an ultrasound method, determination of a fetal heart rate is based on analyzing the Doppler effect, i.e., a frequency change of US beam reflected from moving parts of a fetal heart. The received returning echo is a signal with a frequency of the US wave modulated by the speed of a given moving structure. After amplification it is demodulated for separation of particular Doppler components. Then the signal is converted into a digital form, and a further processing is carried out— leading to the FHR determination. When all phases of a cardiac cycle, originating both from valve and wall movements are present (Figure 1) much more complex signal structure (US1, US2) in relation to a simultaneously recorded FECG can be noticed. The amplitude and shape of individual components change during long-term monitoring, as the object (fetus) can change its position relative to a transducer [31]. The FECG additionally presented in Figure 1, allows for detailed identification of cardiac phases observed in signals of mechanical heart activity. In addition to an occurrence rate, particular cardiac phases are characterized by different bandwidths. The lower frequencies correspond to the heart wall, while the valve movements are represented by higher frequency components (Figure 2). The US signal envelope describing the valves is characterized by higher dynamics, and narrow peaks enabling accurate cardiac cycle length estimation. In turn, the envelope signal from the walls is much more smoothed, with wider peaks which reduces periodicity measurements accuracy. Additionally, a risk of erroneously doubling the heart rate increases, since the signal envelope fragments during the systole and diastole phases are very similar. Figure 1. An example of the one-second direct electrocardiogram (FECG) with Doppler ultrasound signals from two transducers (US1, US2) attached to maternal abdomen in different position relative to a fetal heart. In the US signal the cardiac phases are observed: opening and closing the crescent— Ao and Ac, atrioventricular valve—Mo and Mc, atrial wall contraction—Atc and walls heart chambers—Vc. Figure 1. An example of the one-second direct electrocardiogram (FECG) with Doppler ultrasound signals from two transducers (US1, US2) attached to maternal abdomen in different position relative to a fetal heart. In the US signal the cardiac phases are observed: opening and closing the crescent—Ao and Ac, atrioventricular valve—Mo and Mc, atrial wall contraction—Atc and walls heart chambers—Vc. In addition to an occurrence rate, particular cardiac phases are characterized by different bandwidths. The lower frequencies correspond to the heart wall, while the valve movements are represented by higher frequency components (Figure 2). The US signal envelope describing the valves is characterized by higher dynamics, and narrow peaks enabling accurate cardiac cycle length estimation. In turn, the envelope signal from the walls is much more smoothed, with wider peaks which reduces periodicity measurements accuracy. Additionally, a risk of erroneously doubling the Sensors 2020,20, 4079 5 of 25 heart rate increases, since the signal envelope fragments during the systole and diastole phases are very similar. Sensors 2020, 20, x FOR PEER REVIEW 5 of 25 Figure 2. Doppler ultrasound signal (US) with separated components relating to the fetal heart wall (USW) and valve (USV) movements. Corresponding envelopes (ENVW and ENVV) are shown together with an envelope of full signal (ENV). 2.1. Data Acquisition The research material used for this study was collected at the Obstetric Clinic of the Medical University of Silesia in Zabrze. It was carried out by investigators from the Medical University of Silesia during development the project “Application of direct electrocardiography for intrapartum assessment of fetal distress “(Bioethics Committee approval number: NN-013-345/02)”. The fetal monitor used in this study enables to acquire a raw Doppler US signal (1 MHz) obtained after the demodulation process. The FHR_CTG measurements calculated by built-in fetal monitor software were accessible on a digital output. Since the additional direct FECG recording (for determining FHR_E) required measuring electrode placed on fetal head, all signals were recorded as part of a routine procedure for labor monitoring. Six sessions (39 ÷ 41 week of pregnancy) lasting from 280 to 1263 s were carried out. The final research material included 58 min of simultaneously acquired: Doppler US signals being the input data for the algorithms investigated, original FHR_CTG signals determined by a fetal monitor, as well as reference information in the conducted research—FHR_E signals. The resolution of both FHR_CTG and FHR_E signals were 0.25 bpm. 2.2. Signal Processing The goal of the research was to develop a complete channel for processing the Doppler ultrasound signal enabling precise periodicity measurement. The basic assumption was that the resulting FHR signal should be in form of time series of events, i.e., each subsequent interval between two heartbeats should be represented only by a single measurement—like in FECG. The first stage of the research focused on optimal representation of the Doppler ultrasound signal (Figure 3). The signal information content was analyzed for the influence of individual cardiac phases on the accuracy of measurement. Two bandwidths of Doppler frequencies were separated using bandpass filtering: one describing wall and one relating to valve movements. The next stage concerned reduction of signal complexity by calculating its envelope. Classic low-pass filtering and the Hilbert transform were considered. For methods evaluation, the FHR values were determined by AF, calculated every 250 ms. Figure 2. Doppler ultrasound signal (US) with separated components relating to the fetal heart wall (US W ) and valve (US V ) movements. Corresponding envelopes (ENV W and ENV V ) are shown together with an envelope of full signal (ENV). 2.1. Data Acquisition The research material used for this study was collected at the Obstetric Clinic of the Medical University of Silesia in Zabrze. It was carried out by investigators from the Medical University of Silesia during development the project “Application of direct electrocardiography for intrapartum assessment of fetal distress “(Bioethics Committee approval number: NN-013-345/02)”. The fetal monitor used in this study enables to acquire a raw Doppler US signal (1 MHz) obtained after the demodulation process. The FHR_CTG measurements calculated by built-in fetal monitor software were accessible on a digital output. Since the additional direct FECG recording (for determining FHR_E) required measuring electrode placed on fetal head, all signals were recorded as part of a routine procedure for labor monitoring. Six sessions (39 ÷ 41 week of pregnancy) lasting from 280 to 1263 s were carried out. The final research material included 58 min of simultaneously acquired: Doppler US signals being the input data for the algorithms investigated, original FHR_CTG signals determined by a fetal monitor, as well as reference information in the conducted research—FHR_E signals. The resolution of both FHR_CTG and FHR_E signals were 0.25 bpm. 2.2. Signal Processing The goal of the research was to develop a complete channel for processing the Doppler ultrasound signal enabling precise periodicity measurement. The basic assumption was that the resulting FHR signal should be in form of time series of events, i.e., each subsequent interval between two heartbeats should be represented only by a single measurement—like in FECG. The first stage of the research focused on optimal representation of the Doppler ultrasound signal (Figure 3). The signal information content was analyzed for the influence of individual cardiac phases on the accuracy of measurement. Two bandwidths of Doppler frequencies were separated using bandpass filtering: one describing wall and one relating to valve movements. The next stage concerned reduction of signal complexity by calculating its envelope. Classic low-pass filtering and Sensors 2020,20, 4079 6 of 25 the Hilbert transform were considered. For methods evaluation, the FHR values were determined by AF, calculated every 250 ms. Sensors 2020, 20, x FOR PEER REVIEW 6 of 25 During the next stage, we focused on tuning algorithms for determining the FHR signal: selection of autocorrelation windows parameters both for the autocorrelation and prediction function. Next, two different methods were developed to provide FHR signal in form of time series of events: using correction of duplicate measurements [65] and basing on segmentation of instantaneous periodicity measurements. During the algorithm development, different factors were recognized as having a major impact on the final Ti measurements. No less important was the final measurements validation, if the determined interval value was correct, due to interferences originating from maternal heart activity or fetal movements. The final stage of the study concerned metrological assessment of the proposed methods in relation to reference values provided by a fetal electrocardiogram. To validate the results, a precise synchronization of FHR signals was very important to ensure comparing the corresponding intervals. Two different procedures were proposed to estimate the inconsistency. A direct comparison of signals by the differences of instantaneous FHR values was accompanied by signal loss analysis and evaluating the signals impact on the indices describing long-term and short-term FHR variability. Figure 3. Block diagram generally describing the procedure of subsequent signal processing stages and methods for final evaluation of the results obtained. 2.2.1. Signal Preprocessing In the first step, the appropriate signal bandwidth was selected to obtain “optimal” signal information content for more accurate interval measurement. Initially, it was assumed that the entire bandwidth, corresponding to both heart valve and wall movements, should be analyzed. It is difficult to measure signal periodicity precisely because each cardiac cycle in the envelope signal is represented by several impulses related to different cardiac phases. The number of events observed in a signal can be reduced by separation of signal components corresponding to different motion speeds of heart structures. In the case of fetal monitoring, it is assumed that using an ultrasound transducer emitting with 1 MHz, wall and valves movements are within bandwidth of 100 ÷ 600 Hz [36] Using band-pass filtering the subranges were separated, covering frequencies from 100 to 300 Hz (wall movements—Figure 4B) and from 300 to 600 Hz (valve movements—Figure 4C). Filtering below 100 Hz enabled to remove components originating mainly from fetal trunk and limbs movements, while above 600 Hz from flows in maternal blood vessels. Figure 3. Block diagram generally describing the procedure of subsequent signal processing stages and methods for final evaluation of the results obtained. During the next stage, we focused on tuning algorithms for determining the FHR signal: selection of autocorrelation windows parameters both for the autocorrelation and prediction function. Next, two different methods were developed to provide FHR signal in form of time series of events: using correction of duplicate measurements [ 65 ] and basing on segmentation of instantaneous periodicity measurements. During the algorithm development, different factors were recognized as having a major impact on the final Ti measurements. No less important was the final measurements validation, if the determined interval value was correct, due to interferences originating from maternal heart activity or fetal movements. The final stage of the study concerned metrological assessment of the proposed methods in relation to reference values provided by a fetal electrocardiogram. To validate the results, a precise synchronization of FHR signals was very important to ensure comparing the corresponding intervals. Two different procedures were proposed to estimate the inconsistency. A direct comparison of signals by the differences of instantaneous FHR values was accompanied by signal loss analysis and evaluating the signals impact on the indices describing long-term and short-term FHR variability. 2.2.1. Signal Preprocessing In the first step, the appropriate signal bandwidth was selected to obtain “optimal” signal information content for more accurate interval measurement. Initially, it was assumed that the entire bandwidth, corresponding to both heart valve and wall movements, should be analyzed. It is difficult to measure signal periodicity precisely because each cardiac cycle in the envelope signal is represented by several impulses related to different cardiac phases. The number of events observed in a signal can be reduced by separation of signal components corresponding to different motion speeds of heart structures. In the case of fetal monitoring, it is assumed that using an ultrasound transducer emitting with 1 MHz, wall and valves movements are within bandwidth of 100 ÷ 600 Hz [ 36 ] Using band-pass filtering the subranges were separated, covering frequencies from 100 to 300 Hz (wall movements—Figure 4B) and from 300 to 600 Hz (valve movements—Figure 4C). Filtering below 100 Hz enabled to remove components originating mainly from fetal trunk and limbs movements, while above 600 Hz from flows in maternal blood vessels. The next step involved determining the Doppler US signal envelope aimed to reduce the signal complexity. The easiest way to obtain an envelope is to convert a bipolar signal into a unipolar one using a module function and next a low-pass filtering. Too low filter cutofffrequency causes blurring of phases observed in a cardiac cycle, which, leads to lowering the AF peak and reducing measurement accuracy. On the other hand, too high cutofffrequency (e.g., 150 Hz) incorrectly shifts the peak due to high AF corrugation. The filter parameters should be adjusted to Doppler frequencies, which are Sensors 2020,20, 4079 7 of 25 proportional to the 1 MHz. During the tests, the envelope signal was determined using different cutoff frequencies (25, 50, 75, 100, 150 Hz), as well as without any additional filtering. Sensors 2020, 20, x FOR PEER REVIEW 7 of 25 Figure 4. Doppler ultrasound signal processing. (A) original signal after demodulation; (B) signal after band-pass filtering 100 ÷ 300 Hz (wall movements) and two envelopes determined using 50 Hz low-pass filter and the Hilbert transform; (C) signal after band-pass filtering 300 ÷ 600 Hz (valve movements) and two envelopes as in B. The next step involved determining the Doppler US signal envelope aimed to reduce the signal complexity. The easiest way to obtain an envelope is to convert a bipolar signal into a unipolar one using a module function and next a low-pass filtering. Too low filter cutoff frequency causes blurring of phases observed in a cardiac cycle, which, leads to lowering the AF peak and reducing measurement accuracy. On the other hand, too high cutoff frequency (e.g., 150 Hz) incorrectly shifts the peak due to high AF corrugation. The filter parameters should be adjusted to Doppler frequencies, which are proportional to the 1 MHz. During the tests, the envelope signal was determined using different cutoff frequencies (25, 50, 75, 100, 150 Hz), as well as without any additional filtering. The Hilbert transform was used as an alternative method for determining the envelope [33]. The envelope determined using the Hilbert transform was characterized by quite high variability, especially in fragments of high dynamics. Therefore, a moving average filter (MA) was additionally used to smooth the envelope. The best smoothing while maintaining the impulse shape, is ensured by averaging over 11 or 21 signal samples (established empirically). For further studies, the envelope was determined using the Hilbert transform without additional filtering and MA with a window of 11 or 21 samples. Figure 4B,C shows exemplary envelopes determined with a low-pass filter with a cutoff frequency of 50 Hz and using the Hilbert transform. To assess the impact of the bandwidth applied on the accuracy of measurements at this stage of study, a simple AF determined every 250 ms within windows of 1.5 s was used. Because at this stage of analysis we just want to compare two methods for envelope determining, so such approach, which is a compromise between ensuring the continuity of measurements and their accuracy, is the most right. Additionally, the AF threshold relating to low level of signal quality and thus defining the FHR signal loss was empirically determined. For this purpose, the most disturbed signal fragments were analyzed, in which in the absence of fetal heart episodes, many local AF peaks were observed, whose amplitude usually did not exceed 0.1. Therefore, if the amplitude of AF peak did not exceed PTH = 0.1, a loss marker was set in the resulting FHR signal. Such AF threshold is more rigid than 0.15 proposed in the study [36]. 2.2.2. Periodicity Determination Figure 4. Doppler ultrasound signal processing. ( A ) original signal after demodulation; ( B ) signal after band-pass filtering 100 ÷ 300 Hz (wall movements) and two envelopes determined using 50 Hz low-pass filter and the Hilbert transform; ( C ) signal after band-pass filtering 300 ÷ 600 Hz (valve movements) and two envelopes as in B. The Hilbert transform was used as an alternative method for determining the envelope [ 33 ]. The envelope determined using the Hilbert transform was characterized by quite high variability, especially in fragments of high dynamics. Therefore, a moving average filter (MA) was additionally used to smooth the envelope. The best smoothing while maintaining the impulse shape, is ensured by averaging over 11 or 21 signal samples (established empirically). For further studies, the envelope was determined using the Hilbert transform without additional filtering and MA with a window of 11 or 21 samples. Figure 4B,C shows exemplary envelopes determined with a low-pass filter with a cutoff frequency of 50 Hz and using the Hilbert transform. To assess the impact of the bandwidth applied on the accuracy of measurements at this stage of study, a simple AF determined every 250 ms within windows of 1.5 s was used. Because at this stage of analysis we just want to compare two methods for envelope determining, so such approach, which is a compromise between ensuring the continuity of measurements and their accuracy, is the most right. Additionally, the AF threshold relating to low level of signal quality and thus defining the FHR signal loss was empirically determined. For this purpose, the most disturbed signal fragments were analyzed, in which in the absence of fetal heart episodes, many local AF peaks were observed, whose amplitude usually did not exceed 0.1. Therefore, if the amplitude of AF peak did not exceed P TH =0.1, a loss marker was set in the resulting FHR signal. Such AF threshold is more rigid than 0.15 proposed in the study [36]. 2.2.2. Periodicity Determination Two different methods for determining the autocorrelation function were used. The first is the detection of the AF peak calculated arbitrarily every 250 ms. In the second case, AF is calculated more frequently with every 25 ms, and the outlier values most often caused by temporary interferences, are rejected from a group of measurements. An important factor affecting both the measurement accuracy of T i interval value and the level of FHR signal loss is the AF window width. It is required Sensors 2020,20, 4079 8 of 25 that the window should cover at least two of the same cardiac phases belonging to two consecutive cardiac cycles. The multiphase motion process observed in US signal means that determination of the T i interval does not require a window including two full cardiac cycles. Even if certain cardiac phase will occur only once in a narrower window, then most likely the other phase will be represented twice, allowing correct measurement. A narrow AF window means minimization of the averaging effect over several cycles, and thus a higher measurement accuracy. However, too narrow of a window may prevent correct measurement, which will result in FHR signal loss. Jezewski [ 36 ] set the minimum window width at 1.5 of the last measured interval duration. The advantage of this adaptive window width selection is reduction of necessary calculations, especially in case of short T i intervals, which may be important in practical implementation in battery-powered mobile fetal monitors. On the other hand, so strong dependence of the calculation process on previously measured values can cause a permanent switch to an incorrect range of Tivalues. The algorithm for the periodicity determination assumes that when the signal quality is satisfactorily high, the maximum of AF found is likely to correspond to the T i interval value in US signal. To prevent erroneous measurements for low-quality signal fragments (AF peak amplitude <0.5), the prediction function (PF) was proposed, where recently determined interval value is taken as an expected value for the next one. Such AF peak amplitude used in prediction function was determined empirically based on analysis of problematic cases observed in available signals. In the absence of interferences, the most frequently observed AF peak amplitudes were usually above 0.6. In the presence of additional disturbing episodes observed due to fetal movement activity, the AF peak amplitude was just below 0.5. The value of the AF is corrected by using a trapezoidal PF window, with a selected width of the lower base, a base ratio of 1:4 and the center located at the point τ =T i−1 (Figure 5). The proposed solution is a modification of the approach from [ 36 ], where a triangular window was used, which however exposed too much the last measured value. In cases of increased heart rate changes, e.g., on acceleration and deceleration slopes, where the differences of successive intervals are significant, the use of a triangular window caused erroneous measurements. Sensors 2020, 20, x FOR PEER REVIEW 9 of 25 Figure 5. Idea of a trapezoidal window applied for the prediction function prediction function (PF). The upper graph (A) shows the autocorrelation function with PF, whose center is determined by the last correctly determined value of the Ti−1 interval. The lower graph (B) shows that correction of autocorrelation function (AF) by trapezoidal PF which enables determination of the correct periodicity value (τ2) and thus the instantaneous measurement Fj, despite the distortions in the US signal, which gives false AF maximum (τ1) before PF application. 2.2.3. Extraction of Time Series of Events The first method (D1) for determining the fetal heart rate signal from the Doppler US signal provides the FHR value every 250 ms (like in classical monitors), after which a correction of duplicate measurements proposed in [65] is applied (Figure 6). The AF was calculated every 250 ms, but sometimes also more often (every 25 ms) and then the representative FHR value was the median of 10 consecutive instantaneous measurements Fj. According to [36], the minimum window width was set at 800 ms, which is 1.5 of maximum interval value (533 ms) in the reference signal, and other tested widths were: 1, 1.5, 2 and 3 s. The impact of different PF window widths (S) on the final measurements was also tested with 125, 250 and 500 ms, which were adopted basing on visual analysis of the AF. Value S = 0 means that PF was not used, and then the position of maximum peak of AF (in the physiologically acceptable range from 30 to 2000 ms) defined periodicity. Figure 6. General block diagram of the D1 algorithm for determining the instantaneous value of fetal heart rate (Ti interval) based on a raw Doppler US signal. The FHR values determined every 250 ms, underwent correction of duplicate measurements to obtain a signal in the form of series of events according to the authors’ algorithm [65]. The aim of the algorithm for identifying and correcting duplicate measurements is determination of the number of true heart cycles represented by a sequence of instantaneous measurements. For each analyzed Figure 5. Idea of a trapezoidal window applied for the prediction function prediction function (PF). The upper graph ( A ) shows the autocorrelation function with PF, whose center is determined by the last correctly determined value of the T i−1 interval. The lower graph ( B ) shows that correction of autocorrelation function (AF) by trapezoidal PF which enables determination of the correct periodicity value ( τ 2) and thus the instantaneous measurement F j , despite the distortions in the US signal, which gives false AF maximum (τ1) before PF application. Sensors 2020,20, 4079 9 of 25 The new trapezoidal window causes that the AF amplitude for shifts τ close to the values of recently determined T i remains unchanged. While for less likely intervals ( τ values significantly different from T i−1 ) the amplitude is reduced according to the applied PF. In addition, it was assumed that if the current interval was determined using a prediction algorithm, then it does not modify the PF position for the next iteration. Finally, the T i takes the value of τ for which the AF( τ ) reaches maximum value. In the next part of investigation, we developed two independent methods based on the autocorrelation function, which allow the resulting FHR signal to get the required form of time series of events. For this purpose, a method of correcting duplicate measurements [ 65 ] was used, as well as a segmentation of instantaneous periodicity measurements. 2.2.3. Extraction of Time Series of Events The first method (D1) for determining the fetal heart rate signal from the Doppler US signal provides the FHR value every 250 ms (like in classical monitors), after which a correction of duplicate measurements proposed in [ 65 ] is applied (Figure 6). The AF was calculated every 250 ms, but sometimes also more often (every 25 ms) and then the representative FHR value was the median of 10 consecutive instantaneous measurements F j . According to [ 36 ], the minimum window width was set at 800 ms, which is 1.5 of maximum interval value (533 ms) in the reference signal, and other tested widths were: 1, 1.5, 2 and 3 s. The impact of different PF window widths (S) on the final measurements was also tested with 125, 250 and 500 ms, which were adopted basing on visual analysis of the AF. Value S=0 means that PF was not used, and then the position of maximum peak of AF (in the physiologically acceptable range from 30 to 2000 ms) defined periodicity. Sensors 2020, 20, x FOR PEER REVIEW 9 of 25 Figure 5. Idea of a trapezoidal window applied for the prediction function prediction function (PF). The upper graph (A) shows the autocorrelation function with PF, whose center is determined by the last correctly determined value of the Ti−1 interval. The lower graph (B) shows that correction of autocorrelation function (AF) by trapezoidal PF which enables determination of the correct periodicity value (τ2) and thus the instantaneous measurement Fj, despite the distortions in the US signal, which gives false AF maximum (τ1) before PF application. 2.2.3. Extraction of Time Series of Events The first method (D1) for determining the fetal heart rate signal from the Doppler US signal provides the FHR value every 250 ms (like in classical monitors), after which a correction of duplicate measurements proposed in [65] is applied (Figure 6). The AF was calculated every 250 ms, but sometimes also more often (every 25 ms) and then the representative FHR value was the median of 10 consecutive instantaneous measurements Fj. According to [36], the minimum window width was set at 800 ms, which is 1.5 of maximum interval value (533 ms) in the reference signal, and other tested widths were: 1, 1.5, 2 and 3 s. The impact of different PF window widths (S) on the final measurements was also tested with 125, 250 and 500 ms, which were adopted basing on visual analysis of the AF. Value S = 0 means that PF was not used, and then the position of maximum peak of AF (in the physiologically acceptable range from 30 to 2000 ms) defined periodicity. Figure 6. General block diagram of the D1 algorithm for determining the instantaneous value of fetal heart rate (Ti interval) based on a raw Doppler US signal. The FHR values determined every 250 ms, underwent correction of duplicate measurements to obtain a signal in the form of series of events according to the authors’ algorithm [65]. The aim of the algorithm for identifying and correcting duplicate measurements is determination of the number of true heart cycles represented by a sequence of instantaneous measurements. For each analyzed Figure 6. General block diagram of the D1 algorithm for determining the instantaneous value of fetal heart rate (Tiinterval) based on a raw Doppler US signal. The FHR values determined every 250 ms, underwent correction of duplicate measurements to obtain a signal in the form of series of events according to the authors’ algorithm [ 65 ]. The aim of the algorithm for identifying and correcting duplicate measurements is determination of the number of true heart cycles represented by a sequence of instantaneous measurements. For each analyzed sequence, the Min and Max values are determined, defining two different numbers of such intervals. If the rounded Min and Max values are equal, then this value represents the only correct solution. In the case of different values, determination of the number of real heart intervals is based on analysis of time dependencies as well as assessment of probability of a given solution existence. In the D2 method, the autocorrelation function is determined in a window shifted with 25 ms (Figure 7). Subsequent values of instantaneous periodicity F j are determined using the prediction function (as in the D1 method). The next stage is segmentation of measurements F j , which identifies a group of subsequent measurements, which will then be used to calculate the resulting interval T i . Segmentation relies on continuous computation of a median value until a time being the product of number of measurements and repetition period (25 ms) exceeds the current median value. Then the median becomes the searched T i value and the procedure is repeated. Unfortunately, the obtained results depend on a starting point determining the phase shift between segments and series of F j . Two approaches are proposed to reduce this unwanted effect. Sensors 2020,20, 4079 16 of 25 width (from 1.2 to 1.26 ms). For a narrower window, some of the erroneous measurements, but fulfilling van Geijn criteria, are considered correct, causing a slight increase in the error value. In turn, a wider window means that the resulting measurement represents a greater number of consecutive cardiac cycles, which leads to averaging of T i values and to error increase. Ultimately, the best results both for |∆Ti| and signal loss, were obtained for AF window of 1 s with repetition of measurements every 25 ms and using an additional prediction function. Finally, the T i interval determination using two different signal envelopes (low-pass filtering with a 50 Hz cutoffand the Hilbert transform), was used in the algorithms to extract the FHR signal in the form of time series of events. After correction of duplicate measurements [ 65 ], for envelope with low-pass filtering, the error |∆Ti| increased from 1.24 ms for K =25 ms (Table 3), to 1.98 ms for the representation as time series of events. For the Hilbert transform, the error values were 1.25 ms and 1.91 ms after the correction algorithm (Table 3). As an input signal for D2, we used two optimal (selected within D1) signal envelopes determined by 50 Hz low-pass filtering and the Hilbert transform. The AF was calculated every 25 ms in one-second window, whereas PF with 500 ms window was used to determine the Tivalue. During the studies it was tested how the obtained results depend on the starting point selection, i.e., the location of the first AF window in the envelope signal. Table 4presents the |∆Ti| and FHR_U signal loss for different starting points: arbitrarily selected at the signal beginning, as well as five other points determined using the RMS method calculated in windows: 100, 200, 250, 500 and 1000 ms. In addition, we estimated usability of the developed matching algorithm (MAA) with continuous correction of the determined sequence of events in relation to the series of instantaneous measurements Fj. Table 4. Interval measurements errors and the FHR_U signal loss for different starting points, with and without additional matching algorithm MAA. RMS Window Width (ms) |∆Ti|(ms) FHR_U Signal Loss (% ) Without MAA MAA Without MAA MAA Filter 50 Hilbert Filter 50 Hilbert Filter 50 Hilbert Filter 50 Hilbert –1.62 #1.63 1.44 1.45 0.43 0.43 0.43 0.44 100 1.64 1.67 1.43 1.45 0.38 0.36 0.40 0.44 200 1.67 1.70 1.48 1.43 0.38 0.42 0.40 0.43 250 1.66 1.69 1.53 1.46 0.41 0.41 0.39 0.43 500 1.68 1.61 1.45 1.35 x0.47 0.43 0.44 0.41 x 1000 1.69 1.72 1.48 1.44 0.41 0.43 0.41 0.42 #mean value xminimal Tierror value determining the optimal values of control parameter. The starting point selection practically did not affect the accuracy of the T i measurement. Without additional matching algorithm, regardless of the envelope type, the mean error |∆Ti| ranged from 1.62 to 1.72 ms. The reason for such small errors for different starting points is that very often all RMS functions provided the same starting point for a given signal. In addition, there were cases where different starting points after some time of algorithm operation converged at one point, and the analysis comprised identical envelope fragments. Regardless of the envelope type, much better results were noted with additional MAA—the |∆Ti| ranged from 1.35 to 1.53 ms. The algorithms used for selection of the starting point and further correction MAA did not affect a level of the signal loss, whose value was around 0.4%. 3.2. Final Assessment Section 3.1 presents the results regarding determination of control parameters for the algorithms D1 and D2. The final results of the algorithm performance were showed in relation to the original data from the fetal monitor. Table 5presents final results obtained for FHR determination methods based Sensors 2020,20, 4079 17 of 25 on the Doppler ultrasound signal. For comparative purposes, the column FHR_CTG contains results obtained for a signal provided by a fetal monitor used in this study and thus calculated by its built-in US signal processing algorithms. Table 5. Final results obtained for various methods for FHR_U signal determination based on a raw US signal in relation to the reference FHR_E from direct FECG. Parameter FHR Determination Method FHR_CTG FHR_U D1 FHR_U D2 ∆Ti(ms) 0.13 0.11 0.04 SD_∆Ti(ms) 4.01 3.53 2.14 |∆Ti|(ms) 2.78 1.91 1.35 FHR signal loss (%) 0.0 0.47 0.41 The mean value of differences ∆Ti for both developed methods did not exceed 0.11 ms. In turn, standard deviation was 3.53 ms for D1 method and 2.14 ms for D2. The D2 method involves multiple AF calculation and segmentation of instantaneous measurements based on measured periodicity value. The method ensures that 95% of measurement errors do not exceed 4.28 ms. The mean |∆Ti| is 1.35 ms and is more than twice lower than for the original FHR_CTG signals determined inside a fetal monitor. The developed methods are characterized by a very low level of the signal loss, which for D2 is only 0.41%. Lack of signal loss for FHR_CTG results from wide AF window used, which is responsible for the measurement averaging, thus increasing the error value (Table 5). A very important criterion for assessing the investigated methods are the errors of instantaneous FHR variability determination. Table 6lists relative error values for the original FHR signal determined by US signal processing algorithms built-in a fetal monitor, the D1 method based on determination of subsequent FHR values every 250 ms and with duplicate measurements correction algorithm [ 65 ], as well as the D2 method in which time series of events is determined already at the stage of US signal analysis. Table 6. Relative errors δ for STV and LTV variability indices, determined for original signals FHR_CTG from a fetal monitor output, and for FHR_U from the new Doppler ultrasound signal analysis methods, in relation to the FHR_E reference signal from FECG. δ(%) FHR Determination Method FHR_CTG FHR_U D1 FHR_U D2 STV δS_DAW −66.2 ±7.99 #−6.4 ±6.25 −5.0 ±4.75 δS_YEH −50.8 ±11.07 −3.0 ±7.69 −2.4 ±6.78 δS_HAA −80.0 ±7.27 −13.3 ±12.32 −8.9 ±12.02 δS_ZUG −51.5 ±9.28 −2.1 ±9.09 −1.9 ±8.15 δS_HUE −91.3 ±8.73 −11.9 ±16.07 −10.1 ±12.00 δS_DAL −65.5 ±7.59 −6.9 ±6.25 −5.4 ±4.72 LTV δL_DAW −10.3 ±9.00 −0.7 ±6.47 −0.1 ±5.76 δL_YEH −6.3 ±6.78 −0.3 ±2.79 −0.1 ±2.36 δL_HAA −2.3 ±11.16 −1.4 ±6.89 −0.9 ±5.09 δL_ZUG −6.7 ±7.26 −0.2 ±2.50 −0.1 ±2.19 δL_HUE 9.80 ±16.95 −1.3 ±10.22 −1.3 ±8.04 δL_DAL −5.9 ±6.93 −0.1 ±2.87 −0.1 ±2.14 #mean value ±SD. Sensors 2020,20, 4079 18 of 25 The values of all variability indices determined for the FHR_CTG signal are significantly lower (as negative values) than those obtained for the FHR_E reference signal. Errors of STV indices are in the range from − 50.8% for the S_YEH index to as much as − 91.3% for the S_HUE index. Such large errors are caused by duplicate measurements and the intervals averaging effect decreasing beat-to-beat variability. This is confirmed by the results of D1 method, in which based on the developed algorithm of US signal analysis, instantaneous FHR values are calculated similarly to the fetal monitor (every 250 ms), and then corrected by detecting and removing duplicate measurements. This leads to significant error reduction, also for LTV indices. The D2 method provides the best representation of the real FHR signal variability, whose reference pattern comes from the signal derived from a direct FECG. For the most clinically relevant STV indices describing beat-to-beat variability, the error value ranged from − 1.9% for S_ZUG to −10.1% for the S_HUE index. 4. Discussion The paper proposed new methods for determining a fetal heart rate signal based on a detailed analysis of the raw Doppler ultrasound signal. The primary requirement for these methods was that the FHR signal should be in the form of time series of events reflecting subsequent cardiac cycles—Tiintervals. Experimental studies were carried out based on simultaneously recorded Doppler ultrasound signal, original FHR_CTG signal calculated on its basis inside a fetal monitor, as well as the FHR_E signal determined from a direct fetal electrocardiogram (being the source of reference information). The effectiveness of the developed methods was evaluated on the basis of measurement accuracy of T i intervals between successive heartbeats, combined with estimation of the resulting FHR signal loss. In turn, the second stage of assessment concerned the accuracy evaluation of the instantaneous FHR variability, based on dedicated indices describing short-term and long-term variability. There are various ways to assess the quality of FHR signals from the outputs of popular ultrasound-based fetal monitors relative to the reference FECG. The comparison of signals in metrological aspect concerns both visual assessment of the FHR traces and estimation of statistical parameters describing the signals, as well as comparison of the synchronized measurements—T i intervals. On the other hand, the verification of the ultrasound method in a clinical aspect consists in comparing diagnostically significant patterns determined on the FHR basis, such as: baseline, acceleration and deceleration patterns or indices describing instantaneous variability. The first reports regarding the FHR signal recorded by an US method were published when fetal monitors became popular-in the seventies. Already in 1976, Lauersen [ 72 ] studied the differences between simultaneously recorded FHR signals via ultrasound and electrocardiography. The standard deviation for differences in corresponding measurements was as much as 9.8 bpm, mainly because of a lack of precise signal synchronization. Dawes [ 73 ] also omitted the synchronization problem—and when compared the FHR mean values representing one-hour signals—he obtained a correlation of 0.9. In addition, he assessed the FHR short-term variability using the S_DAW index, for which a correlation was 0.63. The difference between mean values of the S_DAW index: 2.3 ms for FECG and 4.6 ms for US, was statistically significant. Lawson [ 74 ] used the same method, but with a different fetal monitor. For FHR signals from the FECG channel, the mean S_DAW value was 2.1 ms, while for the US it was three times higher − 6.2 ms. The obtained errors [ 72 – 74 ] resulted from using monitors with US channel of the so-called first generation, where T i intervals were determined solely on the peak detection. In another work [ 75 ], Lawson has already used a monitor with a US channel of second generation, in which AF has already been applied to measure a periodicity. Furthermore, in this case, the difference was statistically significant, although much smaller and with the opposite trend—the S_DAW index was 1.29 ms for the US, which was 35% lower than its value of 1.98 ms for the FECG channel. The reason for underestimation was the autocorrelation function leading to averaging of Tiintervals. Sensors 2020,20, 4079 19 of 25 Jezewski made a comprehensive assessment of the accuracy of FHR measurements [ 62 ], comparing 21,941 pairs of T i intervals recorded simultaneously with US channel of fetal monitor MT-430 (Toitu) and direct fetal electrocardiography. Reference FECG was used to eliminate duplicate measurements, which after synchronization served as a clock signal for the FHR signal from a fetal monitor. This correction method has enabled a reliable assessment, but it was impossible to apply in clinical practice. For higher reliability of results, the FHR reference signal was determined by two independent algorithms from the FECG, and during final analysis the measurements for which the difference between them exceeded one millisecond were rejected. The mean error |∆Ti| between US and FECG was 2.98 ms. In turn, the errors of variability indices determination ranged from − 2.0 to − 9.3% for LTV and from − 5.9 to as much as − 39.5% for STV. No results were provided for uncorrected data from a fetal monitor, which would allow estimating the impact of duplicate measurements on errors of the instantaneous variability assessment. Other studies concerned methods for determining FHR signals in the form of time series of events from measurements occurring every 250 ms and recorded with a classical fetal monitor [ 63 – 65 ]. Research conducted in [ 65 ] clearly showed that, in relation to the FECG the fetal heart rate signal containing duplicate measurements introduces errors, especially when determining the STV indices (from − 51.5% to as much as − 91.3%). Effective removal of duplicate measurements—conversion of the FHR signal to the form of time series of events—allowed to reduce these values (from − 28.0% to − 74.0%), which unfortunately still does not allow for reliable analysis of the FHR variability. This was confirmed by the results in [ 62 ], but although errors of the algorithm for duplicate measurements correction were eliminated by FECG signal timing, the errors of STV indices still had unacceptable values. The main reason was a use of too wide AF windows being responsible for measurements averaging. Therefore, only the development of more effective algorithms to analyze the raw US signal may improve the accuracy of FHR determination and thus the quantitative assessment of instantaneous variability. Therefore, the authors developed two new methods for extraction of time series of events—the occurrences of fetal heartbeats. The first is D1 which is based on the duplicate measurements correction previously developed by the authors [ 65 ]. The obtained FHR signal and STV indices are much accurate than results obtained for the FHR_CTG signal from a fetal monitor. The mean error of interval measurements |∆Ti| was 1.91 ms, against 2.78 ms for the signal from a monitor output. For LTV indices, errors ranged from − 0.1% to − 1.4%. In contrast, for STV indices errors ranged from − 2.1% to − 13.3%. It is worth to note that, errors obtained for the original FHR_CTG signals were from − 50.8% to − 91.3%. In the second D2 method, required events are results of segmentation of instantaneous FHR measurements, based on the analysis of measured values. This method provides a significantly lower interval measurement error of 1.35 ms. The increase in measurement accuracy does not occur at the expense of signal loss level equal to 0.41%. Errors of determining LTV indices using D2 were very similar to the D1 method and ranged from − 0.1% to − 1.3%. In the case of STV indices, errors were always lower than for D1 and ranged from − 1.9% to − 10.1%. The obtained results clearly confirm that the D2 method is the best in terms of both the accuracy of measurement and reliability of estimation of FHR instantaneous variability. In general, the research carried out in the work covered a very wide range of issues, from an initial analysis of raw Doppler ultrasound signal after demodulation to the final measurement of FHR values. Therefore, the obtained results were influenced by many different control parameters. Some of them as: threshold of the AF peak amplitude defining the FHR signal loss and threshold deciding about the measurements support by an additional prediction function, were determined empirically basing on our previous experience [ 36 , 70 ]. Significantly larger group are parameters whose optimal values have been determined on the basis of results of research carried out as part of this study. These are: Doppler signal bandwidths, control parameters in determining the envelope signal, window parameters of the autocorrelation and prediction functions, as well as the width of the RMS window used in the final Sensors 2020,20, 4079 20 of 25 segmentation process. The process of determining all these parameters and evaluation of their impact on final results is documented in Section 3.1 Study of sensitivity. The analyzed signal bandwidth has a significant impact on the results obtained. It should cover the range from 300 to 600 Hz corresponding to fetal heart valves. The use of short AF window (1 s) and measurement support with a trapezoidal prediction function turned out to be very important, which ensured high accuracy of T i measurement and low signal loss in the resulting FHR. Nevertheless, for reliable assessment of indices describing beat-to-beat variability, a signal in the form of time series of events is required. In this case, better results were obtained using segmentation of instantaneous measurements Fj. As the final results showed, it was possible to obtain a high accuracy of T i measurement and what is very important at a low level of FHR signal loss. However, we are far from claiming that the proposed algorithms are ready-to-use as completely reliable solutions. The work is primarily of a research nature, thanks to which it was shown that it is not the Doppler ultrasound signal itself (as representing less accurate mechanical activity of the fetal heart), that is the source of errors, but the methods of its analysis used. The problem of determining the fetal heart rate in the form of time series of events for the US signal, which source was the second-generation HP 8040A monitor, was raised by Peters [ 66 ]. He applied a low-pass filter with a cutofffrequency of two hertz, and then the obtained markers representing the subsequent fetal heartbeats were used to separate fragments of an envelope signal, for which the Ti intervals were determined using AF. The 100-s signal was processed, and the results were referred to a direct FECG. The correlation coefficient between FHR signals was high and reached 0.977. Unfortunately, there was no assessment of the interval measurement accuracy, only a graph illustrating differences between the measurements for both methods against their means was provided. A detailed assessment of an interval measurement error and STV indices determination based on an US signal analysis is provided in [ 67 ]. Referring to a simultaneously recorded FECG signal, in the US signal an AF window with the width of two full cardiac cycles was positioned. The tests concerned different values of window shift relative to an initial position, determined by the occurrence of the R wave. The minimum error |∆Ti| was 1.59 ms, while the minimum error of the S_DAW index was 7.12%. Unfortunately, the work has not any practical value because the US signal analysis requires to be supported by the FECG signal. Jezewski also dealt with analysis of raw US signal to determine FHR in the form of time series of events [ 36 ]. The signals were recorded during labor using a prototype instrumentation, while the reference was a direct FECG signal. An AF window was used, whose width and step were changed adaptively depending on the last determined interval. The best results were obtained for a window width equal to twice the width of a previous interval, and for a step ensuring five measurements of each interval. The mean error |∆Ti| was 1.91 ms, with FHR signal loss of 1.6%. In addition, an error of −6.9% was determined for the S_HAA variability index. Although investigations [ 36 ] were conducted using other database than in this study, in both cases a direct FECG was the reference method, and the number of analyzed intervals was similar—8945 against 7498 intervals examined by the authors. Therefore, the statistical significance of both studies represents similar level. This, in turn, allows to state that the results obtained in this work show higher accuracy of periodicity measurement and lower level of FHR signal loss. For the D2 method, the error in interval measuring was only 1.35 ms, with a signal loss of 0.41%, which confirms a high efficiency of the developed solution. However, lower interval error did not lead to decrease of the S_HAA variability index error, as it was − 8.9%, while in [ 36 ] only − 6.9%. Unfortunately, in [ 36 ] only S_HAA index was analyzed, whereas it is obvious that various indices react in different ways to the error in measuring the interval. Sensors 2020,20, 4079 21 of 25 5. Conclusions Using a popular fetal monitor based on the Doppler ultrasound technique, it is possible to analyze the fetal heart rate signal only in a classical approach, which consists in determining the FHR baseline for recognizing the patterns of tachycardia/bradycardia or acceleration/deceleration. However, it is almost impossible to reliably assess the instantaneous FHR variability that is a part of automated analysis carried out in computer-aided fetal monitoring system. The presented improvement was ensured by FHR signal correction consisting in removing of duplicate measurements. Thanks to it, the signal took the form of time series of events, as in electrocardiography. Unfortunately, from the point of view of fetal state assessment, errors of determining indices for instantaneous FHR variability remained at an unacceptable level. This means that another factor affects fetal monitors—the built-in algorithm for determining subsequent FHR values, focusing too much on being resistant to disturbances from fetal movements. Our research revealed that the Doppler ultrasound signal itself is not a source of errors. It is too long window in correlation methods that causes an effect of “averaging” several consecutive fetal heart cycles. A comprehensive study of the Doppler US signal describing the mechanical activity of a fetal heart was carried out. The best results in terms of interval measurement accuracy and evaluation of instantaneous variability were obtained using the authors’ methods analyzing the raw Doppler ultrasound signal. The D2 method based on segmentation of instantaneous measurements particularly allows obtaining accuracy close to the electrocardiographic method. The results unambiguously confirm that it is possible to analyze the Doppler ultrasound signal in a way that allows effective determination of fetal heart rate in the form of time series of events—significantly improving the reliability of determining the indices for FHR instantaneous variability description. Author Contributions: Conceptualization, T.K., A.M. and M.J.; data curation, J.J., K.H. and R.M.; methodology, T.K., J.J. and R.C.; software, A.M., M.J. and J.W.; supervision, J.J. and R.M.; validation, A.M., J.W. and R.C.; writing—original draft, T.K., A.M. and K.H.; writing—review & editing, T.K., M.J., J.J. and K.H. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Acknowledgments: This scientific research work was supported by the National Science Center in Poland under the grant 2017/27/B/ST6/01989, as well as from the Ministry of Science and Higher Education funding for statutory activities BK-2020. Conflicts of Interest: The authors declare no conflicts of interest. References 1. Horoba, K.; Wrobel, J.; Jezewski, J.; Kupka, T.; Roj, D.; Jezewski, M. Automated detection of uterine contractions in tocography signals–comparison of algorithms. Biocybern. Biomed. Eng. 2016 ,36, 610–618. [CrossRef] 2. Jezewski, J.; Matonia, A.; Czabanski, R.; Horoba, K.; Kupka, T. Classification of Uterine Electrical Activity Patterns for Early Detection of Preterm Birth. In Computer Recognition Systems 8–CORES 2013; Advances in Intelligent Systems and Computing; Burduk, R., Jackowski, K., Kurzynski, M., Wozniak, M., Zolnierek, A., Eds.; Springer: Heidelberg, Germany, 2013; Volume 226, pp. 559–568. 3. Jezewski, J.; Wrobel, J.; Horoba, K.; Gacek, A.; Sikora, J. Fetal heart rate variability: Clinical experts versus computerized system interpretation. In Proceedings of the 24th International Conference of IEEE Engineering in Medicine and Biology Society, Houston, TX, USA, 23–26 October 2002; pp. 1617–1618. 4. Horoba, K.; Jezewski, J.; Matonia, A.; Wrobel, J.; Czabanski, R.; Jezewski, M. Early predicting a risk of preterm labour by analysis of antepartum electrohysterograhic signals. Biocybern. Biomed. Eng. 2016 ,36, 574–583. [CrossRef] 5. Sikora, J.; Matonia, A.; Czabanski, R.; Horoba, K.; Jezewski, J.; Kupka, T. Recognition of Premature Threatening Labour Symptoms from Bioelectrical Uterine Activity Signals. Arch. Perinat Med. 2011,17, 97–103. Sensors 2020,20, 4079 22 of 25 6. Kording, F.; Schoennagel, B.; Lund, G.; Ueberle, F.; Jung, C.; Adam, G.; Yamamura, J. Doppler ultrasound compared with electrocardiogram and pulse oximetry cardiac triggering: A pilot study. Magn. Reson. Med. 2014,74, 1257–1265. [CrossRef] 7. Jaros, R.; Martinek, R.; Kahankova, R. Non-Adaptive Methods for Fetal ECG Signal Processing: A Review and Appraisal. Sensors 2018,18, 3648. [CrossRef] 8. Van Scheepen, J.A.M.; Koster, M.P.H.; Vasak, B.; Redman, C.; Franx, A.; Georgieva, A. Effect of signal acquisition method on the fetal heart rate analysis with phase rectified signal averaging. Physiol. Meas. 2016 , 37, 2245–2259. [CrossRef] 9. Jezewski, J.; Horoba, K.; Roj, D.; Wrobel, J.; Kupka, T.; Matonia, A. Evaluating the fetal heart rate baseline estimation algorithms by their influence on detection of clinically important patterns. Biocybern. Biomed. Eng. 2016,36, 562–573. [CrossRef] 10. Jezewski, J.; Horoba, K.; Matonia, A.; Gacek, A.; Berny´s, M. A new approach to cardiotocographic fetal monitoring based on analysis of bioelectrical signals. In Proceedings of the 25th International Conference of IEEE Engineering in Medicine and Biology Society, Cancum, Mexico, 1–4 September 2003; pp. 3145–3149. 11. Jezewski, J.; Wrobel, J.; Matonia, A.; Horoba, K.; Martinek, R.; Kupka, T.; Jezewski, M. Is abdominal fetal electrocardiography an alternative to Doppler ultrasound for FHR variability evaluation? Front. Physiol. 2017,8, 305. [CrossRef] 12. Kahankova, R.; Jezewski, J.; Nedoma, J.; Jezewski, M.; Fajkus, M.; Kawala-Janik, A.; Wen, H.; Martinek, R. Influence of gestation age on the performance of adaptive systems for fetal ECG extraction. Adv. Electr. Electron. Eng. 2017,15, 491–501. [CrossRef] 13. Jezewski, J.; Matonia, A.; Kupka, T.; Roj, D.; Czaba´nski, R. Determination of the fetal heart rate from abdominal signals: Evaluation of beat-to-beat accuracy in relation to the direct fetal electrocardiogram. Biomed. Eng. 2012,57, 383–394. [CrossRef] 14. Kahankova, R.; Jaros, R.; Martinek, R.; Jezewski, J.; Wen, H.; Jezewski, M.; Kawala-Janik, A. Non-Adaptive Methods of Fetal ECG Signal Processing. Adv. Electr. Electron. Eng. 2017,15, 476–490. [CrossRef] 15. Martinek, R.; Kahankova, R.; Jezewski, J.; Jaros, R.; Mohylova, J.; Fajkus, M.; Nedoma, J.; Janku, P.; Nazeran, H. Comparative Effectiveness of ICA and PCA in Extraction of fECG from aECG Signals: Towards Multichannel Non-Invasive fHR Monitoring. Front. Physiol. 2018,9, 648. [CrossRef] [PubMed] 16. Docker, M.F. Doppler ultrasound monitoring technology. Br. J. Obstet. Gynaecol. 1993 ,100, 18–20. [CrossRef] [PubMed] 17. Khandoker, A.H.; Kimura, Y.; Ito, T.; Palaniswami, M. Non-invasive determination of electromechanical time intervals of cardiac cycle using abdominal ECG and Doppler ultrasound signals from fetal hearts. Comput. Cardiol. 2007,34, 657–660. 18. Abdulhay, E.W.; Oweis, R.J.; Alhaddad, A.M.; Sublaban, F.N.; Radwan, M.A.; Almasaeed, H.M. Non-Invasive Fetal Heart Rate Monitoring Techniques: Review Article. Biomed. Sci. Eng. 2014,2, 53–67. 19. Alnuaimi, S.; Jimaa, S.; Khandoker, A.H. A Review of Fetal cardiac Doppler Signal Processing for Screening Foetal Well Being. Front. Bioeng. Biotechnol. 2017,5, 82. [CrossRef] 20. Kording, F.; Schoennagel, B.; Tavares de Sousa, M.; Fehrs, K.; Adam, G.; Yamamura, J.; Ruprecht, C. Evaluation of a Portable Doppler Ultrasound Gating Device for Fetal Cardiac MR Imaging: Initial Results at 1.5T and 3T. Magn. Reson. Med. Sci. 2018,17, 308. [CrossRef] 21. Mert, A.; Sezdi, M.; Akan, A. A test and simulation device for Doppler-based fetal heart rate monitoring. Turk. J. Electr. Eng. Comput. Sci. 2015,23, 1187–1194. [CrossRef] 22. Karjadi, M.; Salahuddin, N.S.; Wibowo, E.P.; Afandi, H. Digital Filter Design of Infinite Impulse Response (IIR) Infrasound to Detect Fetal Heart Rate. Int. J. Eng. Res. Sci. 2016,2, 25–38. 23. Kupka, T.; Je˙zewski, J.; Matonia, A.; Horoba, K.; Wr ó bel, J. Timing events in Doppler ultrasound signal of fetal heart activity. In Proceedings of the 26th International Conference of IEEE Engineering in Medicine and Biology Society, San Francisco, CA, USA, 1–4 September 2004; pp. 337–340. 24. Voicu, I.; Girault, J.M.; Roussel, C.; Decock, A.; Kouame, D. Robust estimation of fetal heart rate from US Doppler signals. Phys. Procedia 2010,3, 691–699. [CrossRef] 25. Marzbanrad, F.; Kimura, Y.; Funamoto, K.; Sugibayashi, R.; Endo, M.; Ito, T.; Palaniswami, M.; Khandoker, A.H. Automated Estimation of Fetal Cardiac Timing Events From Doppler Ultrasound Signal Using Hybrid Models. IEEE J. Biomed. Health Inform. 2014,18, 1169–1177. [CrossRef] [PubMed] Sensors 2020,20, 4079 23 of 25 26. Marzbanrad, F.; Kimura, Y.; Funamoto, K.; Oshio, S.; Endo, M.; Sato, N.; Palaniswami, M.; Khandoker, A.H. Model-Based Estimation of Aortic and Mitral Valves Opening and Closing Timings in Developing Human Fetuses. IEEE J. Biomed. Health Inform. 2016,20, 240–248. [CrossRef] [PubMed] 27. Marzbanrad, F.; Stroux, L.; Clifford, G.D. Cardiotocography and beyond: A review of one-dimensional Doppler ultrasound application in fetal monitoring. Physiol. Meas. 2018,39, 08TR01. [CrossRef] [PubMed] 28. Khandoker, A.; Kimura, Y.; Ito, T.; Sato, N.; Okamura, K.; Palaniswami, M. Antepartum non-invasive evaluation of opening and closing timings of the cardiac valves in fetal cardiac cycle. Med. Biol. Eng. Comput. 2009,47, 1075–1082. [CrossRef] 29. Lee, C.S.; Masek, M.; Lam, C.P.; Tan, K.T. Towards Higher Accuracy and Better Noise-Tolerance for Fetal Heart Rate Monitoring using Doppler Ultrasound. In Proceedings of the IEEE Region 10 Conference (TENCON), Singapore, 23–26 November 2009; pp. 1–6. 30. Voicu, I.; Menigot, S.; Kouame, D.; Girault, J.M. New estimators and guidelines for better use of fetal heart rate estimators with doppler ultrasound devices. Comput. Math. Method Med. 2014 ,2014. [CrossRef] [PubMed] 31. Shakespeare, S.A.; Crowe, J.A.; Hayes-Gill, B.R.; Bhogal, K.; James, D.K. The information content of Doppler ultrasound signals from the fetal heart. Med. Biol. Eng. Comput. 2001,39, 619–626. [CrossRef] 32. Wrobel, J.; Roj, D.; Jezewski, J.; Horoba, K.; Kupka, T.; Jezewski, M. Evaluation of the robustness of fetal heart rate variability measures to low signal quality. J. Med. Imaging Health Inform. 2015 ,5, 1311–1318. [CrossRef] 33. Hamelmann, P.; Mischi, M.; Kolen, A.F.; van Laar, J.O.E.H.; Vulings, R.; Bergmans, J.W.M. Fetal Heart Rate Monitoring Implemented by Dynamic Adaptation of Transmission Power of a Flexible Ultrasound Transducer Array. Sensors 2019,19, 1195. [CrossRef] 34. Vlachos, M.; Yu, P.; Castelli, V. On Periodicity Detection and Structural Periodic Similarity. In Proceedings of the 5th SIAM International Conference on Data Mining, Newport Beach, CA, USA, 21–23 April 2005; pp. 449–460. 35. Hamelmann, P.; Vulings, R.; Kolen, A.F.; Bergmans, J.W.M.; van Laar, J.O.E.H.; Tortoli, P.; Mischi, M. Doppler ultrasound yechnology for fetal heart rate monitoring: A review. IEEE Trans. Ultrason. Ferroelectr. 2019 ,67, 226–238. [CrossRef] 36. Jezewski, J.; Roj, D.; Wr ó bel, J.; Horoba, K. A nowel technique for fetal heart rate estimation from Doppler ultrasound signal. Biomed. Eng. Online 2011,10, 1–17. [CrossRef] 37. Hamelmann, P.; Vullings, R.; Mischi, M.; Kolen, A.F.; Schmitt, L.; Bergmans, J.W.M. An Extended Kalman Filter for Fetal Heart Location Estimation During Doppler-Based Heart Rate Monitoring. IEEE Trans. Instrum. Meas. 2019,68, 3221–3231. [CrossRef] 38. Valderrama, C.E.; Marzbanrad, F.; Stroux, L.; Clifford, G.D. Template-based Quality Assessment of the Doppler Ultrasound Signal for Fetal Monitoring. Front. Physiol. 2017,8, 511. [CrossRef] [PubMed] 39. Divon, M.Y. Autocorrelation techniques in fetal monitoring. Am. J. Obstet. Gynecol. 1985 ,151, 2–6. [CrossRef] 40. Fukushima, T.; Flores, C.A. Limitations of autocorrelation in fetal heart rate monitoring. Am. J. Obstet. Gynecol. 1985,153, 685–692. [CrossRef] 41. Murrills, A.J.; Wilmshurst, T.H.; Wheeler, T. Antenatal measurement of beat-to-beat fetal heart rate variation: Accuracy of the Hewlett-Packard ultrasound autocorrelation technique. Fetal Physiol. Meas. 1986,6, 36–44. 42. Taylor, J.; Paull, C.J.; Hayes-Gill, B.R.; Crowe, J.A. Data compression of fetal Doppler ultrasound audio signals using zero-crossing analysis. Med. Eng. Phys. 1997,19, 572–580. [CrossRef] 43. Spilka, J.; Chudacek, V.; Bursa, M.; Zach, L.; Huptych, M.; Lhotska, L.; Janku, P.; Hruban, L. Stability of Variability Features Computed from Fetal Heart Rate with Artificially Infused Missing Data. Comput. Cardiol. 2012,39, 917–920. 44. Zhang, L.; Huang, M.J.; Wang, H.J. A Novel Technique for Fetal Heart Rate Estimation Based on Ensemble Learning. Mod. Appl. Sci. 2019,13, 137–147. [CrossRef] 45. Tuck, D.L. Improvement in Doppler ultrasound human foetal heart rate records by signal correlation. Med. Biol. Eng. Comput. 1982,20, 357–360. [CrossRef] 46. Romano, M.; Bifulco, P.; Ruffo, M.; Improta, G.; Clemente, F.; Cesarelli, M. Software for computerised analysis of cardiotocographic traces. Comput. Meth. Prog. Biomed. 2016,124, 121–137. [CrossRef] 47. Romano, M.; Iuppariello, L.; Ponsiglione, A.M.; Improta, G.; Bifulco, P.; Cesarelli, M. Frequency and Time Domain Analysis of Foetal Heart Rate Variability with Traditional Indexes: A Critical Survey. Comput. Math. Method. Med. 2016,2016, 9585431. [CrossRef] [PubMed] Sensors 2020,20, 4079 24 of 25 48. Wrobel, J.; Kupka, T.; Horoba, K.; Matonia, A.; Roj, D.; Jezewski, J. Recognition of fetal movements–automated detection from Doppler ultrasound signals compared to maternal perception. J. Med. Imaging Health Inform. 2015,5, 1319–1326. [CrossRef] 49. Czabanski, R.; Jezewski, J.; Wrobel, J.; Sikora, J.; Jezewski, M. Application of fuzzy inference system for classification of fetal heart rate tracings in relation to neonatal outcome. Gin. Pol. 2013 ,84, 38–43. [CrossRef] [PubMed] 50. Jezewski, M.; Czabanski, R.; Horoba, K.; Leski, J.M. Clustering with pairs of prototypes to support automated assessment of the fetal state. Appl. Artif. Intell. 2016,30, 572–589. [CrossRef] 51. Jezewski, M.; Leski, J.M.; Czabanski, R. Classification based on incremental fuzzy (1+p) -means clustering. In Man-Machine Interactions 4; Advances in Intelligent Systems and Computing; Gruca, A., Brachman, A., Kozielski, S., Czachorski, T., Eds.; Springer: Cham, Germany, 2016; Volume 391, pp. 563–572. ISBN 978-3-319-23437-3. 52. Importa, G.; Romano, M.; Ponsiglione, A.; Bifulco, P.; Faiella, G.; Cesarelli, M. Computerized Cardiotocography: A Software to Generate Synthetic Signals. J. Health Med. Inf. 2014,5, 1–6. 53. Mantel, R.; Ververs, I.; Colenbrander, G.J.; van Geijn, H.P. Automated antepartum baseline FHR determination and detection of accelerations and decelerations. In A Critical Appraisal of Fetal Surveillance; van Geijn, H.P., Copray, F.J.A., Eds.; Elsevier Science, B.V.: Amsterdam, The Netherlands, 1994; pp. 333–348. 54. Jezewski, J.; Pawlak, A.; Horoba, K.; Wrobel, J.; Czabanski, R.; Jezewski, M. Selected Design Issues of the Medical Cyber-Physical System for Telemonitoring Pregnancy at Home. Microprocess. Microsyst. 2016 ,46, 35–43. [CrossRef] 55. Wrobel, J.; Jezewski, J.; Horoba, K.; Pawlak, A.; Czabanski, R.; Jezewski, M.; Porwik, P. Medical cyber-physical system for home telecare of high-risk pregnancy–design challenges and requirements. J. Med. Imaging Health Inform. 2015,5, 1295–1301. [CrossRef] 56. Wrobel, J.; Matonia, A.; Horoba, K.; Jezewski, J.; Czabanski, R.; Pawlak, A.; Porwik, P. Pregnancy telemonitoring with smart control of algorithms for signal analysis. J. Med. Imaging Health Inform. 2015 ,5, 1302–1310. [CrossRef] 57. Sankhe, M.S.; Desai, K.D.; Gadam, M.A. Estimate of Fetal Autonomic State by Time Spectral and Nonlinear Analysis of Fetal Heart Rate Variability. Int. J. Comput. Inf. Syst. 2016,8, 312–325. 58. Kubo, T.; Inaba, J.; Shigemitsu, S.; Akatsuka, T. Fetal heart variability indices and the accuracy of variability measurements. Am. J. Perinat 1987,4, 179–186. [CrossRef] 59. Cesarelli, M.; Romano, M.; Bifulco, P. Comparison of short term variability indexes in cardiotocographic foetal monitoring. Comput. Biol. Med. 2009,39, 106–118. [CrossRef] [PubMed] 60. Goncalves, H.; Chaves, J.; Costa, A.; Ayres-de-Campos, D.; Bernardes, J. Comparison of the effect of different sampling modes on computer analysis of cardiotocograms. Comput. Biol. Med. 2015 ,64, 62–66. [CrossRef] [PubMed] 61. Goncalves, H.; Costa, A.; Ayres-de-Campos, D.; Costa-Santos, C.; Rocha, A.P.; Bernardes, J. Comparison of real beat-to-beat signals with commercially available 4 Hz sampling on the evaluation of foetal heart rate variability, available 4 Hz sampling on the evaluation of foetal heart rate variability. Med. Biol. Eng. Comput. 2013,51, 665–676. [CrossRef] [PubMed] 62. Jezewski, J.; Wr ó bel, J.; Horoba, K. Comparison of Doppler ultrasound and direct electrocardiography acquisition techniques for quantification of fetal heart variability. IEEE Trans. Biomed. Eng. 2006 ,53, 855–864. [CrossRef] [PubMed] 63. Cesarelli, M.; Romano, M.; Bifulco, P.; Fedele, F.; Bracale, M. An algorithm for the recovery of fetal heart rate series from CTG data. Comput. Biol. Med. 2007,37, 663–669. [CrossRef] [PubMed] 64. Jezewski, J.; Kupka, T.; Horoba, K. Extraction of Fetal Heart Rate Signal as Time Event Series from Evenly Sampled Data Acquired Using Doppler Ultrasound Technique. IEEE Trans. Biomed. Eng. 2008 ,55, 805–810. [CrossRef] 65. Kupka, T.; Matonia, A.; Jezewski, M.; Horoba, K.; Wrobel, J.; Jezewski, J. Coping with limitations of fetal monitoring instrumentation to improve heart rhythm variability assessment. Biocybern. Biomed. Eng. 2020 , 40, 388–403. [CrossRef] 66. Peters, C.H.L.; ten Broeke, E.D.; Andriessen, P.; Vermeulen, B.; Berendsen, R.C.M.; Wijn, P.F.F.; Oei, S.G. Beat-to-beat detection of fetal heart rate: Doppler ultrasound cardiotocography compared to direct ECG cardiotocography in time and frequency domain. Physiol. Meas. 2004,25, 585–593. [CrossRef] Sensors 2020,20, 4079 25 of 25 67. Roj, D.; Kupka, T.; Czaba´nski, R.; Pander, T.; Je˙zewski, J. Improvement in fetal heart periodicity measurement using Doppler ultrasound signal. In Proceedings of the 5th European Conference of the International Federation for Medical and Biological Engineering, Budapest, Hungary, 14–18 September 2011; Volume 37, pp. 133–136. 68. Al-Angari, H.M.; Kimura, Y.; Hadjileontiadis, J.; Khandoker, A.H. A Hybrid EMD-Kurtosis Method for Estimating Fetal Heart Rate from Continuous Doppler Signals. Front. Physiol. 2017,8, 641. [CrossRef] 69. van Geijn, H.P. Analysis of heart rate and beat-to-beat variability: Interval difference index. Am. J. Obstet. Gynecol. 1980,138, 246–252. [CrossRef] 70. Jezewski, M.; Czabanski, R.; Wrobel, J.; Horoba, K. Analysis of Extracted Cardiotocographic Signal Reatures to Improve Automated Prediction of Fetal Outcome. Biocybern. Biomed. Eng. 2010,30, 29–47. 71. Frigo, G.; Giorgi, G. Comparative Evaluation of On-Line Missing Data Regression Techniques in Intrapartum FHR Measurements. In Proceedings of the IEEE IMTC Conference, Torino, Italy, 22–25 May 2017. [CrossRef] 72. Lauersen, N.H.; Hochberg, H.M.; George, M.E.D. Evaluation of the accuracy of a new ultrasonic fetal heart rate monitor. Am. J. Obstet. Gynecol. 1976,125, 1125–1135. [CrossRef] 73. Dawes, G.S.; Visser, G.H.A.; Goodman, J.D.S.; Redman, C.W.G. Numerical analysis of the human fetal heart rate: The quality of ultrasound records. Am. J. Obstet. Gynecol. 1981,141, 43–52. [CrossRef] 74. Lawson, G.W.; Dawes, G.S.; Redman, C.W.G. A comparison of two fetal heart rate ultrasound detector systems. Am. J. Obstet. Gynecol. 1982,143, 840–842. [CrossRef] 75. Lawson, G.W.; Belcher, R.; Dawes, G.S.; Redman, C.W.G. A comparison of ultrasound (with autocorrelation) and direct electrocardiogram fetal heart rate detector systems. Am. J. Obstet. Gynecol. 1983 ,147, 721–722. [CrossRef] © 2020 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 (http://creativecommons.org/licenses/by/4.0/).