Gap detection, interpolation, outlier detection, and detrending are frequently utilized in the preprocessing stage of FHR analysis to achieve reliable signals (Cesarelli et al., 2007; Romano et al., 2013).

Feature transforms (also called feature engineering in the ML field) have great importance for signal representation. The basic morphological features such as baseline (Dawes et al., 1982), the number of acceleration and deceleration patterns (Mantel et al., 1990) and variability in the short-term and long-term (Cesarelli et al., 2009) are fundamental parts of the computerized FHR analysis. Additionally, the linear and statistical features coming from the time-domain and frequency-domain are extracted to support the automated analysis (Czabanski et al., 2012; Dash et al., 2014; Spilka et al., 2014). Further, using non-linear parameters (e.g., entropy, complexity, and fractal dimension) in fetal state assessment have been proposed and tested (Spilka et al., 2012; Fanelli et al., 2013; Doret et al., 2015; Comert and Kocamaz, 2016). Recently, an image-based time-frequency (IBTF) feature analysis approach comprised of a combination of short term Fourier transform (STFT) and a gray level co-occurrence matrix (GLCM) have been employed as diagnostic indices for fetal hypoxia detection (Comert et al., 2018). Moreover, transform-based analysis methods, such as a discrete wavelet transform (DWT) (Chrelias et al., 2008) and empirical mode decomposition (EMD) (Krupa et al., 2011), have been used for both noise reduction and feature extraction. Especially, heart rate variability (HRV) has become a valuable resource in FHR signal characterization (Romano et al., 2018). On the other hand, since not all extracted features are valuable for classification, feature selection algorithms and dimensional reduction methods have been utilized to select an optimal feature set to improve the performance, including a genetic algorithm (GA) (Xu et al., 2014), information gain (IG) (Spilka et al., 2012), and principal component analysis (PCA) (Zhang and Zhao, 2017).

Last, after signal preprocessing and feature transformation, the computerized systems employed ML algorithms to perform two or multiclass classification tasks and thereby discriminate a pathological fetus from a normal fetus. For example, Czabanski et al. (2012) designed an expert system to predict neonatal pathology using a two-stage analysis method based on weighted fuzzy scoring (WFS) and least square support machine (LS-SVM) and obtained good performance with an accuracy of 92.0%. Comert and Kocamaz (2016) applied an artificial neural network (ANN) and performed a classification task with an accuracy, sensitivity, and specificity of 92.40, 95.89, and 74.75%, respectively.

Clearly, although the previous studies concerning computerized FHR analysis presented good classification performance in assessing the fetal state with an accuracy of 90–95%, the conventional ML method needs to extract informative features and select optimal features, which can then be fed into classifiers. Therefore, this approach requires a heavy workload and detailed physiological information regarding the fetus that may be lost during the entire procedure.

Deep learning (DL) has become a highly useful tool for image processing in recent years (LeCun et al., 2015). Especially, convolutional neural networks (CNNs), which involve many layers, have been found to be quite efficient for most image classification problems (Krizhevsky et al., 2012). A CNN essentially self-learns high-level informative features spontaneously by constructing hidden multilayer networks and training data to improve performance in biomedical data processing (Bursa and Lhotska, 2017). Thus, CNN greatly simplifies the subjectivity of feature extraction and discovers the intricate patterns contained in input data compared to manual feature engineering in traditional ML methods. Due to the attractive advantages, CNNs are extensively utilized in the medical field for the purpose of designing diagnostic tools that automatically assist clinicians (Archarya et al., 2017; Comert and Kocamaz, 2018; Li et al., 2018; Petrozziello et al., 2018). For example, Archarya et al. (2017) designed the CNN structure to diagnosis coronary artery disease using an electrocardiogram (ECG) signal and achieved high accuracy of 95.11% (Archarya et al., 2017). In addition, Comert and Kocamaz (2018) also proposed a novel approach to detect fetal hypoxia based on a deep CNN with transfer learning using the FHR signal and Petrozziello et al. (2018) reported their investigation of long short term memory (LSTM) and CNN in analyzing continuous EFM traces from over 35,000 labors for the prediction of fetal compromise. In this study, motivated by previous research, an end-to-end classification strategy (which means that feature transform procedures are ignored) is performed to determine the fetal state using a CNN algorithm, which possesses the ability to self-learn useful features from the input FHR signals.

The remainder of the paper is organized as follows: Section “Materials and Methods” describes the database and overall methodology of the proposed system; Section “Results” presents the experimental setup, results, and relevant discussion; and last, the study’s conclusions are presented and future research is outlined in Section “Discussion.”

A publicly accessible intrapartum Czech Technical University - University Hospital in Brno (CTU-UHB) CTG database was used in the experiment. Chudacek et al. (2014) collected the whole signals between 2010 and 2012 using STAN S21/S31 and Avalon FM40/50 EFM devices. All signals were sampled at 4 Hz and stored in electronic form in an OB TraceVue^® system provided by Philips. Furthermore, the authors selected a total of 552 intrapartum CTG recordings from a subset of 9164 recordings to constitute this database after considering several technical and clinical criteria, such as the woman’s age, week of gestation, type of gravidity, type of delivery, signal quality, and labor outcome measures. Table 1 shows the main parameters and their respective distributions of this database. Due to a space restriction, interested readers can refer to more detailed information about the database in Chudacek et al. (2014). The maximum duration of the recordings was 90 min and each recording started at a maximum of 90 min before delivery. The database is open access and can be freely downloaded from Physionet (Goldberger et al., 2000).

All signals were visually interpreted as being divided into four parts by nine experienced obstetricians (Hruban et al., 2015). Moreover, subjective evaluation criteria, Apgar’s scores, were provided for the 1st and 5th min. Conversely, additional biochemical markers measured after delivery, such as the umbilical artery pH, the base excess (BE), and the base deficit in extracellular fluid (BDecf), were provided for an objective categorization. In this work, to make the comparison among different methods more effective (Spilka et al., 2012; Dash et al., 2014; Hruban et al., 2015; Comert et al., 2018), we adopted the suggestion of these studies and the pH value was adjusted to 7.15, which served as a borderline to separate the fetal state into two classes after comprehensive consideration. The signals that had a greater or equal pH than the specified value were assessed as normal. As a result, 447 normal and 105 abnormal (hypoxia) FHR signals were obtained. To avoid the class imbalance problem, we further randomly selected 105 of 447 normal fetuses to keep the number of abnormal fetuses equal.

Preprocessing is an important step in biomedical applications and can affect not only the subsequently extracted features but also the final classification performance. The FHR signal has two typical acquisition methods: the Doppler ultrasound probe placed on the abdomen of a pregnant woman externally and the direct ECG measured on an electrode connected to the fetal scalp internally (Kunzel, 2012). From this point of view, the FHR may be contaminated by many aspects of noise, such as the movement of pregnant women and fetuses, improper placement of sensors and other equipment, and external environmental factors. For the sake of content integrity, we briefly describe the preprocessing algorithm used in this work.

There are two manifestations of noise contained in the FHR signal: missing values and spiky artifacts (Cesarelli et al., 2007). For the former, a signal segment (FHR is equal to 0) that lasts longer than 15 s is removed directly; otherwise, it is linearly interpolated. Then, for the FHR signal that is not stable (the absolute value of two adjacent points is greater than 25 bpm), interpolation is again performed between the initial sampling point and the first point of the next stable part. Finally, the extreme points are also removed; that is, the signal value is greater than 200 bpm or less than 50 bpm, and the segment is then filled in with Hermite spline interpolation. Overall, the preprocessing algorithm can be summarized as having three steps: gap detection, interpolation, and outlier detection (Romano et al., 2013). Figure 2 shows the original noisy signal and preprocessed signal to be further analyzed (10 min in length).

As a non-stationary and non-linear time series, phase space reconstruction is the first and most important step in the analysis of the FHR signal based on the dynamics theorem. Packard et al. proposed the method of utilizing the time delay coordinates to reconstruct the phase space trajectory of the signal (Packard et al., 1980). Then, Takens (1981) also presented an approach for reconstructing a non-linear system which only requires a time-ordered sequence. Based on Equation (1), a scalar 1-dimensional time series u_i (i = 1, 2, 3, …, L) is usually embedded into an m-dimensional space using this method involving time delay.

where the vector x_k represents the k-th point on the orbit (k = 1, 2, …, L-(m-1)τ). The τ is the delay time, m is the embedding dimension and m is ≥2, which are both difficult to be determined.

Furthermore, there is a prior condition in Takens’ theory that assumes that the signal is infinite and does not contain noise, where the delay time τ can be chosen almost arbitrarily (Packard et al., 1980). Unfortunately, the real FHR signals recorded by clinical equipment are mostly finite and noisy; therefore, we needed to consider the choice of the delay time carefully. If τ is too large, it may incur irrelevance, where the dynamics of an attractor are independent and become causally disconnected. In contrast, if τ is too small, it may lead to redundancy, where the reconstructed attractor is compressed along the identity line (Kim et al., 1999).

Similar to the delay time, the embedding dimension vitally needs to be determined. If m is too small, the geometry of the phase space is partly folded. If m is too large, it may result in massive calculations and increasing the contamination because of the rounding or instrumental error (Kennel et al., 1992). The dimensions of the various portions of the non-linear signal may differ from each other and tend to increase the required dimension influenced by noise. Thus, the choice of a large enough dimension that consists of the relevant dynamics for a noisy signal is considerable (Fraser and Swinney, 1986).

Based on the above-mentioned phase space reconstruction, RP is a qualitative analysis approach for biomedical signals. The concept of RP was proposed by Eckmann et al. (1987), which can depict the recurrence property of a deterministic dynamical system; i.e., visualizing the time dependent behavior of orbits x_i in a phase space. The RP can intuitively reflect the high dimensional phase space motion law of non-linear signals in 2- dimensional space. The key step of an RP is to define the distance between any two vectors in the phase space and its mathematical expression is:

where 𝜀 is a predefined distance threshold; x_i and x_j represent the i-th and j-th phase space vectors, respectively, which can be reconstructed by using the Takens’ time delay approach mentioned above; N is the total number of phase points and N = L-(m-1)τ, ||…|| represents the norm (e.g., Euclidean distance), and Θ(x) is a Heaviside function.

We can generally explain Equation (2) as follows: if x_j is sufficiently close to x_i, which means that x_j falls within a ball of the cutoff distance 𝜀 centered at x_i, R_i,j is 1, then a white dot is placed at a point (i, j); otherwise, R_i,j is 0 and the dot is black. Then, the binary values of R_i,j can be simply denoted as the white (1) and black (0) points, respectively. Hence, we can consider the RP as a visual inspection of a high-dimensional phase space trajectory (Eckmann et al., 2008). An N × N distance matrix can be converted into a 0-1 matrix and this allows the recursive property of the time series graphically. The 𝜀 is chosen by k’s nearest neighbors method in this work (Ouyang et al., 2014).

In summary, there are three optional parameters of RP that can be adjusted: m, τ, and k, which indicate the embedding dimension, delay time, and distance threshold, respectively. Figure 3 shows the RPs of normal and pathological fetuses.

As a brand-new field that has been developing rapidly for more than a decade, DL has received increasing attention from researchers, and it has obvious advantages over shallow models (i.e., ML) in feature engineering and model building (Zhang and Zhao, 2017). DL is good for excavating an increasing number of abstract features that possess better generalization ability from the original input data, especially the image (Schmidhuber, 2015). It overcomes some of the problems that were thought to be difficult to solve in artificial intelligence (AI) in the past. Moreover, with the significant increase in the number of training datasets and the dramatic increase in chip computational power, DL has achieved remarkable success in object detection, computer vision, natural language processing, and voice recognition, so it has also promoted the development of AI (Deng and Yu, 2014; Bengio, 2015). DL is a hierarchical ML method that includes multilevel non-linear transformations and the deep neural network (DNN) is currently the main form (Canziani et al., 2017). The connection mode between neurons in the DNN is inspired by animal visual cortex organization, and CNN is one of the classical and widely used structures (Bouvrie, 2006). The local connections, weight sharing and pooling operations of the CNN can significantly reduce the complexity of the network. Additionally, the number of training parameters makes the model have a certain degree of invariability to translation, distortion, and scaling (Krizhevsky et al., 2012). In addition, CNN has strong robustness and fault tolerance, and it is also easy to train and optimize. Based on these superior characteristics, CNN effectively outperforms the standard fully connected neural network in various signal and information processing tasks (Bursa and Lhotska, 2017).

The basic structure of the CNN is composed of an input layer, a convolution layer, a pooling layer (also called the subsampling layer), a fully connected layer, and an output layer (Bouvrie, 2006). Convolution and pooling layers generally involve several alternate convolution and pooling layers; that is, one convolution layer is connected to one pooling layer, a pooling layer is followed by one convolution layer, and so on (Liu et al., 2015). Because each neuron of the output feature map in the convolution layer is locally connected to its input, the weighted sum of the corresponding connection weights and local input plus the bias value determine the input value of the neuron. The principle is equivalent to the convolution process in mathematical meaning, and it was therefore given the name CNN. Figure 4 shows a graphical representation of our designed CNN architecture consisting of an image input layer and two convolutions; that is, normalization, ReLU layers, two average pooling layers, two fully connected layers, a dropout layer and a final softmax output layer.

After the convolution layer, a batch normalization layer is useful in reducing sensitivity to variations with the input data and an activation layer is usually applied to introduce the non-linearities into the network through mapping to the input data (Krizhevsky et al., 2012). Compared with traditional sigmoid and tanh functions, we employed the rectified linear unit (ReLU) as the activation function because of several attractive advantages (Hara et al., 2015): (i) They possess stronger ability to represent the learned features; (ii) No gradient vanishment phenomenon and the convergence speed of the model is maintained at a stable state; and (iii) The neurons in the neural network have a sparse activation property and the sparse model enables better mining of relevant features and fitting of training data. The form of the function is shown as Equation (5).

To alleviate the influence of class imbalance, we randomly selected the same number of normal and pathological fetuses, which were both 105. The 1-dimensional preprocessed FHR signal was then transformed into a 2-dimensional RP image. By changing three optional parameters (m = 2, 3; τ = 1 to 10 in step of 1; k = 1 to 10 in step of 1), the ultimate dataset contained 21, 000 images for each class; i.e., 21,000 normal and 21,000 abnormal cases, which were considered to be sufficient for DL training.

In this work, the training strategy of 10-fold cross-validation was employed in the performance evaluation to obtain more reliable results. The total RP images were separated equally into 10 segments, thereby 9 of 10 images were used in the training and validation of the CNN model while the remainder (1 of 10) of the images were used to test the performance of the proposed system. This process was iterated 10 times by randomly shifting the training data, as shown in Figure 5. The final results recorded in all 10 iterations were averaged and considered to be the overall classification performance.

As with the conventional ML methods, we need some indicators to measure the performance of the proposed system. The confusion matrix is commonly used in a binary classification problem and it consists of four parameters: True Positive (TP), False Positive (FP), False Negative (FN), and True Negative (TN). The literal meaning of above mentioned positive and negative terms denote normal and hypoxia fetuses, respectively. In clinical practice, the accuracy (Acc), sensitivity (Se) and specificity (Sp) are specifically employed to assess the classification algorithms, as shown in Equations (7) to (9). The Se points out the performance of the model on the detection of normal fetus whereas the Sp measures the performance of the system on detecting a hypoxic fetus. In addition, the area under (AUC) the receiver operating characteristic (ROC) curve is also useful for obstetricians.

In this study, we trained our CNN model on a workstation with an Intel Core 3.70 GHz (i3-4170) processor and a 4GB RAM. The entire process of the proposed approach was implemented with MATLAB^® (2018a).

It is generally acknowledged that the CNN algorithm has many tuning parameters that can influence the classification performance to different degrees. In this work, the initial learning rate was set to 1 × 10^-3, which controls the relatively stable learning speed. L2 regularization was applied to overcome overfitting with a factor of 1 × 10^-4. The training and validation process of the CNN model is presented in Figure 6. It can be clearly observed that, as the iteration (epoch) progresses, the Acc increases and the loss decreases for both training and validation. To obtain better results, an optimization experiment was employed using the validation set in this work divided into three primary aspects.

First, we experimented with the effect of different layer parameters on performance, which were adjusted in turn. For the convolution layer of the CNN model, the size and number of kernels had a certain relationship with the classification performance, as shown in Figure 7. Based on the validation set, we discovered that the Acc of the size of the convolution kernel of 5 × 5 and 7 × 7 remained approximately the same but was higher than that of 1 × 1 and 3 × 3. The validation Acc increased with the number of kernels before reaching 8. Thereby, we set the size and number of kernels to be 5 × 5 and 8 for higher Acc and less training time (as indicated in Figure 7), respectively. In addition, we found that the parameters of the stride and padding factor in the convolution layer, the stride, padding factor and size of pooling kernel in the pooling layer made the Acc basically stable. In summary, Table 2 presents the determined parameters for each specific layer corresponding to Figure 4.

Second, we experimented with the effect of size of the mini-batch and max epoch on the training process, as shown in Figure 8. It can be seen that regardless of the max epoch, the validation Acc was higher when the size of the mini-batch was equal to 64. On the other hand, we observed the trained CNN model might not fully learn the image information when the max epoch was 5. Once the max epoch reached 10, the model learned as many features as possible about the fetal state from the input data and the Acc remained essentially constant with increases of the max epoch. Thus, the size of the mini-batch and max epoch were selected as 64 and 10 (as indicated in Figure 8), respectively.

Third, the training of kinds of multilayer neural networks requires the back propagation (BP) algorithm to train the feed-forward neural network better and faster and to obtain the weight and bias parameters of each layer of the neural network. In recent years, researchers have proposed several powerful BP algorithms for the training of a CNN, including stochastic gradient descent (SGD), adaptive moment estimation (ADAM), and root mean square propagation (RMSP). Table 3 presents the three algorithms together with individual parameter settings and the corresponding results. It was obvious that the Acc of the ADAM algorithm was higher than the other two.

Based on the same optimal experimental method described in Experiment 1, we investigated the influence of different layers of the CNN model with their corresponding CNN parameters on the classification performance and Table 4 presents the experimental results using the testing set. Architectures with more than 12 layers demonstrated overfitting or underfitting and therefore were not considered. Clearly, in the beginning, as the number of layers increased, the performance of our proposed system improved. When the number of layers reached 8 layers, consisting of one input layer, two convolution layers, two pooling layers, two fully connected layers, and one output layer (i.e., 13 layers including the normalization and ReLU layers followed by each convolution layer and the dropout layer located in the successive fully connected layers), the best performance was achieved: of 2,100 normal fetuses, only approximately 0.95% (20) were wrongly classified as pathological. Likewise, for hypoxia fetuses, a total of 2.33% (49) were wrongly classified as normal. Moreover, as the number of layers continued to increase, the system performance began to decline instead. Figure 9 shows the ROC curve with the optimal AUC value (CI: 0.9797–0.9875) when the number of layers was 8. In addition, it can be observed from Table 4 that different layers of CNN architecture had approximately the same testing time for one-fold. In other words, once the proposed CNN algorithm was successfully trained, the corresponding CAD system could immediately identify an unknown fetus regardless of the number of layers.

Furthermore, we conducted the experiment concerning the impact of the signal length on the classification performance. According to the optimal CNN configuration, Figure 10 shows that when the signal length is less than 8 min, the Acc increases gradually; and when the signal length is greater than 8 min, the Acc fluctuates within a small range. Therefore, we discovered that 8 min was the minimum length of the FHR signal containing relevant information regarding the fetus, and 13 min achieved better performance (Acc = 98.69%) using the CNN for fetal state classification.