Performance Analysis of Fetal-Phonocardiogram Signal Denoising Using The Discrete Wavelet Transform

— The obligation for comprehensive fetal heart rate investigation had driven to improve the passive and non-invasive diagnostic instruments despite the USG or CTG method. Fetal phonocardiography (f-PCG) utilizing the auscultation method met the above criteria, but its interpretation frequently disturbed by the presence of noise. For instance, maternal heart and body organ sounds, fetal movements noise, and ambient noise from the environment where it is recording are the noise that corrupted the f-PCG signal. In this work, the use of discrete wavelet transforms (DWT) to eliminate noise in the f-PCG signal with SNR as the performance parameters observed. It was observing the effect of changes in wavelet type and threshold type on the SNR value. The test was carried out on f-PCG data taken from physio.net. Initial SNR values ranged from -26.7 dB to -4.4 dB; after application of DWT procedure to f-PCG, SNR increased significantly. Based on the test results obtained, wavelet type coif1 with the soft threshold gave the best result with 11.69 dB in SNR value. The coif1 had a superior result than other mother wavelets that use in this work, so the fPCG signal analysis for fetal heart rate investigation suggested to use it.


INTRODUCTION
Long-term observation of fetal heart rate (FHR) in the womb can be used to provide the health conditions of the fetus and mother required for pregnancy diagnosis. There are several generally known devices to diagnose fetal conditions, e.g., ultrasonography (USG), cardiotocography (CTG), and fetalelectrocardiography (f-ECG) [1], [2], [3]. The disadvantage of the use of USG or CTG for long-term observation is that the instrument is expensive, and well-skilled expertise is needed for obtaining and evaluating data. So, the more effective and more efficient device is needed to resolve such auscultation method [4]. The auscultation method uses an acoustic signal, a passive and non-invasive method, in the abdomen of maternal women to observe the pregnancy. Measurement of the acoustic signals in the abdominal area of maternal women can provide the FHR condition through recorded passive energy (known as the fetal-phonocardiography (fPCG) method) [2]. Typically, f-PCG has two fetal heart sounds, the first heart sound (S1) is produced by the closure of atrioventricular valves, and the second (S2) represents the closure of semilunar valves [4]. Considering its feature, f-PCG is a potential clinically useful test to evaluate fetal heart rate (FHR) during pregnancy and to assess fetal wellbeing during pregnancy, labor, and delivery [4]. However, f-PCG have an issue about the existence of other signals that coincide with the fetus heart sound such as maternal heart sounds (mHS), maternal organs sound (i.e., maternal digestion, respiratory movements), fetal movements, sensor movements, and ambient (environment) noise [5]. As it is heavily contaminated by noise, measurement fPCG using auscultation method implies mandatory to eliminate those signals from the fetus heart sound signal. Thus, appropriate filtering has to be applied in order to make fPCG clinically usable and to extract important diagnostic information, such as FHR [4]. The development of a method to eliminate noise on fPCG signals had been carried out by several researchers before. Varady et al. [2] showed the use of the wavelet method on two PCG channels and signal processing using the adaptive wavelet approach to obtain FHR. Kovacs et al. [6] developed the method that combines autocorrelation techniques, wavelet transforms, and matching pursuit to evaluate FHR. The use of other wavelet transform methods was carried out by Jaros et al. that observe wavelet performance compared to FIR [7]. In this work, the observation results of a discrete wavelet-based transformer noise (DWT) system by varying several wavelet-transform parameters such as wavelet type and threshold selection. Performance measurement was performed on the SNR value of the decomposed signal. Heriana and Misbah [8] stated that DWT has more advantages in the process of eliminating noise because it can analyze signals in the time domain and frequency domain simultaneously and can analyze in several structures. The present work explores wavelet transforms to reduce noise in fPCG signals.
Performance observations were made of thirty-seven results of recording fPCG signals taken from physio.net secondary data providers [9][10].

II. RESEARCH METHOD
Denoising is the process of basic signal extracting from the mixed signals. Regardless of the frequency information content, all parts of the signal that are considered noise will be eliminated, so a new signal is obtained with the required features (such as S1 and S2 nodes in fPCG). The fPCG signal used in this study was taken from the Physionet.org website [11], which provides datasets that are quite complete and freely accessible. This data is called Simulated Fetal Phonocardiograms (simfpcgdb), which consists of 37 fPCG signals with details in Table 1 [12] that have been mixed with noise with SNR variations from -26.7 dB to -4.4 dB. The fPCG signal is a signal of the measurement of fetal heart rate (non-invasive) through a maternal woman. It is a fetal heartbeat signal that is being conceived so that the mother's organs such as the skin mediate between the sensor and the fetus. There are several types of noise-producing activities for fPCG signals, namely [2] [14] maternal heart sounds (mHS), maternal organs sound (i.e., maternal digestion, respiratory movements), fetal movements, sensor movements, and ambient (environment) noise.
Thus, the denoising procedures to obtain fPCG signal are required before a data from fPCG measurement results is analyzed, so that fPCG nodes such as S1 and S2 can be detected with high accuracy. In this research, a denoising process has been carried out using Discrete Wavelet Transform with variations of the mother wavelet to obtain the appropriate type of mother wavelet for denoising fPCG signals. The denoising system was developed using Matlab 2018a. Thus, the amplitude and sampling frequency information can be read directly by Matlab. In this work, the procedure of denoising fPCG signals is shown in Fig. 1. Meanwhile, the stages of the Signal Denoising process are shown by the flowchart in Fig.  2. To extract the fPCG signal, this work had several steps started by preprocessing using normalized and bandpass filter procedure. To make a threshold, the noise level estimation procedure had carried out and the next step was generating the global threshold used to denoising the process. After the denoising process, its performance was described by SNR values.
Signal denoising procedure is described in Fig. 2 that consists of three steps: forward wavelet transform, wavelet denoising, and inverse wavelet transform. The initial parameter for this procedure was obtained from the bandpass filtering process and the initial parameter that declares at the beginning, such as wavelet type, wavelet level, and threshold.
All fPCG signals were processed by employing a discrete wavelet transform (DWT). The fPCG signal of measurement results containing noise is defined by ( ), the signal without noise is defined by ( ), and the noise is symbolized by ( ) . The relationship between the three variables is expressed by (1).  where σ is the noise level, the ξ(n) signal is a normally distributed random signal with an average of 0, and a standard deviation of 1 or is usually written as N (0, 1). The s(n) signal is the signal used in this study and then added with the noise ξ(n) to get the x(n) signal. By denoising, it is expected that the signal s(n) can be retrieved, which is symbolized by s^' (n) with the quality of the results measured based on the SNR parameter.

A. Preprocessing
The preprocessing phase was used to make the signal in the conditioning stage. The signal will sequentially process in two-step: a) Normalize Signal Firstly, the fPCG signal was normalized to get a signal in the range of -1 v to + 1 v. This process minimized the effect of anxiety signal amplitude. The fPCG signal was measured using sensors influenced by the activity of maternal and fetal organs. All of those activities were mixed in one recording/measurement that causes greater signal variability. For this reason, normalization was needed to reduce this variability. The normalization is defined as (2).
where ̃[ ] is a normalized signal. Although the simfpcgdb signal has been normalized, this step was a procedure always used as the beginning of processing.

b) Bandpass Filter
The noise in fPCG signals is usually below 20 Hz [13] (originating from the mother's internal organs) and above 70 Hz [2], [14], [15], [16] (originating from external factors). In research [17], the noise of your internal organs occupied a frequency of less than 25 Hz. Therefore, only certain frequency ranges will be taken using a bandpass filter. In this work, the Butterworth IIR bandpass filter was used which the filter order (M) was 10, and the frequency range corresponding to the simfpcgdb data was 30-80 Hz. This filter can also eliminate the influence of low frequency caused by signal bias (offset), so the expected value of the signal was 0. The frequency response of this filter is shown in Fig. 3 (a) and (b).
In Fig. 3 (a), the filter used to maintain the gain in the passband region (0 dB in the dotted line region), which means there is no change in the signal magnitude in the frequency range. While the signal outside the desired frequency range (stopband) has weakened to -50 dB. Through Fig. 3 (b), it appears that the phase shift is almost flat (flat) shown by the dotted line, which means that each frequency in the passband does not experience an angle shift.
Mathematically, the Butterworth IIR filter differential equation is formulated by (3).
where y[n] is the filtering signal. Variables a_i and b_i, where i,j=0,1,…,M is the feedforward and feedback filter coefficients.

B. Noise Level Estimation
Thresholding in the wavelet transform was divided into two, namely global thresholding and leveldependent thresholding method. The thresholding process and its value make a difference between global and level-dependent thresholding. In the first method, the threshold value is calculated based on the entire data, including the threshold value. The second method, the threshold value, is different at each level (level), including the threshold carried out for each level. In this study, the method used was global thresholding.
Both methods involve the noise effect described in the noise level for the thresholding process. The noise level can be calculated with a statistical parameter such as standard deviation (σ), but it is not recommended because it requires the flat initial signal [ ] . In [18], another approach was shown for calculating the level of noise, as expressed in (4).
where σ ̃ is the estimated noise level, Q_nis the coefficient of detail at the first level (finest level) wavelet decomposition. The aim of using the first level wavelet was that the noise composition at this level was more dominant than the level afterward. Assuming the most dominant noise was in the highfrequency spectrum, this coefficient was used in the calculation. Consistently with the fact, this coefficient was the result of the wavelet bandpass filter. To calculate the noise level in this study, we used the db1 mother wavelet type. If the data length of y[n] was N then the number of coefficients involved in this calculation was N/2 because the decomposition of the first level wavelet divides the signal into half of the original signal length.

C. Calculate Global Threshold (Calculate Global
Threshold) Determining the threshold values are used by the Universal Rule principle [19] expressed as (5).

=̃√2
( ) is the threshold value calculated based on the level of noise that had been found through (4)? When the signal was normalized to the noise deviation standard, the variability of noise level was nonessential value. This threshold value was used for all wavelet detail coefficients.

D. Signal Denoising
The denoising process using wavelets consists of several stages: a) Feedforward Wavelet Transform The noisily fPCG signal [ ] was processed to wavelet transformation with decomposition level (L) was four, and mother wavelets were Daubechies, Symlets, Coiflets family. The aim of choosing the type of mother wavelet was the investigation of the most suitable for denoising in this study. These three wavelet families were chosen because they are most suitable for denoising adult [20] and fetus signals [21], [22]. Wavelet as a filter bank consists of a low pass filter (LPF) and a high pass filter (HPF). The output of LPF was the approximation coefficient (cA), while the output of HPF was the coefficient of detail (cD). The relationship between the two groups of coefficients was illustrated by (6). Noise usually occupies the high-frequency spectrum, so the denoising process was carried out in this region. Using the 3dB-frequency was 30-80 Hz, whereas the sampling frequency was 1 kHz. By following the Nyquist theory, the maximum frequency of fPCG recordings was 500 Hz. Respectively value of these frequencies, the value above 80 Hz, can be said as the high-frequency spectrum.

b) Denoising
In this study, the denoising process was applied for all the detail coefficients (cD) of wavelet decomposition products. Besides the approximation coefficient (cA) was processed only for the last level coefficient J. The threshold value used was and applied to all levels of decomposition, which was also called the independent thresholding level [19]. The threshold technique used consists of 2 namely hard thresholding and soft thresholding. The hard thresholding technique based on one threshold value shown by (7) [19].

E. Signal to Noise Ratio
The denoising system performance was measured through Signal to Noise Ratio (SNR). In [11] did not provide the original signal contaminated by noise, so it was difficult to compare the denoising signal to the reference signal. Because of the absence of the reference signal, SNR is calculated using the mathematical expression used by [23] in (10).

III. RESULT
The fPCG signal used in this study was a simulated fPCG signal (simfpcgdb) with a total of 37 signals mixed noise with SNR variations. These signals were implemented in the denoising process following the stages that have been designed. In Fig. 5, an fPCG signal with the file name fetalPCG_simulatedSNR-8_1dB recorded at a sampling frequency (fs) = 1 kHz and duration of 30 seconds. In this figure, we just showed the data for the 1-second duration to facilitate observation in the time zone. The simfpcgdb data was downloaded online and automatically using the wfdb toolbox [24]. The downloaded signal looks like a coherence normalized rectangular signal. Based on the Fourier theory, a summary of infinite number signal with different frequencies will form a square signal when it was added, so it is very reasonable when the signal was charged with a bandpass filter at the initial stage. A Butterworth IIR filter has been designed to filter the signal so that it produces the signal in Fig. 6.
The output signal from the IIR filter showed a different signal from the original one. Even, Fig. 6 displayed some vertices S1 and S2, but there are still high-frequency components that need to be attenuate or remove. After the noise removing, the transitions between nodes (S1 to S2 or S2 to S1) became clearer. For this reason, wavelet transformation was used at the denoising stage.  Fig. 7 shows the comparison of fPCG signals after reconstruction using the two threshold techniques and the comparison with the results of the Butterworth IIR bandpass filter. In the time domain, the results of soft thresholding were smoother than hard thresholding in some parts of the signal. It happened because the nonzero value of the detail coefficient had changed mathematically in soft thresholding while it not happened in hard thresholding. The area that highlights with a box in Fig. 7 shows the signal part that was still clear on the hard thresholding results, but it attenuates in the soft thresholding result. However, soft thresholding was always better than hard thresholding because the aims of denoising fPCG signals and other biomedical signals were not only to eliminate noise but also to produce new representations of signals such as the clarity of the position of S1 and S2 signals. The use of soft thresholding tended to weaken the signal power, so S2 became vulnerable to losses due to shrinkage, which did not occur in the hard thresholding, as shown by the two-way arrows in Fig. 9 (a) and (b).
For the next stage, the SNR parameters were calculated for each decomposed signal reconstructed by the wavelet transform with different mother wavelet types and the thresholding method. The results of the SNR parameters are shown in Fig. 8. Based on Fig. 8, it can be concluded that the type of wavelet that had the largest SNR was coif1 from the coiflet family. The best SNR value for this wavelet type was 8,5365 dB for the hard thresholding and 11.6940 dB for the soft thresholding. Even though in Fig. 6 (a), the coif2 wavelet type had a better SNR value about four times, which is circled in red than coif1, but principally the coif1 had better SNR value than coif2. Hence, the coif1 mother wavelet was used as a wavelet type to reconstruct the denoising signal in the next discussion.
The coefficient of detail at all levels before and after the threshold was shown in Fig. 9. The two-way arrow in Fig.9 showed a sample of significant differences between the results of hard thresholding and soft thresholding. It was caused by the amplitude of the reconstruction signal from soft thresholding smaller than the hard thresholding. Hence, a combination of the highest levels approximation coefficient gave the smoother signal (not jaggy) in some samples. The filtering output can also be observed through the frequency spectrum. In Fig. 10 (a) -(d), the frequency spectrum of the original fPCG signal, and after the denoising process were displayed. The original fPCG signal spectrum contained a very dense and distributed frequency spectrum in the range of 0 -500 Hz, shown in Fig. 10 (a). The filtering process also proved the previous explanation related to the basic principles of Fourier analysis. By using the Butterworth IIR bandpass filter, only the frequency spectrum in the range 30 -80 Hz was passed. There was almost no frequency leak (spectrum leakage), as in Fig. 10 (b). The interesting results were the distribution of the frequency spectrum after denoising by hard thresholding and soft thresholding. As shown in the frequency spectrum in Fig. 10 (c) and (d), the frequency spectrum appears above 80 Hz, which was almost flat. The magnitude is indeed not too large compared to the fundamental frequency and harmony, but this confirms that the thresholding process carried out on all the coefficients of detail greatly affects the high-frequency spectrum. This is very reasonable (reasonable) because the coefficient of detail is the result of filtering the wavelet high-frequency escaping filter. When the coefficient values change -which was already in a certain frequency range -resulting in the creation of another high-frequency spectrum, which is also called a frequency leakage (frequency leakage). The spectrum can be removed by setting a certain threshold value (hard thresholding) as part of postprocessing but requires more in-depth and thorough research. Also, through this frequency spectrum, it can be concluded that the hard and soft thresholding techniques result in a weakening of the IIR filter spectrum. Frequency leakage resulting from soft thresholding is greater than hard thresholding because almost all of the detail coefficient values experience changes in value that are not experienced by hard thresholding.

V. CONCLUSION
The denoising process using wavelets was carried out in this study. A total of 10 types from three wavelet families had been tested to determine the proper type in terms of SNR values. Based on the results, mother wavelet coif1 from the coiflet family produced the best SNR value of 8,5365 dB for hard thresholding and 11.6940 dB for soft thresholding. This wavelet type dominated the largest SNR value compared to other wavelet types, including the other coiflet family. This wavelet type was recommended to use in fPCG signal analysis especially in the denoising process, to get the S1 and S2 nodes that are useful for calculating FHR (Fetal Heart Rate).