Noisy Speech Based Temporal Decomposition to Improve Fundamental Frequency Estimation

The first step is devoted to the Ensemble Empirical Mode Decomposition (EEMD). The general idea of the EMD [25] is to analyze a signal $x(t)$ between two consecutive extrema (minima or maxima), and define a local high-frequency part, also called detail $d(t)$, and a local trend $a(t)$, such that $x(t)=d(t)+a(t)$. An oscillatory IMF is derived from the detail function $d(t)$. The high- versus low-frequency separation procedure is iteratively repeated over the residual $a(t)$, leading to a new detail and a new residual. Thus, the decomposition leads to a series of IMFs and a residual, such that

where $\mathrm{IMF}_{k}\left(t\right)$ is the $k$-th mode of $x(t)$ and $r(t)$ is the residual. In contrast to other signal decomposition methods, a set of basis functions is not demanded for the EMD. Besides, this strategy leads to fully data-driven decomposition modes and does not require the stationarity of the target signal.

The EEMD was introduced in [36] to overcome the mode mixing problem that generally occurs in the original EMD [25]. The key idea is to average IMFs obtained after corrupting the original signal using several realizations of White Gaussian Noise (WGN). Thus, EEMD algorithm can be described as:

1. 1.

   Generate $x^{n}\left(t\right)=x\left(t\right)+w^{n}\left(t\right)$ where $w^{n}\left(t\right)$, $n=1,\cdots,N$, are different realizations of WGN;

2. 2.

   Apply EMD to decompose $x^{n}\left(t\right)$, $n=1,\cdots,N$, into a series of components $\mathrm{IMF}_{k}^{n}\left(t\right)$, $k=1,\cdots,K$;

3. 3.

   Assign the $k$-th mode of $x(t)$ as

   $$\mathrm{IMF}_{k}\left(t\right)=\frac{1}{N}\sum_{n=1}^{N}\mathrm{IMF}_{k}^{n}\left(t\right);$$

   (2)

4. 4.

   Finally, $x\left(t\right)=\sum_{k=1}^{K}\mathrm{IMF}_{k}\left(t\right)+r\left(t\right)$, where $r(t)$ is the residual.

Although the lack of mathematical formulation, it is interesting to point out that EEMD is a very powerful tool for analyzing non-stationary real signals and has been successfully applied in several research areas [28], [37].

In this step, the F0 of frames of each IMF is estimated by the selected PEFAC [16] algorithm. This solution detects the F0 by convolving the power spectrum of the frame in the log-frequency domain, using a filter that sums the energy of the pitch harmonics. This filter highlights the harmonic components, improving the noise-robustness of the algorithm. Thus, let $\hat{F}0_{k,q}$ denote the F0 value estimated from frame $q$ of $\mathrm{IMF}_{k}\left(t\right)$, the $\hat{F}0_{q}$ vector is composed as

to express the tendency that the frame is placed in a low or high frequency region. Although the first IMFs are composed of fastest oscillations, these components also present some nonnegligible low-frequency content of the speech signal [24], [38]. Therefore, in this paper it is considered only the first four IMFs ($K=4$) in order to avoid the acoustic noise masking effect, whose energy is mostly concentrated at low frequencies [26], [38], [39].

Fig. 2 shows the F0 estimated with PEFAC for the first four IMFs of a voiced speech segment collected from the CSTR database [32] and corrupted by Babble [40] noise with SNR = 0 dB. Default parameters of PEFAC are considered, whose voiced segment is split into overlapping frames with 90 ms length and 10 ms frame shift. Note that for the analyzed IMFs, the F0 contours are close to the ground truth. In most frames the frequencies have similar values for all the IMFs, particularly when the decomposition attenuates the noisy component of speech. However, some noisy masking effect should appears on IMFs, leading to a difference in the F0, as seen in Fig. 2(a),(d) (IMF₁ and IMF₄) in the frames between 150-200 ms.

A normalized distance is computed between IMFs for the successive frames, in order to detect and overcome the differences in the estimated F0. Let $k$ and $k^{\prime}$ denote IMF indexes, the distance $\delta^{q}_{\hat{F}0}$ is described as

This distance is determined for different values of $k$ and $k^{\prime}$, resulting in the following matrix

The row components of the matrix are summed, and the resulting values obtained express the variation property for the $k$-th IMF. The two IMFs with the smallest variation scores are selected, and the frequency region is defined as the mean value of PEFAC F0 estimates ($\bar{F}0_{q}$). Finally, a low-frequency to high-frequency threshold $\gamma$ is adopted, and the separation is performed such as

In [41] it is shown that the variability of speech F0 is between 50-200 Hz for men and 120-350 Hz for women. Therefore, in this study, a threshold of $\gamma$ = 200 Hz is considered, since this is an average value for both genders of speakers.

The F0 estimation techniques are prone to different types of errors. Typical errors in F0 estimation, e.g., halving and doubling errors, are caused by the detection of harmonics other than the first, resulting in F0 estimates that are multiples of the true F0 [29]. In order to overcome this issue, the proposed method aims to correct the frequency candidates extracted by a F0 detection algorithm according to these types of errors. The PRO method states that, the F0 candidates ($f_{cand}$) must lie in [50,200] Hz or [200,400] Hz, for low-frequency and high-frequency frames, respectively. Thus, considering a low-frequency frame, $f_{cand}$ is corrected following the criteria

where $f_{c}$ is the corrected F0 candidate. Finally, the high-frequency frame is adjusted as follows:

Fig. 3 illustrates the F0 candidates correction procedure for a speech segment corrupted by Babble noise with SNR = 0 dB, with 200 ms duration. Fig. 3(a) denotes the F0 candidates extracted using the HHT-Amp [24] estimation technique, which computes a set of three candidates (one candidate from each of the first three IMFs) for a 10 ms time interval. Ground truth indicates that the entire speech segment is classified as high-frequency. Note that most candidates are located in the low-frequency region ($f_{cand}\leq$ 200 Hz). Fig. 3(b) shows the new candidates, which are corrected according to the criteria described in (7) and (8). Moreover, observe that the adjusted F0 candidates take place in high-frequency region, matching the ground truth.

This technique [30] consists on training a DCNN, in order to classify a voiced speech frame into low ($F0\leq 200$ Hz) or high ($F0>200$ Hz) frequency. The DCNN architecture adopted is based on VGGNet [42], with six convolutional layers, three Fully-Connected (FC) layers and an output classification layer (Softmax). The DCNN input sequence is a 60 ms extracted directly from the voiced speech samples with a sampling rate of 16 kHz, or 960 samples.

According the low- or high-frequency classification, new F0 candidates $f_{j}$ are extracted based on initial estimation attained by a conventional method $f_{p}$. The improved F0 estimate value is selected from the set of new candidates by exploiting a restrained selection procedure. Two spectral attributes are introduced to assist in selecting the improved F0 among candidates: Weighted Euclidean Deviation ($d_{j,q}$) and Weighted Comb Filtering ($y_{j,q}$). The first attribute is computed observing the first five frequency peaks positions in each frame, resulting in a peak vector ${P}_{i,q}=[P_{1,q},P_{2,q},\dots,P_{5,q}]$. For a pitch candidate $f_{j,q}$, a candidate vector is defined as $V_{j,q}=[v_{j1,q},v_{j2,q},\dots,v_{j5,q}]$, where $v_{ji,q}$ is the point-wise ratio between the peak vector and the candidate $f_{j,q}$. Thus, the Weighted Euclidean Deviation is described as

where $T=[1,2,3,4,5]$ denotes the multiple harmonics of pitch, $U_{i}=1/T_{i}$ and $\odot$ is the point-wise multiplication. The Weighted Comb Filtering for $f_{j,q}$ is defined as

where $X_{q}\left(f\right)$ is the power spectrum for the $q$th frame, and $C\left(f/f_{j,q}\right)=\left(\frac{1}{2}+\frac{1}{2}\mbox{cos}\left(2\pi f/f_{j,q}\right)\mbox{exp}\left(-f/f_{j,q}\right)\right)$. Finally, a cost function is defined as

where $\left|\mbox{log}f_{j,q}-\mbox{log}f_{i,q+1}\right|$ is a F0 smoothness feature, $\lambda>0$ and $\alpha$ are regularization parameters, $pr_{q}$ is the low/high frequency probability of the softmax layer of DCNN, and $\varepsilon_{q}=1$ if $f_{j}=f_{p}$, or zero otherwise. The smaller cost function value for the F0 candidates is more likely to be the true F0 value.

SHR method [14] is based on the definition of a parameter measure called Sub-Harmonic-Harmonic Ratio, designed to describe the amplitude ratio between subharmonics and harmonics. Let $A(f)$ denote the amplitude spectrum for each short-term signal, the sum of harmonic amplitude is described as

where $N$ and $F0$ are the maximum number of harmonics considered in the spectrum and fundamental frequency, respectively. In other hand, the sum of subharmonic amplitude is achieved assuming at one half of F0

Finally, $SHR$ is the ratio between $SS$ and $SH$:

In [14], this ratio is computed on a logarithmic frequency scale. If the obtained $SHR$ value is greater than a certain threshold, subharmonics frequencies are considered in the analysis. Otherwise, the harmonic frequencies are adopted in the F0 detection.

The core idea of SWIPE [15] is that if a signal is periodic with fundamental frequency $f$, its spectrum must contain peaks at multiples of $f$ and valleys in between. Since each peak is surrounded by two valleys, the Average Peak-to-Valley Distance (APVD) for the $n$-th peak is defined as

where $X$ is the estimated spectrum of the signal which takes frequency as input and outputs corresponding density. The global APVD is achieved by averaging over the first $p$ peaks

The F0 estimated is the $f$ value that maximizes this function, searching $f$ in the range [50 500]Hz candidates, with samples distributed every $1/48$ units on a base-2 logarithmic scale.

The HHT-Amp method is summarized as follows:

• •

  Apply the EEMD as described in Section II-A to decompose the voiced sample sequence $x_{q}(t)$.

• •

  Compute the instantaneous amplitude functions $a_{k,q}(t)=|Z_{k,q}(t)|,k=1,\ldots,K,$ from the analytic signals defined as

  $$Z_{k,q}(t)=\mbox{IMF}_{k,q}(t)+j\,H\{\mbox{IMF}_{k,q}(t)\},$$

  (17)

  where $H\{\mbox{IMF}_{k,q}(t)\}$ refers to the Hilbert transform of $\mbox{IMF}_{k,q}(t)$.

• •

  Calculate the ACF $r_{k,q}(\tau)=\sum_{t}a_{k}(t)\,a_{k}(t+\tau)$ of the amplitude functions $a_{k,q}(t),k=1,\ldots,K$.

• •

  For each decomposition mode $k$, let $\tau_{0}$ be the lowest $\tau$ value that correspond to an ACF peak, subject to $\tau_{min}\leq\tau_{0}\leq\tau_{max}$. The restriction is applied according to the range $[F_{min},F_{max}]$ of possible $F_{0}$ values. The $k$-th pitch candidate is defined as $\tau_{0}/{f_{s}}$, where $f_{s}$ refers to the sampling rate.

• •

  Apply the decision criterion defined in [24] to select the best pitch candidate $\hat{T}_{0}$. The estimated $F_{0}$ is given by $\hat{F}_{0}=1/\hat{T}_{0}$.

In [24], it was shown that the HHT-Amp method achieves interesting results in estimating the fundamental frequency of noisy speech signals. HHT-Amp was evaluated in a wide range of noisy scenarios, including five acoustic noises, outperforming four competing estimators in terms of GE and MAE.

The experiments consider the CSTR [32] and a subset of TIMIT [33] databases to evaluate the competitive methods. CSTR is composed of 100 English utterances spoken by male (50) and female (50) speakers, sampled at 20 KHz. The reference F0 values are available based on the recordings of laryngograph data. The TIMIT subset is composed of 128 speech signals spoken by 8 male and 8 female speakers, sampled at 16 KHz and with 3 s average duration. The reference F0 values are obtained from [43].

Six noises are used to corrupt the speech utterances: acoustic Babble and Volvo attained from RSG-10 [40], Cafeteria, Train and Helicopter from Freesound.org¹¹1[Online]. Available: https://freesound.org., and Speech Shaped Noise (SSN) from DEMAND [44] database. Experiments are conducted considering noisy speech signals with five SNR values.

For the evaluation, it is adopted the Gross Error rate (GE) and Mean Absolute Error (MAE). GE [34] is defined as:

where $P$ denotes the total number of voiced frames, and $P_{error}$ is the number of voiced frames for which the deviation estimated F0 from the ground truth is more than 20%. MAE [35] is computed as

where $Q$ is the total number of frames, $\hat{F}0\left(q\right)$ the estimate and $F0\left(q\right)$ is the reference. This metric provides a greater perception of the error, since it indicates an absolute distance (in Hz) between F0 reference and estimation.

Table I presents the error results with the PRO and the competing DCNN technique for low/high frequency separation of voiced speech frames. These results denotes the mean values for the utterances of CSTR and TIMIT databases, respectively, considering the speech signals corrupted by four noises and five SNR values. Note that PRO outperforms the DCNN for all of the noisy scenarios. The proposed method achieves the lowest error values in all the 20 noisy conditions and the two databases. For instance, it can be seen that PRO attains interesting results for the most severely conditions, for SNR = -15 dB. In this case, the error scores are about 28 percentage points (p.p.) smaller than those achieved by DCNN approach in CSTR database. The overall errors obtained with PRO are 14.7% and 11.6% for CSTR and TIMIT databases, against 40.2% and 17.0% for DCNN, respectively.

The interesting accuracy results attained by PRO method are particularly important, and can be justified by the fact that it is considered only the first four IMFs of EEMD. This attenuates the noise noise masking effects, since the most part of noisy energy is concentrated in low frequencies, e.g., IMF${}_{\text{6}}$, IMF${}_{\text{7}}$ and IMF${}_{\text{8}}$. Moreover, PEFAC also contributes to the separation accuracy, due to the filter introduced in the algorithm, which rejects high-level narrow-band noise and favors the correct F0 estimates, even in severe noisy environments.

The proposed method and DCNN baseline technique for improvement of F0 detection accuracy are compared in terms of GE and MAE, considering three F0 estimation approaches: SHR, SWIPE and HHT-Amp. In this study, it is assumed that errors in separation of the voiced speech into low-frequency and high-frequency signals can generate some error into the whole system. Furthermore, voiced/unvoiced detection has been done by the strategy based on the Zero Crossing (ZC) rate and energy of the speech signal as described in [45].

Fig. 4 depicts the GE values for the F0 estimation approaches improved by PRO method, in comparison with the performance of DCNN based competitive techniques, averaged over signals of CSTR database. It is interesting to mention that PRO + HHT-Amp achieves the lowest error results, being the most accurate method when compared to baseline solutions. In addition, note that PRO method attains interesting F0 improvement for SHR and HHT-Amp F0 estimators, even assuming the low/high frequency separation errors. For instance, the original HHT-Amp approach attained a GE of 59.6% for Helicopter noise (Fig. 4(e)) with SNR = -15 dB, in contrast with 39.7% of PRO + HHT-Amp scheme. For the SHR method, DCNN improved the F0 estimates in some cases, e.g., Train and Helicopter noises, but it is still surpassed by PRO. However, for SWIPE and HHT-Amp, the errors of DCNN separation affect whole system, causing an increasing in the GE scores. Note that, in severe noisy condition, the SWIPE method obtained the highest GE values. This fact may be explained by its inner voiced/unvoiced detection that is impaired by the background situation.

Table II presents GE results of PRO and DCNN fundamental frequency improvement methods obtained with speech signals of TIMIT database. Note that PRO outperforms DCNN for SHR and HHT-Amp detectors. For example, considering the Helicopter noise with SNR = -15 dB, the GE rate for the SHR F0 estimation scheme decreased from 86.0% with the DCNN [30] to 71.4% with PRO, i.e., a reduction of 14.6 p.p. Likewise the results exposed for CSTR, PRO + HHT-Amp achieves the best accuracy results for TIMIT database, with lowest GE values in all the 30 noisy conditions. For the Babble noise, this strategy attains the average GE of 31.3% against 46.9% for PRO + SHR. On overall average, PRO + HHT-Amp presents a GE score of 23.2%, which is 18.0 p.p., 18.2 p.p and 32.6 p.p smaller than DCNN + HHT-Amp, PRO + SHR and DCNN + SHR baseline approaches, respectively.

Fig. 5 and Fig. 6 illustrate the MAE results obtained for six noisy conditions with the CSTR and TIMIT databases, respectively. Each box-plot denotes the MAE values achieved with five SNR values: -15 dB, -10 dB, -5 dB, 0 dB and 5 dB. Note that again the proposed method achieved the lowest MAE values. For CSTR database, the DCNN strategy attains improvement in original F0 estimation in some cases of SHR and SWIPE, particularly for SSN, Cafeteria, and Helicopter noises. The MAE values for the TIMIT database indicate that PRO improved the F0 estimation accuracy, overcoming the DCNN baseline solution. Once again, the PRO method leads the HHT-Amp estimator to achieve best accuracy results for both databases, with exception for Volvo noise with TIMIT (Fig. 6(f)).

In summary, the proposed method leads to the best results when compared to the solution DCNN, in terms of low/high frequency separation accuracy, and improvement of fundamental frequency estimation accuracy. In addition, the combination of PRO with the HHT-Amp algorithm (PRO + HHT-Amp) attained the lowest GE and MAE scores, for both the CSTR and TIMIT databases. Furthermore, SHR baseline approach are also outperformed by the HHT-Amp. The fact that SHR adopts two F0 candidates, against the three candidates considered in HHT-Amp, can favor this last one in the true F0 selection. In contrast, baseline techniques based on the SWIPE estimation reach the highest error values. This occurs because SWIPE is very sensitive to the severe noise disturbances, which may cause expressive voiced/unvoiced detection errors. Finally, it is interesting to observe that even considering the classification errors in the estimation accuracy results, the proposed method improves F0 estimation of speech signals corrupted by different types of noise, and hence it can be used in real world applications.

Table III indicates the computational complexity which refers to the normalized processing time required for each scheme evaluated for 512 samples per frame. These values are obtained with an Intel (R) Core (TM) i7-9700 CPU, 8 GB RAM, and are normalized by the execution time of the most accurate method PRO + HHT-Amp. The processing time required for the training process by DCNN is not considered here. Therefore, note that the HHT-Amp baseline scheme and the proposed method present a longer processing time, since they are based on the EEMD, and demands a relevant computational cost.