Frequency and Time Domain Analysis of Foetal Heart Rate Variability with Traditional Indexes: A Critical Survey

Monitoring of foetal heart rate and its variability (FHRV) covers an important role in assessing health of foetus. Many analysis methods have been used to get quantitative measures of FHRV. FHRV has been studied in time and in frequency domain and interesting clinical results have been obtained. Nevertheless, a standardized definition of FHRV and a precise methodology to be used for its evaluation are lacking. We carried out a literature overview about both frequency domain analysis (FDA) and time domain analysis (TDA). Then, by using simulated FHR signals, we defined the methodology for FDA. Further, employing more than 400 real FHR signals, we analysed some of the most common indexes, Short Term Variability for TDA and power content of the spectrum bands and sympathovagal balance for FDA, and evaluated their ranges of values, which in many cases are a novelty. Finally, we verified the relationship between these indexes and two important parameters: week of gestation, indicator of foetal growth, and foetal state, classified as active or at rest. Our results indicate that, according to literature, it is necessary to standardize the procedure for FHRV evaluation and to consider week of gestation and foetal state before FHR analysis.


Introduction
Foetal heart rate (FHR) monitoring is of great importance to obtain information about foetal health during pregnancy and labour. In particular, electronic monitoring of the FHR (EFM), commonly named cardiotocography when an external Doppler probe is used, is the most employed method to detect foetal distress and prevent neurologic damage or even foetal death [1,2]. However, despite its usefulness for obstetricians, the problem of EFM is its poor predictive value. It often lacks specificity leading to unnecessary interventions that increase caesarean delivery and operative vaginal delivery rates [1][2][3][4]. Moreover, there is often disagreement between obstetricians analysing FHR traces, since the interpretation is usually performed by means of a visual inspection, which obviously lacks objectivity and reproducibility [5,6]. In recent years, hence, interest has grown in how to recognize changes in FHR that might predict more accurately foetal distress. For example, in order to overcome the subjective nature of FHR interpretation, several attempts have been made to automate the diagnosis of the foetal status and many computerised algorithms have been developed to assess FHR parameters [1,[6][7][8][9][10].
Independently of the recording methods, the main FHR morphologic characteristics and parameters observed by physicians for foetal health evaluation are FHR mean value (which is related to week of gestation), baseline of the FHR, acceleration's rate and shape, deceleration's rate and shape, and FHR variability (FHRV) [11][12][13][14][15]. Among these, FHRV is probably the most important one, since it reflects the activity of autonomic nervous system (ANS) in the foetus who is growing and developing [1,[16][17][18][19], although the exact contributions of the two branches of the ANS are still object of investigation even in adult subjects [20]. The study of FHRV, often referred to as the beat-to-beat fluctuations of the FHR signal, could be, like for heart rate variability (HRV) of adult subjects, a base for a more powerful, detailed, and objective FHR analysis and for better knowledge of ANS reactions and its development [11,16,[21][22][23][24]. With respect to the first studies [25,26], nowadays, the knowledge of FHRV is improved, so that different ranges of variability can be identified in order to classify FHR recordings [27] and its assessment employed to evaluate foetal reactivity and wellbeing in nonstress condition [28]. Even though the presence of good variability may not always be, by itself, a certain sign of reassuring FHR signal (corresponding to a well-oxygenated foetus), most clinicians agree that minimal or absent variability could be an indicator of foetal distress [28][29][30].
Due to its recognized importance, a large number of new analysis methods have been enforced to obtain more objective and quantitative measures of FHRV [31,32]. Traditionally, as HRV of adult subjects, FHRV can be studied both in time (statistical indexes) and in frequency (spectral indexes) domain. In the time domain, Short Term Variability (STV) indexes and Long Term Variability (LTV) indexes are usually distinguished. In the frequency domain, different methods have been employed to estimate the power spectral density (PSD), which is widely considered the index that best covers all the information of the heart rate series. However, since FHRV signal shows a nonstationary behaviour, the time-frequency analysis of the FHRV is generally employed [1,22,[33][34][35][36][37][38][39]. As far as time domain analysis (TDA) is concerned, it could present some limitations because it mainly relies on statistical measurements, so it can only describe the magnitude of the variability around an average value, without providing further information about the physiological mechanisms involved [40,41]. Some limitations have been also shown for the frequency domain analysis (FDA) since it is generally sensitive to artifacts [42] and can provide information only about periodical fluctuations of the heart rate rhythm, without inspecting other possible nonperiodic trends embedded in the variability signal [40,41]. In order to overcome these limitations and also to investigate and improve risk stratification, during the last decades, techniques analysing nonlinear dynamics, such as symbolic dynamics, approximate entropy, and fractal analysis, have been employed both in adults and in foetuses, even if, at the moment, none appears to be predominant and completely satisfying [26,43,44].
TDA and FDA hence remain the most used methods of HRV and FHRV analysis also due to their simplicity and higher acceptance in clinical environments. Nevertheless, a standardized definition of FHRV is still lacking and subject to changes and updates. Traditional analysis lack of standardized methods for computing time or frequency domain indexes and the relationship between FHRV and foetal growth is actually not completely clear and most of the studies inspecting this issue are by now rather old [52,54,71,72].
This works aims to compare some common indexes of TDA and FDA, employing both real and artificial FHR signals.
To get this objective, we firstly introduce a brief report on the most relevant literature works about traditional FHRV analysis, focusing the attention on those studies based on the computation of the STV, as a time domain index, and on the power of the FHRV in different frequency bands.
Then, through the use of simulated FHR signals, we define the methodology to be employed for the estimation of the chosen indexes for FDA.
Finally, we apply the time and frequency indexes to the analysis of real FHR recordings in order to assess the capability of these indexes to correlate with foetal development during gestation and with foetal state and to provide an overview of their reference values.

Spectral Analysis of FHRV:
A Brief Literature Report. The study of biologic signals in the frequency domain can offer deeper knowledge of their behaviour. In adults and foetuses, the spectral analysis permits estimating the power of the periodic HR fluctuations and it represents a noninvasive and powerful tool to understand ANS functional state and reactions [28].
After the study of Akselrod et al. [73], which introduced the power spectral analysis of short term heart rate fluctuations as a noninvasive quantitative probe of beat-tobeat cardiovascular autonomic control, FDA has been widely performed to point out the relation between the ANS activity and low frequency and high frequency bands, whose power content reflects changes of sympathetic and vagal activity [24,26,61]. Besides, changes in power distribution have been recognized as predictors of foetal distress, both in antepartum and in intrapartum periods [65,74]. The work of Padhye et al. [19] investigated the correlation of power in LF (computed using the Lomb periodogram) and HF bands and noted an increasing trend of the power with gestational age. In their study, van Laar et al. [63] used fast Fourier transform to calculate the FHRV spectral power in LF and HF bands in order to compare spectral values between near term and postterm foetuses and found a sympathetic predominance in foetuses near term during active state, but an increased vagal modulation in postterm foetuses during rest state. According to Kwon et al. [4], changes in spectral power corresponding to a low pH are different between term and preterm foetuses, confirming a correlation between frequency indexes and gestational age.
Despite these interesting and important clinical results, problems in interpretation and comparison arise because literature works employ different frequency analysis methodologies, use discordant measure units for the PSD, and show disagreement about the frequency bands of the FHRV spectrum [36,75], even if most of the literature agrees that, like the case for adult subjects, three bands can be detected in the FHR power spectrum: a very low frequency (VLF)band, which seems to be related to thermoregulation mechanisms [12,67]; a low frequency (LF) band, which is mainly associated with the sympathetic branch activity and is an indicator of foetal development and wellbeing [54,67]; a high frequency (HF) band, which reflects the respiratory activity and the vagal stimulation [38,45,47,49,61].
In order to clarify the definition of frequencies bands, a study of literature was conducted, involving about eight Results from the studied literature works (listed in chronological order) on FHRV spectral bands (l and u indicate, resp., the lower and upper limit of each band).
hundred literature works concerning foetal monitoring, published between 1983 and 2013. Among these, only a hundred works are directly related to frequency analysis and only about thirty works gave details about the three bands; the major disagreement is about the VLF band. Some researchers identifies the VLF band in the range from 0 to 0.03 Hz, while others consider VLF band ranging from 0 to 0.04 Hz or from 0 to 0.05 Hz [12,50,56,67]. Furthermore, some authors introduce a middle frequency (MF) band in the range of 0.15-0.5 Hz or 0.2-0.4 Hz [36,43]. A concise overview of different literature works focused on the computation of the FHRV spectral bands and corresponding power content is shown, respectively, in Tables 1 and  2. Empty cells are due to the absence of data in the original papers.

STV in Foetal
Monitoring: A Brief Literature Report. As mentioned, in the time domain, Short Term Variability (STV) indexes and Long Term Variability (LTV) indexes are usually distinguished. The former, also according to FIGO guidelines [76], refer to the continuous small changes in difference between successive interbeat intervals, which occur under physiological conditions. These minimal oscillations cannot be reliably interpreted by the naked eye; furthermore, a correspondent shared mathematical definition is lacking, so that this important parameter has lost part of its relevance and in some more recent guidelines it is not even considered but they speak broadly of variability, referring implicitly to the amplitude of FHR signal and without differentiating by LTV, since in practice they are visually determined as a unit [14,77,78]. When a computerised system is available, different indexes are used, many of which are borrowed from studies concerning adult heart rate. Among them there are Root Mean Square Successive Difference (RMSSD), that is, the Results from the studied literature works (listed in chronological order) on FHRV power estimation.
square root of the mean squared differences of successive RR intervals, and pNN50, that is, the percentage of differences between following RR intervals greater than 50 ms [26,79]. The LTV, instead, refer to fluctuations in the FHR over seconds, such as SDNN-Index, that is, the mean of the 5minute standard deviation of the NN interval (normal to normal interval) calculated over 24 h, and SDANN, that is, the standard deviation of the average NN interval calculated over short periods [26].
In clinical practice, beat-to-beat indexes are often preferred since a good beat-to-beat variability is widely accepted as a significant index to assess foetal wellbeing, since a good beat-to-beat variability is a reliable indicator of a healthy foetal ANS [80]. Hence, many studies in time domain attempted to compute indexes for quantifying STV in foetuses by using very different techniques and methods (modification of the mean, standard deviation (SD), slope changes, and varying epoch lengths) [81]. The lack of a unique standardized methodology along with the fact that STV formulas, usually based on ECG, are often applied without any adaptation to the ultrasound technique makes the comprehension of the measure and the comparison between two or more indexes very difficult [79].
Despite the lack of standardization, STV is broadly employed. For example, Short Term Variability was found to be a good predictor of Apgar scores by Ayres-De-Campos et al. [82]. D'Elia et al. [83] analysed healthy term foetuses subjected to vibroacoustic stimulation by means of computerised CTG and found a statistically significant increase in foetal movements, acceleration rates, and STV with foetal activity. In their work, Serra et al. [84] examined the clinical value of the STV in the timing of the delivery of severely growthretarded foetuses and confirmed that the STV can assess the condition of foetuses with severe intrauterine growth restriction (IUGR) and that it is an important marker of perinatal outcome in severely growth-retarded foetuses. Also Galazios et al. [85] have rather recently observed that STV value is associated with foetal distress and, more recently, the study of Annunziata et al. [86] evaluated the impact of vibroacoustic stimulation on STV of CTG recordings in low and high risk pregnancies and noted that an increase in STV is significantly associated with good perinatal outcome.
Finally, Cesarelli et al. [79] proposed a comprehensive study on nine different mathematical indexes utilised to compute STV from CTG recordings, testing their robustness, sensitivity, and dependence on other parameters (FHR storage rate, FHR mean value, etc.) and demonstrated that the SD index, computed after floating line extraction, provides efficient information and is independent of the considered variables.

Data Collection.
Real CTG traces were recorded by healthy pregnant women during the clinical practice, using commercially available cardiotocographs (HP-135x or Sonicaid). Five hundred and eighty recordings, lasting on average more than 25 minutes, recorded from women between the 24th and the 42nd gestation weeks, were considered for the study. Gestational age was determined from the last menstruation date or from ultrasound measurements executed in the Computational and Mathematical Methods in Medicine 5 first trimester of pregnancy. The database was completed with other pieces of clinical information of patients and newborns. All patients gave their informed consent to participate in the research concerning foetal monitoring. CTG signals were processed by software previously developed by the authors [7,84]. Firstly, they were preprocessed by means of a software [10,84] which processes signals in output from the cardiotocographs in order to recognize signal tracts having good and bad quality (these last including tracts of signal loss); for each segment of good quality, recover the real uneven FHR series when CTG output is evenly spaced (case of HP/Philips cardiotocographs) [87] and detect and process outliers [88]; interpolate signal tracts of poor quality (according to an index provided by the equipment) or signal loss which last maximum 3 s, in order to avoid an excessive fragmentation of the signal.
In this way, all CTG signals available for further processing have the same characteristics and are unevenly sampled (in correspondence with the real heart beat) regardless of the equipment employed for their recording, hence regardless of the acquisition and sampling mode.

CTG Simulation.
Like it is known, real FHR signals are affected by a considerable amount of variability and complexity, because of relationships, complicated and not yet fully known, among the different physiological mechanisms involved in heart rate regulation. Hence, in order to have available signals with characteristics known a priori, artificial CTG signals were employed for the analysis carried out in order to define an adequate methodology for FHRV evaluation.
According to previous works [79,89,90], an artificial uneven RR series with specific power spectrum characteristics, proper for FHR, was firstly generated. In particular, we set central frequencies of the spectral bands, the bandwidths, and the ratio between power of low and high frequency bands. Then, in line with the operations made by the cardiotocographs which detect heart beats (by means of Doppler technique), compute interbeats distances series, and then reversing it get the FHR signal, this last one was computed as FHR = 60/RR, fixing FHR mean value and standard deviation of its peak-to-peak amplitude. After that, acceleration rates and deceleration rates with different amplitude and duration were simulated by using Gaussian-like signal tracts (for a more detailed explanation about how these synthetic signals were generated, please refer to previous works of the authors [79,89,90] and see Figure 1 for an example).
Our final simulated signals, hence, are in principle similar to those obtained with ultrasound technique [90].
The software for generating artificial CTG signals, as well as that for signal processing employed for his work, was developed using MATLAB R2011a [10].

FHRV Estimation.
Previously, we proposed definition of FHRV as difference between FHR signal and its floating line [10,12]. Of course, for FHRV assessment, the floating line has to be correctly estimated, so that we proposed also a procedure, using spline nonlinear filtering, developed to this aim (the methodology has been recently updated, and results are submitted but not yet published).
Here, as further test, we compared the mean frequency spectrum computed on 30 simulated FHRV signals obtained with two different methodologies: by means of the application of our procedure for floating line estimation or simply as a result of the detrend operation, often employed in the literature [16,55,59,63,67]. (Let us remember that detrend is an operation which removes the linear trend from a signal; in Matlab it is a default function.)

Frequency Domain Indexes.
We estimated FHRV as explained in Section 3.3 and considered for this signal LF (0.05-0.2 Hz) and HF bands (0.2-1 Hz); then, as VLF band, we computed the power spectral density (PSD) of the floating line.
Because of the nonstationarity of FHR and hence of FHRV, PSD was estimated by means of the Short Time Fourier Transformation [10], the methodology still more used for its simplicity, using a sliding Hamming window of 32 s [34,47]. This window is shifted sample by sample and a new PSD is computed each time [35]. To be able to use the STFT, the FHR signal was previously interpolated (4 Hz sampling rate) by means of cubic interpolation, which has been demonstrated to reduce the error introduced [37]. Then, by means of a simple integral rule, we have computed the power values, absolute (called LF , HF ), and percentage, with respect to the total power (called LF%, HF%), in the bands defined above. Further, the sympathovagal balance (SVB) index, which is an important index that reflects the relations between vagal and sympathetic branches of the ANS [26], was computed as ratio between LF and HF . Finally, using the same methodology, we computed VLF (and VLF%) as power of PSD of the floating line (using a commercial personal computer with a processor i5-3337U@1.8 GHz and RAM of 8 GB, the computation requires about 4 s for a CTG signal of more than 2 hours of length and about 10 s when also a 3D representation of the PSD is required).

Time Domain Indexes.
As index of variability in time domain, we chose STV for its importance in monitoring of foetal health, since it is related to regulation mechanisms elicited by ANS activity [84,86], as already mentioned in the previous sections. In agreement with a previous study of the authors, we assessed STV as standard deviation of the FHRV, obtained after subtraction of the floating line from the FHR signal [79]. In order to enrich the available information, without worsening the computational complexity, the calculation is carried out on sliding windows of length (with covering 30 s), with an overlap of − 1 samples [10]; then, to provide an overall STV index of the signal, the mean of all STV values is computed (the computation and the graphical representation require about 8 s for a CTG signal of more than 2 hours of length, using the same personal computer utilised for FDA).

Week of Gestation and Foetal State.
To test time and frequency indexes, we chose to verify their relationship with the week of gestation, the most simple indicator of foetal growth, and the foetal state.
It is known that foetal behavioural states have a great importance in influencing FHR patterns. However, the exact recognition of the foetal state is not a simple task, since respiratory acts, eyes closure and opening, just to name some aspect, should be simultaneously detected by ultrasound imaging. Hence, according to the literature [36,63], for CTG recorded from 30th week of gestation onward, we defined active or resting foetal state. In particular, we classified a FHR signal as recorded from a foetus in active state if, at visual inspection, it showed a good variability (at least equal to 5 bpm) and at least two acceleration rates in 20 minutes of recording. If the signal showed low variability and absence of acceleration rates, the foetus was classified as at rest. Doubt cases were excluded from the analysis.
In Figure 2, examples of FHR signals classified as active and rest are shown ((a) and (b), resp.).

Statistical Analysis.
For all indexes studied, we analysed the correlation with week of gestation by means of regression lines (a polynomial quadratic curve was considered).
Besides, we tested the ability of these indexes to differentiate active from rest foetal state by means of the Mann-Whitney test, because of the non-Gaussian distribution of data.

FHRV Mean Spectrum.
In order to verify which is the more appropriate operation for FHRV estimation, detrend or floating line subtraction, in Figures 3(a) and 3(b), we show results of the frequency analysis conducted on simulated FHR signals.
It is clear that when the detrend operation is used (Figure 3(a)), a residual of VLF component appears superimposed on LF component. Of course, that can alter power computation and results of sympathovagal balance estimation. Table 3, the results of the Mann-Whitney test carried out for all tested indexes are shown; they concern Week of gestation the comparison between foetuses at rest (55 CTG recordings) and active foetuses (384 CTG).

FHRV Indexes.
Since all indexes, except SVB, resulted significantly differently in active state with respect to rest state, in Tables 4 and 5 we report ranges (minimum and maximum value) of all indexes here analysed and their mean and SD, computed on real FHR, separately for the two foetal states. Figures 4-7, as examples, the trends with the pregnancy course of average values of some analysed indexes are depicted and in Table 6 values of coefficient of determination ( 2 ) are shown for all indexes here analysed. In this case, the analysis is carried out without distinguishing foetal states that are not defined early in pregnancy. Week of gestation FD-P VLF Week of gestation

Conclusions
In this paper, firstly we provided an overview of the literature concerning the use of some traditional indexes. It is a concise overview but, at the best of our knowledge, it is the first   time that so many quantitative indications and results are compared (Tables 1 and 2).
In the review section, we reported just some clinical results (being this aspect out of the main aim of the work and already treated in the literature [75]) and we neglected the choice of the methodology to compute PSD, since in a previous work we compared three methodologies and obtained not so different results. So, we focused on the definition of FHRV (since a FHRV definition shared and mathematically translatable is still missing) and frequency bands.
As here shown, about the frequency spectrum, although there is agreement in considering three main bands, VLF, LF and HF, and a MF when a more detailed analysis is necessary, even the same author or authors of the same research group employ different bands limits in their works (Table 1). Besides, frequency analysis is often carried out without providing power values and, in case, a large variability is present among different papers ( Table 2), so that a comparison is quite difficult.
Furthermore, it is worth underlying that the analyses are very often carried out starting from RR intervals, maybe for historical reasons or for "continuity" with studies involving adult subjects. Nevertheless, ECG is not yet so diffused in clinical routine for difficulty in signal processing (in case of recording through maternal abdomen) or ethical reasons (in case of direct fECG), whereas clinicians are used to analyse FHR signal, measured in beats per minute (bpm).
For these reasons, we prefer to process FHR signal (without converting it into RR signal) and propose definition of FHRV as the difference between FHR and floating line. Preliminary results here shown (Figure 3) confirmed the usefulness of this methodology. About PSD computation, in previous works [34,37], we used the Lomb method which can be applied directly on uneven series; however, it has too long computational time by losing an advantage of the frequency analysis that can be used, if desired, also for real-time analysis; therefore, for successive as well as for this work we employed STFT. Furthermore, comparing different papers and taking into account experimental results obtained on our database, we decided to define for FHR signals the following bands limits: 0-0.05 Hz for VLF; 0.05-0.2 Hz for LF; and 0.2-1 Hz for HF.
Once defined details of the procedure are to be employed, through a retrospective study on a very large amount of CTG data recorded during physiological pregnancies (439 signals whereas most papers report results computed on a much smaller number of signals), we provided range of values of all frequency indexes here considered (differentiated on the basis of the foetal state, Tables 4 and 5). Some of them show high values of SD that put in evidence the high intersubject variability, even in healthy foetuses.
About STV index, the literature is even in more disagreement and we did not find reference values to report in our brief summary. However, we limited our analysis to a literature overview of the main results obtained in its clinical use since previously we already analysed different mathematical formula employed for its assessment and proposed a new evaluation methodology. By means of this methodology, we computed its value in healthy foetuses and for week of gestation. It, on average, resulted in increase with gestational age (Figure 4), according to the literature [91], and is very different in rest or active foetal state.
Then, we carried out a statistical analysis to test the correlation between the different indexes employed and the foetal growth (week of gestation).
Regression analysis showed that, on average, only some indexes can represent foetal growth in a satisfying way ( 2 about 0.7, Table 6). Moreover, although a reliable comparison is rather difficult, because of the significant differences in employed methodologies and in a general lack of details provided, as proved here and in the literature [75], they are Computational and Mathematical Methods in Medicine 9 substantially consistent with others. Absolute LF power and HF power, for example, (Figure 5), increase with gestational age but decrease their growth rate towards the end of pregnancy [58,61]. With regard to the VLF, we do not have clear literature results with which we make a comparison; however, we may observe that the power in this band increases with week of gestation ( Figure 6) and this is coherent with the increase of foetal movements and acceleration rates [61,92] (let us remember that we computed VLF as power of the floating line which, in turn, includes acceleration rates and deceleration rates). Besides, its value represents the most of the total power so that the two indexes have an equivalent behaviour both with gestational age and with regard to foetal state. We can observe also that in the last weeks of gestation all powers change really little (both in absolute and percentage values), confirming the almost ended development of the SNA (further modifications will concern the adaptation to postnatal life).
The trend of SVB ( Figure 7) is an exception but it is not a surprising result. Its value decreases with gestational age, coherently with the literature [36,56,61], and with the increase of HF band, related to development of the vagal branch.
Finally, about the foetal state, almost all indexes appeared to be able to separate the two groups ( Table 3).
The findings we get in this work lead us to say, in accordance with the literature [63,75], that standardization of FHRV assessment is necessary and that, since foetal state and gestational age can strongly affect results, it is not possible to process FHR signals regardless of these conditions.