Hidden Markov Convolutive Mixture Modelfor Pitch Contour Analysis of Speech
Kota Yoshizato1 , Hirokazu Kameoka1,2 , Daisuke Saito1 , Shigeki Sagayama1 ,1Graduate School of Information Science and Technology, The University of Tokyo, Japan
2
NTT Communication Science Laboratories, NTT Corporation, Japan{yoshizato,kameoka,dsaito,sagayama}@hil.t.u-tokyo.ac.jp

Abstract
This paper proposes a stochastic model of speech F0 contours, based on the stochastic formulation of the Fujisaki
model. Our motivation for the stochastic formulation is
twofold. Firstly, it allows us to derive a well-behaved
algorithm for estimating the Fujisaki model parametersfrom a raw F0 contour. Secondly, it will open the door
to incorporating the well-founded F0 contour model into
various statistical speech processing problems. We quantitatively evaluated the performance of our method interms of an Fujisaki-model parameter estimation accuracy using real speech data. Experimental results revealed that our method was superior to a state-of-the-art
Fujisaki model parameter extractor.
Index Terms: speech F0 contours, statistical model, Fujisaki model, hidden Markov model, EM algorithm

1. Introduction
The fundamental frequency (F0 ) contours in normal
speech contains various types of non-linguistic information such as speaker’s identity, emotion and attention. Modeling the F0 contours of speech utterances canthus be potentially very useful for many speech applications, including speech recognition, speaker recognition,
speech synthesis, and dialogue systems.
The Fujisaki model [1] is a well-founded mathematical model, which describes the generating process of thewhole F0 contour of a speech utterance. The remarkablefeature of the Fujisaki model is that it consists of physiologically and physically meaningful parameters (called
the phrase and accent commands) and is able to ﬁt F0
contours of real speech well when they are chosen appropriately. For this reason, the Fujisaki model has often been used to manually design an F0 contour for synthesizing natural speech. Thus, to enable speech synthesizers to automatically generate natural-sounding F0
contours, one way would be to incorporate the Fujisaki
model into the generative model of phonemic sequences
so that its parameters can be learned from a speech corpus in a uniﬁed manner. However, estimating (learning)
the Fujisaki model parameters from raw F0 contour observations has been a difﬁcult task. Several techniques
have already been developed [2, 3, 4], for the purposeof incorporating the extracted parameters in automatic
speech/emotion recognition systems to improve their per-formance, but so far with limited success due to the analytical complexity of the Fujisaki model.
We have previously derived a stochastic model ofspeech F0 contours by translating the Fujisaki model into
a probabilistic generative model [5, 6]. The stochastic
reformulation of the Fujisaki model has allowed us to derive an efﬁcient algorithm for automatically estimating
the Fujisaki-model parameters from raw F0 contours. In
this paper, we present an improved version of our previous models described in [5, 6]. The rest of this paper is organized as follows. Section 2 brieﬂy reviews the original
Fujisaki model. Section 3 formulates a probabilistic generative model of speech F0 contours based on the Fujisaki
model. Section 4 derives an algorithm for ﬁnding the
maximum a posteriori (MAP) estimates of the Fujisaki
model parameters from an observed F0 contour. Section
5 presents results of a quantitative evaluation conductedusing real speech data excerpted from the ATR speech
database. Section 6 concludes this paper.

2. Original Fujisaki Model
The Fujisaki model [1] assumes that an F0 contour on alogarithmic scale, y(t), where t is time, is the superposition of three components: a phrase component yp (t), an
accent component ya (t), and a base component yb :y(t) =yp (t) + ya (t) + yb .

(1)

The phrase component yp (t) consists of the major-scale
pitch variations over the duration of the prosodic units,
and the accent component ya (t) consists of the smallerscale pitch variations in accented syllables. These twocomponents are modeled as the outputs of second-order
critically damped ﬁlters, one being excited with a command function up (t) consisting of Dirac deltas (phrasecommands), and the other with ua (t) consisting of rectangular pulses (accent commands):
yp (t) =Gp (t) ∗ up (t),{ 2 −αt
α teGp (t) =0ya (t) =Ga (t) ∗ ua (t),{ 2 −βt
β te
Ga (t) =0(2)(t ≥ 0),(t < 0)(3)(4)(t ≥ 0),(t < 0)(5)sition. To parameterize the durations of the self transitions, each state is split into a certain number of substatessuch that they all have exactly the same emission densities. Fig. 2 shows an example of the splitting of statean . The number of substates is set at a sufﬁciently large
value and the transition probability from substate an,m tosubstate an,m+1 is set at 1 for m ̸= 0. This state splitting
allows us to ﬂexibly control the durations for which theprocess stays in state an through the settings of the transition probability. The transition probability from substate an,0 to substate an,m (m ≥ 1) corresponds to theprobability of the present HMM generating a rectangular pulse that has a particular duration. In the same way,
we split states p0 and a0 to parameterize the probability of the spacing between phrase and accent commands.
Henceforth, we use the notation p0 = {p0,0 , p0,1 , . . .},a0 = {a0,0 , a0,1 , . . .}, and an = {an,0 , an,1 , . . .}.
The present HMM is now deﬁned as follows:

Figure 1: Command function modeling with HMM.

Figure 2: The splitting of state an into 4 substates an,0 ,
an,1 , an,2 , and an,3 . ϕan,0 ,an,1 corresponds to the probability of staying at state an with 4 consecutive times.
where ∗ denotes convolution over time. The baseline
component yb is a constant value related to the lower
bound of the speaker’s F0 , below which no regular vocal fold vibration can be maintained. α and β are natural angular frequencies of the two second-order systems,
which are known to be almost constant within an utterance as well as across utterances for a particular speaker.
It has been shown that α = 3 rad/s and β = 20 rad/s canbe used as default values.

3. Stochastic model of speech F0 contours
Here, we model the generative process of an entire F0
contour of speech based on the discrete-time version of
the Fujisaki model.
We ﬁrst describe the process for generating the phrase
and accent command functions, up [k] and ua [k], where kdenotes the discrete-time index. In the original Fujisaki
model, it is required that the phrase commands must consist of Dirac deltas and the accent commands must consist of rectangular pulses. In addition, they are not allowed to overlap each other. To incorporate these requirements, we ﬁnd it convenient to model the up [k] and ua [k]pair, i.e., o[k] = (up [k], ua [k])T , using a Hidden Markov
Model (HMM). Speciﬁcally, we assume that {o[k]}K
k=1
is a sequence of outputs emitted from an HMM with thespeciﬁc topology illustrated in Fig. 1. The output distribution of each state is a Gaussian distribution
o[k] ∼N (o[k]; ν[k], Υ[k]) ,
[
]
[ 2
]
υp [k]
0µp [k]
ν[k] =, Υ[k] =,µa [k]0υa2 [k](6)(7)

where the mean vector ν[k] and variance matrix Υ[k] areconsidered to evolve in time as a result of the state tran-

Output sequence: {o[k]}K
k=1
Set of states: S = {p0 , p1 , a0 , . . . , aN }State sequence: {sk }K
k=1
Output distribution: P (o[k]|s8`k ) =´TN (o[k]; ‌[k], Υ[k])> 0, 0(sk ∈ p0 , a0 )<`´T‌[k] = Ap [k], 0 (sk = p1 )>:`(n) ´TAa(s–k ∈ an )
» 0,
2
0
σp,s
k
Υ[k] =
20σa,skTransition probability: ϕi′ ,i = log P (sk = i|sk−1 = i′ )Given the state sequence s = {sk }K
k=1 , the above
HMM generates the up [k] and ua [k] pair. From (2) and(4), up [k] and ua [k] are then fed through different critically damped ﬁlters, Gp [k] and Ga [k], to generate thephrase and accent components, xp [k] and xa [k]:xp [k] = up [k] ∗ Gp [k],xa [k] = ua [k] ∗ Ga [k],

(8)
(9)

where ∗ denotes convolution over k. The entire F0 contour is then given byx[k] = xp [k] + xa [k] + ub ,(10)

where ub denotes the baseline component. In non-tonal
languages such as standard Japanese, the pharse and accent commands should be non-negative. In our previous
model [5, 6], we treated up [k], ua [k] and ub as latent variables (i.e., parameters to be marginalized out), and did
not explicitly take the non-negativity constraints on up [k]and ua [k] into consideration. To explicitly impose thenon-negativity constraints, it is convenient to treat thesevariables as model parameters instead of latent variables.
For real speech F0 contours, observed F0 s should
not always be considered reliable. For example, F0 estimates obtained with a pitch extractor in unvoiced regions would be totally unreliable. When performing parameter inference, we would want to trust only reliable
observations and neglect unreliable ones. To incorporate
the degree of uncertainty of F0 observations, we consider
modeling an observed F0 contour y[k] as a superposition

of the “ideal”
x[k] and a noise component( F0 contour)
xn [k] ∼ N 0, υn2 [k] , where υn2 [k] represents the degreeof uncertainty of the F0 observation at time k, which isassumed to be given.
Overall, an observed F0 contour y[k] is described asy[k] = x[k] + xn [k]. For simplicity, we henceforth treat22, υn2 [k], α, and β as constants. By
ϕi′ ,i , ub , σp,i
, σa,imarginalizing xn [k] out, we obtain the probability densityKfunction of y = {y[k]}K
k=1 , given o = {o[k]}k=1 , asFigure 3: An example of command sequence matching.

K∏

P (y|o) =N (y[k]; x[k], υn2 [k]),k=1

x[k] = Gp [k] ∗ up [k] + Ga [k] ∗ ua [k] + ub .

(11)

Recall from (6) that given a state sequence s = {sk }K
k=1
(n) N
and θ = {{Ap [k]}Ka }n=1 , }, o is generated
k=1 , {A∏K
according to P (o|s, θ) =k=1 N (o[k]; ν[k], Υ[k]).
P (s) is given by the product of the state transition prob∏Kabilities: P (s) = ϕs1 k=2 ϕsk ,sk−1 . Furthermore, we
assume that θ is uniformly distributed.

4. Parameter Optimization Process
In this section, we describe an iterative algorithm thatsearches for the maximum a posteriori estimates of o andθ by locally maximizing P (o, θ|y) given y using thegeneralized Expectation-Maximization (EM) algorithm.
We treat s as a latent variable and consider marginalizing
P (o, θ, s|y) ∝ P (y|o)P (o|s, θ)P (s) with respect to s
to obtain the objective P (o, θ|y). The auxiliary function(as known as the “Q-function”) can be written as∑
Q(o, θ, o′ , θ ′ ) =P (s|y, o′ , θ ′ ) log P (o, θ, s|y)sc

= log P (y|o) +∑P (s|y, o′ , θ ′ ) log P (o|s, θ)P (s),scwhere = denotes equality up to constant terms.
An iterative algorithm that consists of computing
P (s|y, o′ , θ ′ ) (via the Forward-Backward algorithm), increasing Q(o, θ, o′ , θ ′ ) with respect to o and θ, and thensubstituting o and θ into o′ and θ ′ locally maximizes theposterior P (o, θ|y). Here, care must be taken that increasing Q(o, θ, o′ , θ ′ ) with respect to o must be performed subject to non-negativity. This can be done by
invoking the idea of [7]. By using the Jensen’s inequality
we obtain an inequality
−

(

∑

∑

i∈{p,a,b}l

Gi [k − l]ui [l]≥−)2∑ ∑ G2 [k − l]u2 [l]ii, (12)
λi,k,li∈{p,a,b} lwhere Gb [k] = δ[k] (Kronecker’s∑ ∑delta), λi,k,l ≥ 0 is an
auxiliary variable satisfying i l λi,k,l = 1. We canuse this inequality to construct a lower bound function

for Q(o, θ, o′ , θ ′ ). The maximization of this lower boundfunction w.r.t. o (subject to non-negativity) and λ can beachieved analytically, which guarantees a certain increaseof Q(o, θ, o′ , θ ′ ).
After convergence, we search for the optimal state sequence s by using the Viterbi algorithm.

5. Experiment
One important contribution of our work is that the Fujisaki model has successfully been translated into a statistical model. We believe that this will open the door to
combining our model and the various statistical speech
applications so that the Fujisaki-model parameters aswell as the spectral parameter sequences can be learned
from a speech corpus in a uniﬁed manner. In this regard, our model is already superior to conventional “nonstatistical” methods such as [4]. However, it is not yet
clear whether our statistical model is able to estimate theFujisaki model parameters from real speech data as accurately as the state-of-the-art technique. Thus, we quantitatively evaluated the parameter estimation accuracy ofthe present algorithm using real speech data, excerpted
from the ATR Japanese speech database B-set [8]. This
database consists of 503 phonetically balanced sentences.
We selected speech samples of one male speaker (MHT).
We used Fujisaki model parameters that had been manually annotated by an expert in the ﬁeld of speech prosody
as the ground truth data, where the baseline component
was set at log 60 Hz. F0 contours were extracted using
a method we had previously developed [9], from which
the Fujisaki model parameters were estimated using thepresent algorithm. The constant parameters were ﬁxed
respectively at N = 10, t0 = 8 ms, α = 3.0 rad/s, β =20.0 rad/s, υp2 [k] = 0.22 , υa2 [k] = 0.12 , υb2 = 0.0012 ,υn2 [k] = 1015 for unvoiced regions and υn2 [k] = 0.22 forvoiced regions. µb was set at the minimum log F0 value
in the voiced regions. The initial values of Θ were set
at the values obtained with the non-statistical method [4].
The EM algorithm was then run for 20 iterations. The
number of substates in the HMM and the transition probability ϕi′ ,i were determined according to the manuallyannotated data of the ﬁrst 200 sentences. The parameter estimation algorithm was then tested on the remaining
303 sentences. We evaluated the accuracy of the parameter estimation based on the following two criteria: log F0
RMSE (root mean squared error) and detection rates. Our
aim was to conﬁrm whether the present model and algorithm can achieve high model reconstruction accuracy

(a): The previous model(b): The proposed model

Figure 4: (1) An observed F0 contour in voiced regions (in solid line) and the estimated F0 contours (in dotted line) alongwith (2) the estimated phrase and accent commands.
Table 1: Detection rates and log F0 RMSE (S=0.3s).
Init
Estimated

Detection rates
0.688
0.695

log F0 [Hz] RMSE
0.17190.0611

while keeping the meaningfulness of the model parameters. log F0 RMSE was used to evaluate the reconstruction accuracy, which measures the root mean squared error between an observed F0 contour and the estimatedF0 contour. The detection rate was used to evaluate themeaningfulness of the parameter estimates, which was
calculated in the following way: We performed matching between the estimated and ground truth commandsequences as illustrated in Fig. 3 on a command-bycommand basis by using an dynamic programming algorithm. If the time difference between the estimatedand ground truth phrase commands was shorter than S
seconds, the estimated phrase command was considered“matched” and the local distance was set at zero. Otherwise the local distance was set at 1. As for the accentcommands, we took the average of the time difference
between the onsets of the estimated and ground truth accent commands and the time difference between the offsets of the estimated and ground truth accent commands.
In the same way, when the average time difference wasshorter than S seconds, the estimated accent command
was considered matched. The magnitudes of the phrase
and accent commands were not taken into account in ourevaluation. This is because the magnitude estimation wasvery sensitive to the baseline F0 value, which was set differently in the present method and in the manual annotation. Let NE , NA be the total numbers of commandsin the estimated and ground truth command sequences,
NM be the number of the matched commands between
the two sequences, NEsum , NAsum , and NMsum be thesum of NE , NA , NM for all 303 sentences. We deﬁned
the insertion error rate EI as (NEsum − NMsum )/NAsum ,the deletion error rate ED as (NAsum − NMsum )/NAsum ,and the detection rate D as 1 − EI − ED .
Tab. 1 shows the result of our quantitative evaluation with S = 0.3 s. The “Init” row shows the detection
rate and log F0 RMSE of the initial command sequence(which was obtained with the non-statistical method [4]),and the “Estimated” row shows that of the estimated com-mand sequence after the EM iterations. From the results,
we conﬁrmed that our method was comparable to a stateof-the-art Fujisaki model extractor in terms of the detection rate. On the other hand, our method was superior
to the conventional method in terms of the model reconstruction accuracy. We can also conﬁrm from Fig. 4 thatthe present model is able to ﬁt an observed F0 contour
more accurately than our previous model.

6. Conclusion
In this paper, we proposed a statistical model of speech
F0 contours and parameter estimation algorithm. We
evaluated the parameter estimation accuracy of the proposed method using real speech data, and conﬁrm the advantage of the proposed method. Future work will include incorporating the present model into the statistical speech applications such as the HMM-based speechsynthesis system (HTS) in such a way that the Fujisakimodel parameters can be learned from a speech corpus ina uniﬁed manner.

7. References[1] H. Fujisaki, In Vocal Physiology: Voice Production, Mechanisms
and Functions, Raven Press, 1988.
[2] H. Fujisaki and K. Hirose, “Analysis of voice fundamental frequency contours for declarative sentences of japanese,” J. Acoust.
Soc. Jpn (E), vol. 5, no. 4, pp. 233–242, 1984.
[3] H. Mixdorf, “A novel approach to the fully automatic extractionof fujisaki model parameters,” in Proc. ICASSP, 2000, vol. 3, pp.
1281–1284.
[4] S. Narusawa, N. Minematsu, K. Hirose, and H. Fujisaki, “A
method for automatic extraction of model parameters from fundamental frequency contours of speech,” in Proc. ICASSP, 2002,
pp. 509–512.
[5] H. Kameoka, J. L. Roux, and Y. Ohishi, “A statistical model ofspeech F0 contours,” in Proc. SAPA, 2010, pp. 43–48.
[6] K. Yoshizato, H. Kameoka, D. Saito, and S. Sagayama, “Statistical approach to fujisaki-model parameter estimation from speech
signals and its quantitative evaluation,” in Proc. Speech Prosody
2012, 2012, pp. 175–178.
[7] H. Kameoka, T. Nakatani, and T. Yoshioka, “Robust speech dereverberation based on non-negativity and sparse nature of speech
spectrograms,” in Proc. ICASSP, 2009, pp. 45–48.
[8] A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri, H. Kuwabara,
and K. Shikano, “ATR japanese speech database as a tool ofspeech recognition and synthesis,” Speech Communication, vol.
27, pp. 187–207, 1999.
[9] H. Kameoka, “Statistical speech spectrum model incorporating all-pole vocal tract model and F0 contour generating processmodel,” in Tech. Rep. IEICE, 2010, in Japanese.