Bayesian Pitch Tracking Based on the Harmonic Model

TL;DR

This paper introduces a Bayesian fundamental frequency tracking method using the harmonic model and Markov processes, improving robustness and accuracy over existing methods, especially in noisy conditions.

Contribution

It proposes a novel Bayesian algorithm that incorporates temporal smoothness priors for fundamental frequency, model order, and voicing, enhancing estimation accuracy and noise robustness.

Findings

Reduces mean absolute errors by 15-36% compared to state-of-the-art methods.

Decreases gross errors by 20-26% in noisy conditions.

Achieves superior performance on speech and voice datasets, including Parkinson's disease voices.

Abstract

Fundamental frequency is one of the most important characteristics of speech and audio signals. Harmonic model-based fundamental frequency estimators offer a higher estimation accuracy and robustness against noise than the widely used autocorrelation-based methods. However, the traditional harmonic model-based estimators do not take the temporal smoothness of the fundamental frequency, the model order, and the voicing into account as they process each data segment independently. In this paper, a fully Bayesian fundamental frequency tracking algorithm based on the harmonic model and a first-order Markov process model is proposed. Smoothness priors are imposed on the fundamental frequencies, model orders, and voicing using first-order Markov process models. Using these Markov models, fundamental frequency estimation and voicing detection errors can be reduced. Using the harmonic model, the proposed fundamental frequency tracker has an improved robustness to noise. An analytical form of the likelihood function, which can be computed efficiently, is derived. Compared to the state-of-the-art neural network and non-parametric approaches, the proposed fundamental frequency tracking algorithm reduces the mean absolute errors and gross errors by 15\% and 20\% on the Keele pitch database and 36\% and 26\% on sustained /a/ sounds from a database of Parkinson's disease voices under 0 dB white Gaussian noise. A MATLAB version of the proposed algorithm is made freely available for reproduction of the results\footnote{An implementation of the proposed algorithm using MATLAB may be found in \url{https://tinyurl.com/yxn4a543}

Figures

Tables

Table 1. TABLE I: Total error ratio in colored noise

SNR		-5.00	0.00	5.00	10.00
PEFAC (90.5 ms)	Babble	0.42	0.38	0.34	0.29
	Factory	0.34	0.30	0.27	0.25
PEFAC (25.6 ms)	Babble	0.41	0.35	0.29	0.24
	Factory	0.30	0.25	0.22	0.21
SWIPE	Babble	0.50	0.42	0.29	0.19
	Factory	0.52	0.49	0.40	0.28
CREPE	Babble	0.40	0.29	0.21	0.16
	Factory	0.39	0.28	0.20	0.15
YIN	Babble	0.50	0.43	0.32	0.22
	Factory	0.50	0.43	0.32	0.22
Prew-Fast NLS	Babble	0.35	0.27	0.18	0.12
	Factory	0.28	0.20	0.14	0.11
Prew-Proposed	Babble	0.34	0.25	0.17	0.12
	Factory	0.28	0.20	0.15	0.12

Table 2. TABLE II: Gross error ratio in colored noise

SNR		-5.00	0.00	5.00	10.00
PEFAC (90.5 ms)	Babble	0.62	0.51	0.44	0.39
	Factory	0.56	0.47	0.41	0.38
PEFAC (25.6 ms)	Babble	0.72	0.65	0.60	0.57
	Factory	0.68	0.61	0.57	0.54
SWIPE	Babble	0.96	0.81	0.55	0.36
	Factory	1.00	0.94	0.76	0.54
CREPE	Babble	0.73	0.50	0.34	0.24
	Factory	0.75	0.53	0.36	0.26
YIN	Babble	0.95	0.83	0.61	0.42
	Factory	0.96	0.83	0.61	0.42
Prew-Fast NLS	Babble	0.57	0.41	0.30	0.24
	Factory	0.55	0.42	0.33	0.28
Prew-Proposed	Babble	0.53	0.36	0.27	0.24
	Factory	0.51	0.37	0.29	0.25

Table 3. TABLE III: Mean absolute value [Hz] in colored noise with oracle voicing detector

SNR		-5.00	0.00	5.00	10.00
PEFAC (90.5 ms)	Babble	49.81	39.15	31.73	27.96
	Factory	36.20	31.24	27.97	26.69
PEFAC (25.6 ms)	Babble	81.49	72.65	65.71	60.54
	Factory	72.61	64.93	57.93	54.20
SWIPE	Babble	31.73	17.94	10.95	8.04
	Factory	43.91	27.02	16.02	10.51
CREPE	Babble	68.95	44.93	30.57	21.89
	Factory	79.00	52.41	34.51	24.70
YIN	Babble	56.25	39.05	23.86	14.96
	Factory	57.37	38.53	23.41	14.97
Prew-Fast NLS	Babble	64.81	45.79	31.45	23.79
	Factory	74.58	57.88	44.93	36.50
Prew-Proposed	Babble	33.33	17.91	12.22	10.81
	Factory	19.32	13.20	11.23	10.48

Equations

$${\mathbf{y}}_n=h({\mathbf{x}}_n,{\mathbf{v}}_n),$$

$${\mathbf{x}}_n=f({\mathbf{x}}_{n-1},{\mathbf{m}}_n),$$

$${\mathbf{Y}}_{n-1}=[{\mathbf{y}}_1,\cdots ,{\mathbf{y}}_{n-1}].$$

$$p({\mathbf{x}}_n|{\mathbf{Y}}_{n-1})$$

$$=$$

$$p({\mathbf{x}}_1)=\int p({\mathbf{x}}_1|{\mathbf{x}}_0)p({\mathbf{x}}_0)d{\mathbf{x}}_0,n=1,$$

$$p({\mathbf{x}}_n|{\mathbf{Y}}_n)$$

$$p({\mathbf{x}}_1|{\mathbf{Y}}_1)=\frac{p({\mathbf{y}}_1|{\mathbf{x}}_1)p({\mathbf{x}}_1)}{p({\mathbf{y}}_1)},n=1,$$

$${\mathbf{y}}_n={\mathbf{s}}_n+{\mathbf{v}}_n,$$

$${\mathbf{y}}_n=[y_{n,1},\cdots ,y_{n,M}{]}^T,$$

$${\text{H}}_1:s_{n,m}$$

$${\text{H}}_0:{\mathbf{y}}_n$$

$${\text{H}}_1:{\mathbf{y}}_n$$

$$\mathbf{Z}({\omega }_0,K_n)$$

$$\mathbf{c}({\omega }_n)$$

$$\mathbf{d}({\omega }_n)$$

$${\mathbf{a}}_{K_n}$$

$${\mathbf{y}}_n$$

$$\lfloor F\frac{{\omega }_{\max }}{2\pi }\rfloor -\lceil F\frac{{\omega }_{\min }}{2\pi }\rceil +1,$$

$$p({\mathbf{x}}_n|{\mathbf{Y}}_{n-1})=$$

$$p({\sigma }_n^2|{\ddot{\mathbf{x}}}_n,{\mathbf{Y}}_{n-1})p({\ddot{\mathbf{x}}}_n|{\mathbf{Y}}_{n-1}).$$

$$p({\sigma }_n^2|{\sigma }_{n-1}^2)=\mathcal{I}\mathcal{G}({\sigma }_n^2|c,d{\sigma }_{n-1}^2).$$

$$p({\sigma }_n^2|{\sigma }_{n-1}^2,{\ddot{\mathbf{x}}}_n,{\mathbf{Y}}_{n-1})\propto 1/{\sigma }_n^2,{\sigma }_n^2>0.$$

$$p({\mathbf{a}}_{K_n}|{\mathbf{a}}_{K_{n-1}},{\sigma }_n^2,{\ddot{\mathbf{x}}}_n,{\mathbf{Y}}_{n-1},g_n)$$

$$=$$

$$p(g_n|\delta )=\frac{\delta -2}{2}(1+g_n{)}^{-\delta /2},g>0,$$

$$p(\mathcal{M}(n)|\mathcal{M}(n-1))$$

$$=$$

$$p(u_n=1|u_{n-1}=1).$$

$$p({\omega }_n,K_n|{\omega }_{n-1},K_{n-1},u_n=1,u_{n-1}=1)$$

$$=$$

$$p(K_n|K_{n-1},u_n=1,u_{n-1}=1),$$

$$p(\mathcal{M}(n)|{\mathcal{M}}_0(n-1))$$

$$=$$

$$p({\omega }_m,K_m|{\mathbf{Y}}_m,u_m=1)=\frac{p({\omega }_m,K_m,u_m=1|{\mathbf{Y}}_m)}{1-p(u_m=0|{\mathbf{Y}}_m)},$$

$$p({\sigma }_n^2|{\ddot{\mathbf{x}}}_n,{\mathbf{Y}}_{n-1})=p({\sigma }_n^2|{\sigma }_{n-1}^2,{\ddot{\mathbf{x}}}_n,{\mathbf{Y}}_{n-1}),$$

$$p({\mathbf{a}}_{K_n}|{\sigma }_n^2,{\ddot{\mathbf{x}}}_n,{\mathbf{Y}}_{n-1})$$

$$=$$

Full text

I Introduction

The problem of estimating the fundamental frequency or pitch information from noisy sound signals occurs in many applications, such as speech synthesis [1], voice disorder detection [2], and automatic speech recognition [3]. Fundamental frequency is a physical feature defined as the lowest frequency component of a periodic signal, while pitch is a perceptual feature, related to human listening [4]. Our objective is to estimate fundamental frequency. But, following [5, 6], we do not distinguish between fundamental frequency and pitch and use them interchangeably. Pitch is usually estimated using a segment of sound signals (a.k.a., frame) with a fixed segment length (e.g., 15-35 ms for speech signals [7]). Numerous pitch estimation algorithms have been proposed in the last fifty years, which can be categorized as unsupervised and supervised approaches. Unsupervised pitch estimation methods can be further categorized as non-parametric and parametric approaches. Examples of non-parametric approaches include the YIN [8], RAPT [9], SWIPE [10] and PEFAC [5] methods. YIN and RAPT compute autocorrelation functions from short frames of sound signals in the time domain. However, they are not robust against noise [11] and suffer from pitch octave errors (that is, a rational multiple of the true pitch) [12]. To reduce the pitch octave errors, SWIPE uses the cross-correlation function against a sawtooth signal combined with the spectrum of the signal, and exploits only the first and prime harmonics of the signal. PEFAC estimates the pitch in the log-frequency domain by convolving each frame’s power spectrum with a filter that sums the energy of the pitch harmonics. Dynamic programming is used to obtain a smooth estimate of the pitch track. Due to the filtering and built-in spectral normalization methods, PEFAC is claimed to work in high levels of noise. However, a long frame length (e.g., 90.5 ms in PEFAC by default) is required to obtain good pitch estimation accuracy which is not practical in many real-time applications. More recently, a single frequency filtering approach based pitch estimation algorithm is proposed, which exploits the high SNR frequency component to overcome the effects of degradations in speech signal [13].

By contrast, parametric methods (e.g., harmonic model-based pitch estimators [14, 6, 15]) have also been proposed for pitch estimation. Compared with non-parametric approaches, harmonic model-based pitch estimators work with a short frame length (e.g., 20 ms), and show higher robustness to additive noise, fewer octave errors, and better time-frequency resolution [16, 17]. Recently, a computationally efficient pitch estimator based on a harmonic model has been proposed, which is referred to as the fast NLS [11]. However, one problem with most of the harmonic model based pitch estimators is that they do not take the temporal smoothness of the pitch, the harmonic order, and voicing into account as they process each frame independently. As a result, outliers, due to octave errors or voicing detection errors, occur. A sample-by-sample Kalman filtering-based pitch tracking algorithm using a time-varying harmonic model is proposed in [18] by assuming that the pitch and weights follow first-order Markov chains. A particle filtering-based pitch tracking algorithm based on the source-filter speech model combining with the harmonic modelling of input source is introduced in [19]. However, the good performance of the algorithms in [18] and [19] requires careful initializations. Moreover, it is difficult to integrate the time-varying model order into these algorithms, see [20] as an example of combing discrete and continuous state spaces. With either a known or estimated model order, a maximum a posteriori (MAP) pitch estimator based on the harmonic model has been developed to exploit the temporal dynamics of the pitch [21]. The model weights and observation noise variance are estimated by maximizing the maximum likelihood function (i.e., a frequentist perspective). Smooth pitch estimates are obtained, and thus the pitch octave errors are reduced. An additional voicing state is further considered in [22] for estimating the pitch and obtaining the voiced-unvoiced decision jointly. However, the pitch tracking approach in [21] and [22] has many drawbacks. First, the assumption of a fixed harmonic order for multiple frames is not valid. In fact, in audio signals, the harmonic order often changes from frame to frame [23]. Second, matrix inversions are required to be stored for each candidate pitch to reduce the computational complexity. Third, errors can be found in transition frames where the voicing changes, because the past pitch information is not exploited when an unvoiced frame occurs. Finally, it is well-known that estimating parameters from a frequentist’s perspective leads to over-fitting [24].

More recently, neural network based supervised pitch tracking algorithms were proposed [25, 26, 27], which show robustness against noise. The method proposed in [26] uses deep stacking network for joint speech separation and pitch estimation. The CREPE [27] discretises the pitch in logarithmic scale and uses a deep convolutional neural network to produce a pitch estimate. However, the unvoiced/silent state is not considered in the model. The maximum value of the output of the neural network is used as a heuristic estimate of the voicing probability. Moreover, to satisfy user’s demand for different frequency resolution or frame length, the whole system is required to be retrained, which is usually time-consuming.

In this paper, we propose a fully Bayesian harmonic model-based pitch tracking approach. By using the harmonic model, as opposed to non-parametric methods, improved robustness against background noise and octave errors can be obtained. First-order Markov processes are used to capture the temporal dynamics of pitch, harmonic order, and voicing. By using information from previous frames, the rate of octave errors and the voicing detection errors can be further reduced. Compared to [21] and [22], we not only consider the temporal dynamics of pitch and voicing, but also of the harmonic order, which enables us to detect if any pitch is present, and estimate the pitch and harmonic order jointly and accurately. Moreover, past information on pitch is exploited to improve robustness against temporal voicing changes. Furthermore, by adopting a fully Bayesian approach to model weights and observation noise variances, the overfitting can be avoided. By assigning a proper transition pdf for the weights, fast NLS [11] can be easily incorporated into the proposed algorithm, leading to low computational and storage complexities.

The rest of the paper is organized as follows. In Section II, we briefly review general Bayesian tracking theory. In Section III and Section IV, we present the proposed harmonic observation and state evolution models, respectively. In Section V, the proposed pitch tracking algorithm is derived based on the harmonic observation and state evolution models. In Section VI, we briefly review the prewhitening step for dealing with non-Gaussian noise. Simulation results are given in Section VII, and the conclusions given in Section VIII.

Notation: Boldface symbols in lowercase and uppercase letters denote column vectors and matrices, respectively.

II Bayesian tracking

In this section, we briefly review Bayesian tracking in general, which forms the fundamental structure of the proposed pitch tracking algorithm. Consider the problem of estimating the state sequence $\{\mathbf{x}_{n}\},1\leq n\leq N$ from noisy observations $\{\mathbf{y}_{n}\},1\leq n\leq N$, related by

[TABLE]

where $h(\cdot)$ denotes a mapping function between the state and observation vectors, $\mathbf{v}_{n}$ denotes an i.i.d. observation noise sequence, and $n$ denotes the time index. The state sequence follows a first-order Markov process:

[TABLE]

where $f(\cdot)$ denotes a mapping function between the current and previous states, and $\mathbf{m}_{n}$ denotes an i.i.d. state noise sequence. The elements in the state vector $\mathbf{x}_{n}$ can either be continuous or discrete. Assume that the posterior pdf $p(\mathbf{x}_{n-1}|\mathbf{Y}_{n-1})$ is available with the initial pdf being defined as $p(\mathbf{x}_{0})$, where $\mathbf{Y}_{n-1}$ denotes a collection of observation vectors from the first observation vector up to the $({n-1})^{\mathrm{th}}$ observation vector, i.e.,

[TABLE]

The objective of Bayesian tracking is to obtain a posterior distribution over the state vector $\mathbf{x}_{n}$ based on the current and previous observations recursively, i.e., $p(\mathbf{x}_{n}|\mathbf{Y}_{n})$. The posterior $p(\mathbf{x}_{n}|\mathbf{Y}_{n})$ can be obtained in two stages: predict and update.

In the prediction stage, we obtain the prediction pdf $p(\mathbf{x}_{n}|\mathbf{Y}_{n-1})$ by using the transition pdf $p(\mathbf{x}_{n}|\mathbf{x}_{n-1})$ from (2), i.e.,

[TABLE]

which is known as the Chapman-Kolmogorov equation. Note that if the elements in $\mathbf{x}_{n}$ are all discrete variables, the integration operator should be replaced with the summation operator.

In the update stage, combining (1) and the prediction pdf from the prediction stage, Bayes’ rule can be applied to obtain the posterior, i.e.,

[TABLE]

where $p(\mathbf{y}_{n}|\mathbf{x}_{n},\mathbf{Y}_{n-1})$ and $p(\mathbf{y}_{1}|\mathbf{x}_{1})$ are the likelihood functions and $p(\mathbf{y}_{n}|\mathbf{Y}_{n-1})$ and $p(\mathbf{y}_{1})$ are the normalization factors, respectively. Closed form solutions can be obtained for (3) and (4) in at least two cases. In the first case, when both $\mathbf{v}_{n}$ and $\mathbf{m}_{n}$ are drawn from Gaussian distributions with known parameters, and both $h(\mathbf{x}_{n},\mathbf{v}_{n})$ and $f(\mathbf{x}_{n-1},\mathbf{m}_{n})$ are linear functions over the variables, (3) and (4) reduce to the well-known Kalman-filter[24]. In the second case, when the state space is discrete and has a limited number of states, (3) and (4) reduce to the forward step of the forward-backward algorithm for hidden Markov model (HMM) inference [24]. In other cases, the inference of the posterior $p(\mathbf{x}_{n}|\mathbf{Y}_{n})$ can be approximated using Monte Carlo approaches, such as particle filtering [28]. Next, we define the mapping function $h(\cdot)$ and formulate the observation equation (1) based on the harmonic model in Section III, and then explain the state evolution model (2) for the proposed pitch tracking algorithm in Section IV.

III Harmonic observation model

III-A The harmonic observation model

Consider the general signal observation model given by

[TABLE]

where the observation vector $\mathbf{y}_{n}$ is a collection of $M$ samples from the $n^{\mathrm{th}}$ frame defined as

[TABLE]

the clean signal vector $\mathbf{s}_{n}$ and noise vector $\mathbf{v}_{n}$ are defined similarly to $\mathbf{y}_{n}$, $M$ is the frame length in samples and $n$ is the frame index. We assume that $\mathbf{v}_{n}$ is a multivariate white noise processes with zero mean and diagonal covariance matrix $\sigma_{n}^{2}\mathbf{I}$, $\sigma_{n}^{2}$ is the noise variance, $\mathbf{I}$ is the identity matrix. When voiced speech or music is present, we assume that the pitch, model weights and model order are constant over a short frame (typically 15 to 35 ms for speech signals and longer for music signals) and $s_{n,m}$ (i.e., the $m^{\text{th}}$ element of $\mathbf{s}_{n}$) follows the harmonic model, i.e.,

[TABLE]

where $\alpha_{k,n}$ and $\beta_{k,n}$ are the linear weights of the $k^{\mathrm{th}}$ harmonic, $\omega_{n}={2\pi f_{n}}/{f_{\text{s}}}$ is the normalized digital radian frequency, $f_{\text{s}}$ is the sampling rate, and $K_{n}$ is the number of harmonics. When voiced speech/music is absent (unvoiced or silent), a null model is used, i.e.,

[TABLE]

Note that, based on the source-filtering model of speech generation, the unvoiced speech can be modelled as a coloured Gaussian process [29]. The observation noise in practice may have non-stationary and non-Gaussian properties, such as babble noise. However, we can deal with this by prewhitening the observation signals [23], which will be described in Section VI. Writing (6) in matrix form and combining (5) and (6) yields

[TABLE]

where

[TABLE]

We can further write (7) and (8) together by introducing a binary voicing indicator variable $u_{n}$, i.e.,

[TABLE]

where $u_{n}\in\{0,1\}$. When $u_{n}=0$ and $u_{n}=1$, (9) reduces to the unvoiced and voiced models (7) and (8), respectively.

We collect all the unknown variables into the state vector $\mathbf{x}_{n}=[\mathbf{a}_{K_{n}},\sigma_{n}^{2},\omega_{n},K_{n},u_{n}]^{T}$. Comparing (9) and (1), we can conclude that the mapping function $h(\cdot)$ is a nonlinear function w.r.t. the state vector $\mathbf{x}_{n}$. Moreover, the state vector $\mathbf{x}_{n}$ contains continuous variables $\mathbf{a}_{K_{n}}$, $\sigma_{n}^{2}$, $\omega_{n}$ and discrete variables $K_{n}$ and $u_{n}$. However, due to the non-linear characteristics of (9) w.r.t. $\omega_{n}$, uniform discretisation over the pitch $\omega_{n}$ is commonly used [11]. An off-grid estimate of $\omega_{n}$ can be obtained by pitch refinement algorithms, such as gradient descent [30]. Our target is to obtain estimates of the fundamental frequency $\omega_{n}$, the harmonic order $K_{n}$, and the voicing indicator $u_{n}$, that is a subset of $\mathbf{x}_{n}$ defined as $\ddot{\mathbf{x}}_{n}=[\omega_{n},K_{n},u_{n}]^{T}$, from the noisy observation $\mathbf{y}_{n}$.

IV The state evolution model

In this section, we derive the state evolution model (2) or more generally the transition probability density/mass function (pdf/pmf) $p(\mathbf{x}_{n}|\mathbf{x}_{n-1},\mathbf{Y}_{n-1})$ for continuous/discrete states of the proposed model. Following the fast NLS pitch estimation approach [11], we uniformly discretize the pitch $\omega_{n}\in\{\omega^{f},1\leq f\leq F\}$ over the range $[\omega_{\min},\omega_{\max}]$, where $\omega_{\min}$ and $\omega_{\max}$ denote the lowest and highest pitches in the searching space, respectively. Prior information can be used to set $\omega_{\min}$ and $\omega_{\max}$. For example, pitch is usually between 70 to 400 Hz for speech signals. The grid size is set to

[TABLE]

where $F$ denotes the DFT size for computing the likelihood function (see Section V and [11] for further details), $\left\lfloor\cdot\right\rfloor$ and $\left\lceil\cdot\right\rceil$ denote the flooring and ceiling operators, respectively. It is also shown that the optimal choice of $F$ depends on the frame length and the harmonic order [11]. However, for simplicity and fast implementation, in this paper, we set $F=2^{14}$. The state space for the discrete variables can be expressed as $\{\mathcal{M}(n):[\omega_{n}=\omega^{f},K_{n}=k,u_{n}=1]^{T},1\leq f\leq F,1\leq k\leq K^{\max}\}\cup\{\mathcal{M}_{0}(n):u_{n}=0\}$. The prediction pdf $p(\mathbf{x}_{n}|\mathbf{Y}_{n-1})$ defined in (3) can be factorized as

[TABLE]

We first explain the transition pdfs for the continuous variables $\sigma_{n}^{2}$ and $\mathbf{a}_{K_{n}}$, and then discuss the transition pmfs for the discrete variables $\omega_{n}$, $K_{n}$ and $u_{n}$. The selection of a state evolution model is a trade-off between being physically accurate and ending up with a practical solution.

IV-A Transition pdfs for the noise variance and weights

To obtain the prediction pdf for the noise variance $p(\sigma_{n}^{2}|\ddot{\mathbf{x}}_{n},\mathbf{Y}_{n-1})$, the transition pdf for the noise variance $p(\sigma_{n}^{2}|\sigma_{n-1}^{2},\ddot{\mathbf{x}}_{n},\mathbf{Y}_{n-1})$ should be defined. A reasonable assumption for the noise variance is that it changes slowly from frame to frame. For example, the unknown parameter $\sigma_{n}^{2}$ can be assumed to evolve according to an inverse Gamma distribution [31], i.e.

[TABLE]

where $\mathcal{IG}(x|\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{-\alpha-1}\exp(-\frac{\beta}{x})$ and $\Gamma(\cdot)$ denotes the gamma function. With this transition pdf, an analytical form of the posterior distribution on $\mathbf{x}_{n}$ cannot be derived. A sequential Monte Carlo approach can be used to approximate the posterior numerically [32]. However, the major drawback of any Monte Carlo filtering strategy is that sampling in high-dimensional spaces can be inefficient [33]. A Rao-blackwellized particle filtering approach [34], which marginalises out some of the variables for statistical variance reduction, can be used to deal with this problem. However, we do not pursue this approach any further in this paper, and leave it for future work. Instead, for simplicity, we assume independence between $\sigma_{n}^{2}$ and $\sigma_{n-1}^{2}$, and use the Jeffery’s prior, i.e.,

[TABLE]

As can be seen, the Jeffery’s prior (12) is a limiting case of (11) with $c\to 0$ and $d\to 0$.

Similarly, we define the transition pdf for the weights as $p(\mathbf{a}_{K_{n}}|\mathbf{a}_{K_{n-1}},\sigma_{n}^{2},\ddot{\mathbf{x}}_{n},\mathbf{Y}_{n-1})$. Imposing smoothness dependency on the weight time evolution can reduce pitch octave errors[35]. However, in order to use the fast algorithm [11], we assume that the model weights between consecutive frames are conditionally independent given previous observations and the rest of unknown variables. Following [36], we use the hierarchical prior

[TABLE]

where $\mathcal{N}(\mathbf{x}|\bm{\mu},\bm{\Sigma})$ denotes that the vector $\mathbf{x}$ has the multivariate normal distribution with mean $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$. The prior distribution for the weights (13) is known as Zellner’s g-prior [37]. As can be seen from (13), given $\omega_{n}$ and $K_{n}$, the prior covariance matrix is a scaled version of the Fisher information matrix. With Zellner’s g-prior, a closed-form calculation of the marginal likelihood can be obtained [38]. Moreover, the fast algorithm in [11] for computing the marginal likelihood can be applied (see Section V for detail).

The graphical model for the proposed method is shown in Fig. 1. Note that, instead of obtaining point estimates of the noise variance and weight parameters using maximum likelihood [21], a Bayesian approach is used to represent the full uncertainty over these parameters.

IV-B Transition pmfs for $\omega_{n},K_{n}$ and $u_{n}$

In [21], to reduce octave errors, a first-order Markov model is used for the pitch evolution provided that the harmonic order is fixed and known/estimated for multiple frames. Another voicing evolution model is further considered in [22] by imposing the so-called "hang-over" scheme [39]. Although in some cases, the harmonic order may not be of interest, it is still necessary to estimate it to obtain correct pitch estimates [40]. In fact, considering the temporal dynamics of the model order helps reducing the octave errors, which will be verified by the simulation results. Moreover, using priors for the model order is also necessary for model comparison [36]. In this paper, we propose to track the pitch $\omega_{n}$, the harmonic order $K_{n}$ and the voicing indicator $u_{n}$ jointly. More specifically, we impose smoothness constraints on $\omega_{n}$ and $K_{n}$, and hang-over on voicing state using first-order Markov processes. The transition probability for the $n^{\mathrm{th}}$ frame to be voiced with pitch $\omega_{n}$ and harmonic order $K_{n}$ when the previous frame is also voiced with $\omega_{n-1}$ and $K_{n-1}$ can be expressed as

[TABLE]

We assume that the pitch $\omega_{n}$ and harmonic order $K_{n}$ evolve according to their own, independent dynamics given $u_{n}=1$ and $u_{n-1}=1$, i.e.,

[TABLE]

which means when both time frame $n-1$ and $n$ are voiced, the pitch and harmonic order only depend on their previous states. In fact, this assumption is only true when the product of the maximum allowed harmonic order and the pitch is less than half of the sampling frequency. However, by using a Bayesian approach, a model with a larger harmonic order is more penalized than with a smaller one. Even if a large value is used for the maximum allowed harmonic order in the proposed approach, the posterior model probability with a large harmonic order can be small [41]. In [42], an infinite number of harmonics is used, and the non-parametric prior distribution is used to penalize the models with large harmonic orders. By assuming the pitch and harmonic order are conditionally independent given $u_{n}=1$ and $u_{n-1}=1$, the Bayesian inference of the model posterior, shown in Section V, can be simplified. The transition probability for the $n^{\mathrm{th}}$ frame to be voiced with pitch $\omega_{n}$ and harmonic order $K_{n}$ when the previous frame is unvoiced/silent can be expressed as

[TABLE]

The priors from an unvoiced frame to a voiced frame $p(\omega_{n},K_{n}|u_{n}=1,u_{n-1}=0)$ are set to $p(\omega_{m},K_{m}|\mathbf{Y}_{m},u_{m}=1)$, which can be calculated as

[TABLE]

where $m$ is the closest frame index to $n$ that satisfies the constraint $p(u_{m}=0|\mathbf{Y}_{m})<0.5$ ($m^{\mathrm{th}}$ frame is voiced). In fact, if the previous frame is not voiced, we exploit the information from the latest frame that is voiced as the prior for the pitches and harmonic orders. The motivation for this choice is that the pitch and harmonic order usually do not change abruptly after a short segment of unvoiced/silent frames. Using the past information as the prior, robustness against the voicing state changes can be improved. The graphical model for the evolution of $\ddot{\mathbf{x}}(n)$ is shown in Fig. 2. Assuming the Markov processes are time-invariant, we can express the transition matrices for the pitch, harmonic order and voicing as $\mathbf{A}^{\omega}$, $\mathbf{A}^{K}$ and $\mathbf{A}^{u}$, respectively.

V Pitch tracking

In this section, a joint pitch and harmonic order tracking, and voicing detection algorithm is derived based on the Bayesian tracking formulas (3) and (4). First, note that, by assuming that $\sigma_{n}^{2}$ and $\sigma_{n-1}^{2}$ are conditionally independent given $\ddot{\mathbf{x}}_{n}$ and $\mathbf{Y}_{n-1}$, and $\mathbf{a}_{K_{n}}$ and $\mathbf{a}_{K_{n-1}}$ are conditionally independent given $\sigma_{n}^{2}$, $\ddot{\mathbf{x}}_{n}$ and $\mathbf{Y}_{n-1}$, the prediction pdfs are equal to the transition pdfs, i.e.,

[TABLE]

Based on (3), prediction pmfs for discrete variables $p(\ddot{\mathbf{x}}(n)|\mathbf{Y}_{n-1})$ can be expressed as

[TABLE]

[TABLE]

With the prediction pdfs and pmfs in hand, we can obtain the update equation based on (4). In order to obtain the posteriors for the pitch, harmonic order and voicing indicators, the weights and noise variance can be integrated out from the update equation, i.e.,

[TABLE]

where $p(\mathbf{y}_{n}|\ddot{\mathbf{x}}_{n},\mathbf{Y}_{n-1})$ denotes a marginal likelihood, defined as

[TABLE]

Using (9), (12), (13), (14), (19) and (20), a closed-form marginal likelihood can be obtained, i.e.,

[TABLE]

where

[TABLE]

$m_{M}(\mathbf{y}_{n})$ denotes the null model likelihood (i.e., $p(\mathbf{y}_{n}|u_{n}=0)$) and ${}_{2}F_{1}$ denotes the Gaussian hypergeometric function [43]. To compute $R^{2}(\omega_{n},K_{n})$ for all the candidate pitches and harmonic orders, the fast algorithm [11] can be applied. Moreover, from a computational point of view, a Laplace approximation of (V) can be derived as an alternative instead of marginalizing w.r.t. $g_{n}$ analytically [36]. Note that, for the discrete vector $\ddot{\mathbf{x}}_{n}$, it should satisfy the normalisation constraint,

[TABLE]

Finally, estimates of the pitch and harmonic order and the voiced/unvoiced state can be jointly obtained using the maximum a posterior (MAP) estimator. More specifically, the $n^{\mathrm{th}}$ frame is labeled as voiced if $p(u_{n}=0|\mathbf{Y}_{n})<0.5$, and the pitch and harmonic order are obtained as

[TABLE]

The proposed Bayesian pitch tracking algorithm is shown in Algorithm 1. To make inferences, we need to specify the transition matrices for the pitch $p(\omega_{n}|\omega_{n-1},u_{n}=1,u_{n-1}=1)$, the harmonic order $p(K_{n}|K_{n-1},u_{n}=1,u_{n-1}=1)$ and $p(u_{n}|u_{n-1})$. Following [21], we set $p(\omega_{n}|\omega_{n-1},u_{n}=1,u_{n-1}=1)=\mathcal{N}(\omega_{n}|\omega_{n-1},\sigma_{\omega}^{2})$. The transition probability for the model order is chosen as $p(K_{n}|K_{n-1},u_{n}=1,u_{n-1}=1)=\mathcal{N}(K_{n}|K_{n-1},\sigma_{K}^{2})$. Smaller $\sigma_{\omega}^{2}$ and $\sigma_{K}^{2}$ lead to smoother estimates of the pitch and harmonic order while larger values make the inference less dependent on the previous estimates. The matrix $\mathbf{A}^{u}$ is controlled by $p(u_{n}=1|u_{n-1}=0)$ and $p(u_{n}=0|u_{n-1}=1)$. In order to reduce the false negative (wrongly classified as unvoiced when a frame is voiced) rate, we set $p(u_{n}=1|u_{n-1}=0)=0.4$, $p(u_{n}=0|u_{n-1}=1)=0.3$, respectively, that is, the transition probability from unvoiced to voiced is higher than from voiced to unvoiced. Note that each row of $\mathbf{A}^{\omega}$, $\mathbf{A}^{K}$, and $\mathbf{A}^{u}$ is normalised to ensure they are proper pmfs. By setting $\sigma_{\omega}^{2}\to\infty$, $\sigma_{K}^{2}\to\infty$, $p(u_{n}=1|u_{n-1}=0)=0.5$ and $p(u_{n}=0|u_{n-1}=1)=0.5$, the proposed algorithm reduces to the fast NLS algorithm [11].

VI Prewhitening

The fast NLS and proposed pitch tracking algorithm are derived under the assumption of white Gaussian noise. However, this assumption is usually violated in practice, for example, babble noise in a conference hall. Therefore, a prewhitening step is required to deal with the inconsistency between the white Gaussian noise model assumption and real life noise model. A linear prediction (LP) based prewhitening step is applied to each frame to deal with the non-white Gaussian noise (see [23, 44] for detail). The power spectral density (PSD) of the noise given noisy signals is estimated using the minimum mean-square error (MMSE) estimator [45]. We refer to the fast NLS and proposed algorithm with prewhitening step as Prew-Fast NLS and Prew-Proposed, respectively.

VII Simulation

In this section, we test the performance of the proposed harmonic model-based pitch tracking algorithm on real speech signals.

VII-A Databases

The databases used for evaluating the performance of different algorithms are as follows:

Keele pitch database: containing 10 spoken sentences from five male and five female speakers at a sampling rate of 20 kHz [46]. the "ground truth" pitch estimates are extracted from electroglottography with 10 ms time frame increment and 25.6 ms frame length. In fact, there are many spikes and wrong estimates in the "ground truth" pitch values, especially in the transient frames. However, we present the results for the Keele database to facilitate comparison with other pitch estimation algorithms that use this database.

Parkinson’s disease database: containing 130 sustained /a/ phonations from patients with Parkinson’s disease [47] at a sampling rate of 44.1 kHz. Each of the phonations is in one second length. The estimated "ground truth" pitches in 10 ms time frame increment are extracted from electroglottography (EGG).

VII-B Performance measures

Three performance measures are considered:

Total error ratio (TER) [22]: voicing detection performance measure. It is calculated based on the ratio between the number of incorrect voicing detection (false positive and true negative) estimates and the number of total estimates.

Gross error ratio (GER) [10]: accuracy measure of pitch estimates. It is computed based on the ratio between the number of pitch estimates that differ by more than 20 percents from the ground truth and the number of total estimates. The unvoiced frames from the ground truth are excluded and the pitch value of the voiced frame that is wrongly labeled as unvoiced frames by different pitch estimation algorithms is set to 0.

Mean absolute error (MAE) [47]: accuracy measure of pitch estimates. It is computed based on mean of the absolute errors between the ground truth and estimates. The unvoiced frames from the ground truth are excluded and the oracle voicing detector is used for all the algorithms.

VII-C Results

First, the proposed approach is tested on concatenated speech signals uttered by a female speaker first, male speaker second, sampled at 16 kHz. The spectrogram of the clean speech signals, pitch estimates, order estimates and the voicing detection results for PEFAC, CREPE, YIN, fast NLS and the proposed algorithm are shown in Fig. 3. The time frames of the spectrograms without red lines on top are unvoiced or silent frames. The variances for the transition matrices $\sigma_{\omega}^{2}$ and $\sigma_{K}^{2}$ are set to $\frac{16\pi^{2}}{f_{\text{s}}^{2}}$ and $1$, respectively. The SNR for white Gaussian noise is set to 0 dB. The candidate pitch $\omega_{0}$ is constrained to the range $2\pi\left[70\ 400\right]/f_{\text{s}}$ for PEFAC, YIN, fast NLS and the proposed algorithm. However, the results for the neural network based approach CREPE is based on the model with the pitch range $2\pi\left[32.7\ 1975.5\right]/f_{\text{s}}$ provided by the authors [27]. To change the settings for CREPE, re-training of the neural network model is required. The maximum allowed harmonic order for the proposed and fast NLS is set to 10. The frame length is $M=400$ samples (25 ms) with $60\%$ overlap. As can be seen from Fig. 3, the voicing detection results of both the fast NLS and the proposed algorithm are better than those of YIN, PEFAC and CREPE. For example, the frames around 2.8 s are correctly classified as voiced by the fast NLS and the proposed, but wrongly labeled as unvoiced by YIN, PEFAC and CREPE. Fast NLS suffers from octave errors, and has outliers particularly in the transition frames where voicing decisions change. In the transition frame around 1.8 s, the pitch and number of harmonics are wrongly estimated to 84.8 Hz and five, respectively, by the fast NLS. In contrast, they are estimated to 248.8 Hz and one, respectively, by the proposed. Clearly, the estimates of the proposed fit better to the spectrogram than the estimates of the fast NLS. The reason for the robustness against transient frames of the proposed algorithm is that the pitch and harmonic order information of the latest voiced frame is used as the prior, i.e. (18). The harmonic order of the frame in 2.2 s is estimated to two by both the fast NLS and the proposed. However, the pitch is wrongly estimated to 288.8 Hz by the fast NLS, but correctly estimated to 150.4 Hz by the proposed. By imposing temporal smoothness prior on the pitch using the Markov process model $p(\omega_{n}|\omega_{n-1},u_{n}=1,u_{n-1}=1)$, smoother estimates of the pitches are obtained. An octave error is produced by the fast NLS in the frame around 3.4 s. The pitch and harmonic order are estimated to 72 and six, respectively, by the fast NLS, but 143.2 and three, respectively, by the proposed. In fact, harmonic orders are estimated to three in the surrounding frames by both the fast NLS and the proposed. By using Bayesian tracking for the pitches and harmonic orders, smoother estimates of the pitches and harmonic orders are obtained. In conclusion, the proposed Bayesian pitch tracking algorithm obtains smooth estimates of the pitch and harmonic orders, and good voicing detection results by exploiting the past information.

The second experiment tests the performance of the proposed algorithm on musical instrument sounds (flute) from the University of Iowa Musical Instrument Samples database, decreasing from note B5 to C5. The spectrogram of the clean signals and the pitch estimates from fast NLS and the proposed algorithm are shown in Fig. 4. The music signal is downsampled to 16 kHz. The SNR for Gaussian white noise is set to -5 dB. The pitch $\omega_{0}$ is constrained to the range $2\pi\left[100\ 1500\right]/f_{\text{s}}$. The other parameters are set to the same as for Fig. 3. As can be seen, the proposed algorithm not only has smoother estimates of the pitch than fast NLS but also better voicing detection results.

In the third experiment, we test the performance of the proposed algorithm on the Keele pitch database with white Gaussian noise. Averages over 20 independent Monte Carlo experiments are used to obtain the experimental results. TER, GER and MAE for different SNRs for PEFAC, SWIPE, YIN, CREPE, fast NLS and the proposed algorithm are shown in Fig. 5, Fig. 6 and Fig. 7, respectively. For YIN, fast NLS and the proposed algorithm, the frame length is set to the same as the reference, i.e., 25.6 ms. Frame lengths 25.6 ms and 90.5 ms (default value) are used for PEFAC. The other parameters are set to the same as for Fig. 3. As can be seen from Fig. 5 and Fig. 6, PEFAC has better performance in terms of both GER and TER than CREPE at -10 dB SNR. Moreover, using a longer frame length (90.5 ms) for PEFAC leads to a lower GER but a higher TER compared with a shorter frame length (25.6 ms). SWIPE and YIN have similar performance in terms of TER and GER. The fast NLS method outperforms the PEFAC, SWIPE, YIN and CREPE. By imposing a smoothing prior on the pitches, harmonic orders and the voicing and using the harmonic model combined, the proposed algorithm achieves lower GER and TER than the fast-NLS. As can be seen from Fig. 7, when the oracle voicing detector is used, the SWIPE has the lowest MAE from 5 to 20 dB while the proposed algorithm achieves the best performance from -10 to 0 dB.

In the fourth experiment, the performance of the proposed algorithm with prewhitening is tested on the Keele pitch database in colored noise, i.e., babble noise and factory noise. The TER, GER and MAE results for Prew-proposed, Prew-fast NLS, PEFAC, Yin and SWIPE are shown in I, II and III, respectively. The linear prediction order for the prewhitening is set to 30. The maximum allowed harmonic order for the proposed and fast NLS is set to 30. The other parameters are set to the same as for Fig. 5. As can be seen from TABLE I and II, PEFAC with 90.5 ms and 25.6 ms have a lower TER and GER than YIN and SWIPE in -5 and 0 SNR conditions. The Prew-Proposed and Prew-Fast NLS have lower voicing detection errors and Gross errors than YIN, PEFAC and SWIPE in both babble and factory noise conditions. Although similar performance in term of TER can be seen for Prew-Proposed and Prew-Fast NLS, the Prew-Proposed has a lower GER than Prew-Fast NLS. As can be seen from TABLE III, when the oracle voicing detector is used, the SWIPE achieves the lowest MAE in babble noise. The Prew-proposed has a comparable performance with the SWIPE in babble noise and has the best performance in factory noise.

In the fifth experiment, we investigate the effect of reverberation on the performance of different pitch estimation algorithms. Reverberation is the process of multi-path propagation and occurs when the speech or audio signals are recorded in an acoustically enclosed space. A commonly used metric to measure the reverberation is the reverberation time (RT60) [48]. The reverberated signals used for testing are generated by filtering the signal by synthetic room impulse responses (RIRs) with RT60 varying from 0.2 to 1 s in 0.1 s step. The dimension of the room is set to $10\times 6\times 4$ m. The distance between the source and microphone is set to 1 m. The RIRs are generated using the image method [49] and implemented using the RIR Generator toolbox [50]. The position of the receiver is fixed while the position of the source is varied randomly from 60 degrees left of the receiver to 60 degrees right of the receiver for each Monte Carlo experiment. The TER, GER and MAE results on the Keele pitch database for the proposed, fast NLS, PEFAC, Yin and SWIPE are shown in Fig. 8, Fig. 9 and Fig. 10, respectively, where the parameters are set to the same as for Fig. 5. As can be seen from Fig. 8, the PEFAC (90.5 ms) has the lowest voicing detection errors in more reverberated conditions (RT60 from 0.5 to 1 s) while the proposed algorithm has a better voicing detection performance in less reverberated conditions. The proposed and the fast NLS has similar performance in terms of TER. However, as can be seen from Fig. 9, the proposed outperforms the PEFAC, SWIPE, CREPE, YIN and fast NLS in terms of GER. From Fig. 10, we can conclude that SWIPE has the best performance while the proposed is the second best one in terms of MAE.

In the final experiment, the performance of the proposed algorithm is tested on sustained /a/ signals (voiced) from the Parkinson’s disease database. The signals are downsampled to 16 kHz. TER, GER and MAE for different SNRs are shown in Fig. 11, Fig. 12 and Fig. 13, respectively. The parameters are set to the same as for Fig. 5. Similar conclusions to Fig. 5 and Fig. 6 can be drawn from Fig. 11 and Fig. 12. The proposed algorithm has the best performance in terms of the TER and GER. Moreover, the proposed algorithm presents the lowest MAE, especially from -15 to 0 dB. The spectrogram of one of the sustained /a/ sounds from the Parkinson’s disease database, pitch estimates of the PEFAC (oracle), YIN (oracle) and the proposed algorithm in 0 dB white Gaussian noise are shown in Fig. 14. The oracle voicing detector from the ground truth (all voiced) is used for both PEFAC and YIN. As can be seen from Fig. 14, the proposed algorithm outperforms the PEFAC (oracle) and YIN (oracle).

Based on the above experiments, PEFAC obtains a better pitch estimation and voicing detection performance than the neural network-based CREPE in low SNR scenarios. The SWIPE offers good performance in terms of MAE in high SNRs. The proposed algorithm obtains superior performance in terms of GER, TER and MAE compared to PEFAC, SWIPE,YIN, CREPE, and the fast NLS in low SNR scenarios for the Keele pitch database and Parkinson’s disease database.

VIII Conclusions

In this paper, a fully Bayesian harmonic model-based pitch tracking algorithm is proposed. Using a parametric harmonic model, the proposed algorithm shows good robustness against noise. The non-stationary evolution of the pitch, harmonic order and voicing state are modelled using first-order Markov chains. A fully Bayesian approach is applied for the noise variance and weights to avoid over-fitting. Using the hierarchical g-prior for the weights, the likelihood function can be easily evaluated using the fast NLS. The computational complexity of the recursive calculation of the predicted and posterior distributions is reduced by exploiting conditional independence between the pitch and harmonic order given the voicing indicators. Simulation results show that the proposed algorithm has good robustness against voicing state changes by carrying past information on pitch over the unvoiced/silent segments. The results of the pitch estimates and voicing detection for spoken sentences and sustained vowels are compared against ground truth estimates in the Keele and Parkinson’s disease databases, showing that the proposed algorithm presents good pitch estimation and voicing detection accuracy even in very noisy conditions (e.g., -15 dB).