INTERSPEECH 2013

Generative modeling of speech F0 contours
Hirokazu Kameoka1,2 , Kota Yoshizato1 , Tatsuma Ishihara1 ,Yasunori Ohishi2 , Kunio Kashino2 , Shigeki Sagayama1
12Graduate School of Information Science and Technology, The University of TokyoNTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation{kameoka,yoshizato,ishihara,sagayama}@hil.t.u-tokyo.ac.jp,{kameoka.hirokazu,ohishi.yasunori,kashino.kunio}@lab.ntt.co.jp

Abstract

the mechanical factors in the control mechanism of phonation,
the Fujisaki model has been widely used with notable successto design F0 contours for synthesizing natural speech. Sincea prosodic feature in speech is predominantly characterized by
the levels and timings of the phrase and accent components, one
important challenge is to solve an inverse problem of estimating the Fujisaki-model parameters automatically from a raw F0
contour.
However, this problem has been a difﬁcult task. Several
techniques have already been developed ([2–4], to name justa few), but so far with limited success due to the difﬁculty insearching for optimal parameters under the constraints imposedin the Fujisaki model. To overcome this difﬁculty, we havebeen concerned with translating the Fujisaki model into a probabilistic model such that one can make the best use of powerful
methods in statistical estimation theory (such as expectationmaximization algorithm and Viterbi algorithm) for the parameter estimation. The other important motivation for this formulation is to construct an automatically trainable version of theFujisaki model so that in future it can be used to develop a statistical speaking style conversion system or incorporated into
existing text-to-speech synthesis systems to improve the naturalness and intelligibility of computer-generated speech. This
paper introduces our ongoing work along with some new ideas
on generative modeling of speech fundamental frequency (F0 )contours based on a probabilistic reformulation of the Fujisaki
model [5–8]. We also brieﬂy introduce our recent work on generative modeling of F0 contours of singing voice [9] developed
under the same spirit.

This paper introduces our ongoing work on generative modeling of speech fundamental frequency (F0 ) contours for estimating prosodic features from raw speech data. The present F0
contour model is formulated by translating the Fujisaki model,
a well-founded mathematical model representing the control
mechanism of vocal fold vibration, into a probabilistic model
described as a discrete-time stochastic process. The motivation
behind this formulation is two fold. One is to derive a generalparameter estimation framework for the Fujisaki model, allowing for the introduction of powerful statistical methods. The
other is to construct an automatically trainable version of theFujisaki model so that in future it can be used to develop a statistical speaking style conversion system or incorporated into existing text-to-speech synthesis systems to improve the naturalness and intelligibility of computer-generated speech. We also
brieﬂy introduce a generative model of F0 contours of singing
voice developed under the same spirit.
Index Terms: speech F0 contour, Fujisaki model, generativemodel, hidden Markov model, EM algorithm

1. Introduction
Prosody aids the listener in interpreting an utterance by grouping words into larger information units and drawing attention
to speciﬁc words. It also plays an important role in conveying
various types of non-linguistic information such as the identity,
intention, attitude and mood of the speaker. Since the voice fundamental frequency (F0 ) contour is an important acoustic correlate of many prosodic constructs, modeling and analyzing F0
contours can be potentially useful for many speech applications
such as speech synthesis, speaker identiﬁcation, speech conversion and dialogue systems, in which prosodic information plays
a signiﬁcant role.
An F0 contour consists of long term pitch variations overthe duration of prosodic units and short term pitch variations inaccented syllables. The former usually contribute in phrasing
while the latter contribute in accentuation during an utterance.
These two types of pitch variations can be interpreted as the
manifestations of two independent movements by the thyroid
cartilage. The Fujisaki model [1] is a well-founded mathematical model that describes an F0 contour as the sum of these twocontributions. This model is known to approximate actual F0
contours of speech surprisingly well when the model parameters are chosen appropriately, and its validity has been shownfor many, typologically diverse languages. For this reason, andthanks to the intuitive association of the model parameters with

Copyright © 2013 ISCA

2. Original Fujisaki Model
The Fujisaki model [1], shown in Fig. 1, assumes that an F0
contour on a logarithmic scale, y(t), where t is time, is the superposition of three components: a phrase component xp (t), an
accent component xa (t), and a base component xb :y(t) =xp (t) + xa (t) + xb .

(1)

The phrase component xp (t) consists of the major-scale pitch
variations over the duration of the prosodic units, and the accent
component xa (t) consists of the smaller-scale pitch variationsin accented syllables. These two components are modeled as
the outputs of second-order critically damped ﬁlters, one being excited with a command function up (t) consisting of Diracdeltas (phrase commands), and the other with ua (t) consisting1826

25- 29 August 2013, Lyon, France

Figure 1: Original Fujisaki model [1].
of rectangular pulses (accent commands):
xp (t) =Gp (t) ∗ up (t),(α2 te−αt
Gp (t) =0xa (t) =Ga (t) ∗ ua (t),(β 2 te−βt
Ga (t) =0(2)(t ≥ 0),(t < 0)(3)
Figure 2: Command function modeling with HMM[5]. In state r0 ,(4)
(t ≥ 0),(t < 0)µp [k] and µa [k] are both constrained to be zero. In state p1 , µp [k]can take a non-zero value, Ap [k], whereas µa [k] is still restricted tozero. In state p1 , no self-transitions are allowed. In state r1 , µp [k]and µa [k] become zero again. This path constraint restricts µp [k] toconsisting of isolated deltas. State a0 leads to states a1 , . . . , aN , in(n)each of which µa [k] can take a different non-zero value Aa , whereasµp [k] is forced to be zero. Direct state transitions from state an to statean′ without passing through state r1 are not allowed. This constraint
restricts µa [k] to consisting of rectangular pulses.

(5)

where ∗ denotes convolution over time. The baseline component xb is a constant value related to the lower bound of the
speaker’s F0 , below which no regular vocal fold vibration canbe maintained. α and β are natural angular frequencies of thetwo 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 can be used as default values.

3. Discretized Fujisaki model
Figure 3: Duration-explicit representation of the hidden states [6].
The splitting of state an into substates an,0 , an,1 , an,2 , and an,3 allows us to parametrize the duration of each hidden state. For example,
ϕan,0 ,an,1 corresponds to the probability of staying at state an with 4
consecutive times.

In this section, we apply a backward difference s-to-z transform to the phrase and accent control mechanisms describedas continuous-time linear systems in order to obtain a discretetime version of the Fujisaki model [5]. The reason for the discretization will be made apparent later. The transfer function(Laplace transform) of the impulse response of the phrase control mechanism is given in the s-domain asˆ˜
Gp (s) = L Gp (t) =α2.
(s + α)24. Generative model of speech F0 contours
Here, we model the generative process of a speech F0 contour
based on the discrete-time version of the Fujisaki model.
We ﬁrst describe the process for generating the phrase andaccent command functions, up [k] and ua [k]. 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
proposed in [5] to model the up [k] and ua [k] pair, i.e., o[k] =(up [k], ua [k])T , using a hidden Markov model (HMM) with thespeciﬁc topology illustrated in Fig. 2. The output distribution ofeach state is a Gaussian distribution(6)

The backward difference transform approximates the time differential operator s by the backward difference operator in thez-domain such that s ≃ (1 − z −1 )/t0 , where t0 is the sampling
period of the discrete-time representation. By undertaking thistransform, the transfer function of the inverse system Gp−1 (s)can be written in the z-domain as
Gp−1 (z) = a2 z −2 + a1 z −1 + a0 ,(7)

with a2 = (ψ − 1) , a1 = −2ψ(ψ − 1), and a0 = ψ 2 ,where ψ = 1 + 1/(αt0 ). Let us use up [k] and xp [k] to denote the discrete-time version of the phrase command function
and phrase component, respectively, where k is the discretetime index. Then, xp [k] can thus be regarded as the output of aconstrained all-pole system whose characteristics are governedby a single parameter ψ (or α), such that2

up [k] = a0 xp [k] + a1 xp [k − 1] + a2 xp [k − 2].

o[k] ∼N (o[k]; ‌[k], Υ[k]) ,"#"µp [k]
υp2 [k]‌[k] =, Υ[k] =µa [k]
0(9)

with b2 = (φ − 1) , b1 = −2φ(φ − 1), and b0 = φ , whereφ = 1 + 1/(βt0 ). Altogether, the discrete-time version of
the Fujisaki model can be expressed as the superposition of thethree components: xp [k] + xa [k] + xb .
20.
υa2 [k](10)(11)

where the mean vector ‌[k] and covariance matrix Υ[k] evolvein time as a result of the state transition. To parameterize thedurations of the self transitions, we proposed in [6] to split eachstate into a certain number of substates such that they all haveexactly the same emission densities. Fig. 3 shows an example ofthe splitting of state an . The number of substates is set at a sufﬁciently large value and the transition probability from substatean,m to substate an,m+1 is set at 1 for m ̸= 0. This state splitting allows us to ﬂexibly control the durations for which the process stays in state an through the settings of the transition probability. The transition probability from substate an,0 to substate(8)

In the same way, the relationship between the accent command
function ua [k] and the accent component xa [k] is described asua [k] = b0 xa [k] + b1 xa [k − 1] + b2 xa [k − 2],

#

2

1827

(n)

N
s and „ = {{Ap [k]}K
k=1 , {Aa }n=1 } are the free parameters
to be estimated. Obviously, estimating s and „ corresponds
to estimating the command sequences, i.e., the Fujisaki model
parameters. To sum up, the likelihood function of the Fujisaki
model parameters s and „ given y is given as

an,m (m ≥ 1) corresponds to the probability of the present
HMM generating a rectangular pulse that has a particular duration. In the same way, we split states r0 and r1 to parameterize
the probability of the spacing between phrase and accent commands. Henceforth, we use the notation r0 = {r0,0 , r0,1 , . . .},r1 = {r1,0 , r1,1 , . . .}, and an = {an,0 , an,1 , . . .}. The HMM
is deﬁned as follows:

P (y|s, „) =Output sequence: {o[k]}K
k=1
Set of states: S = {r0 , p1 , r1 , a1 , . . . , aN }State sequence: {sk }K
k=1
Output distribution: P (o[k]|sk ) = N (o[k]; ‌[k], Υ[k])8` ´T>(sk ∈ r0 , r1 )<`0, 0´T‌[k] = Ap [k], 0 (sk = p1 )>:`(n) ´T(s#k ∈ an )" 0, Aa2σp,s0kΥ[k] =
2
0
σa,s
k

ȷ
ﬀ
|Σ−1 |1/2
1
T −1
exp
−
(y
−
—)
Σ
(y
−
—)
,
2
(2π)K/2— = A−1 —p + B −1 —a + xb 1,(16)` T ´−1` T ´−1−1−1
Σ = A Σp A
+ B Σa B+ Σn ,
whereΣp = diag(υp2 [1], . . . , υp2 [K]),(17)Σa = diag(υa2 [1], . . . , υa2 [K]),(18)Σn = diag(υn2 [1], . . . , υn2 [K]).

(19)

′

Transition probability: ϕi′ ,i = log P (sk = i|sk−1 = i )P (s) is given by the product of the state transition probabilities:
Q
P (s) = ϕs1 K
k=2 ϕsk ,sk−1 .
Readers are referred to [5–8] for detailed derivations of
the parameter inference algorithms: [5, 6] describe an iterativealgorithm for maximizing the posterior density P (s, „|y) ∝P (y, s, „) with respect to s and „. The deviced algorithm isbased on an expectation-maximization (EM) algorithm consisting in performing a Viterbi algorithm in the M-step. [7] describes a parameter estimation algorithm under non-negativity
constraints on the phrase and accent components.

In [8], we further proposed designing the transition networkof the HMM under the hypothesis that phrase and accent command sequences are governed by a vocabulary model. Refer to[8] for more details.
Given the state sequence s = {sk }K
k=1 , the above HMM
generates the up [k] and ua [k] pair. From Eq. (8) and Eq. (9),up [k] and ua [k] are then fed through different all-pole systems to generate the phrase and accent components, xp [k] andxa [k]. It should be noted that in non-tonal languages such as
standard Japanese, the phrase and accent commands must be
non-negative. The treatment of the non-negativity constraint on
xp [k] and xa [k] is discussed in [6, 7].
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, wewould 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” F0 contour,` i.e., x´p [k]+xa [k]+xb , and a noise component xn [k] ∼ N 0, υn2 [k] , where υn2 [k]represents the degree of uncertainty of the F0 observation attime k, which is assumed to be given. The entire F0 contour isthus given byy[k] = xp [k] + xa [k] + xb + xn [k],5. Singing voice F0 contour modeling
So far, we have focused on the modeling of F0 contours in normal speech. Here, we brieﬂy introduce our recent work on thegenerative modeling of singing voice F0 contours.
The F0 contours of a singing voice consist of two types ofdynamic components. One is called a note component, which isinﬂuenced by the physical constraints of the vocal folds, such as
overshoot, preparation and ﬁne ﬂuctuations [12, 13]. The otheris called an expression component, which corresponds to a mixture of singer’s musical expressive intentions, such as vibrato
and portamento [14]. Most previous papers have reported thatthese dynamic components strongly affect singing-voice perception, and that the former relates to the naturalness and individuality of a singing voice while the latter relates to singing
styles and skills [15]. Hence, automatic decomposition of a raw
F0 contour into these components can be potentially beneﬁcial
for many applications such as singer identiﬁcation, singing skill
evaluation, singing voice synthesis and singer style conversion.
Motivated by the above, we have been concerned with modifying the aforementioned F0 contour model so as to adapt tosinging voice F0 contours. The proposed model assumes thatthe process of generating singing voice F0 contours is structurally similar to the Fujisaki model (Fig. 4). The note and expression components are assumed to be the outputs of secondorder linear systems driven by the note and expression commands, corresponding to the musical note sequence and themusical expressive intentions, respectively. The note component contains the note transition and overshoot, whereas the expression component contains vibrato and portamento. The F0
contour is then modeled as a superposition of these two components. With the same strategy described in Sec. 4, we cantranslate this model into a probabilistic model and derive a parameter inference algorithm. For further details, refer to [9].

(12)where xb denotes the baseline component.
Now, let us deﬁne
up = (up [1], . . . , up [K])T , ua = (ua [1], . . . , ua [K])T ,—p = (µp [1], . . . , µp [K])T , —a = (µa [1], . . . , µa [K])T ,xp = (xp [1], . . . , xp [K])T , xa = (xa [1], . . . , xa [K])T ,xn = (xn [1], . . . , xn [K])T , y = (y[1], · · · , y[K])T .
Then, we can write up and ua as
up = Axp ,ua = Bxa ,

(13)
(14)

where

3
2
3
2
b0
O
a0
O
7
6b1 b07
6a1 a07
6
7
6
7
6b2 b1 a0
7
6a2 a1 a0
A=6
7.
7, B = 67
6 .. .. ..
7
6 .. .. ..
4
4. . . 5
. . . 5O
b2 b1 b0Oa 2 a1 a0

(15)

1828

Figure 5: An example of estimated phrase and accent commands (inblue and green, respectively), along with an observed F0 contour insolid red line and the optimized F0 contour model in dotted line [7].

Figure 4: The proposed model of a singing voice F0 contour [9].

6. Experimental evaluation
We quantitatively evaluated the parameter estimation accuracyof the present algorithm using real speech data, excerpted fromthe ATR Japanese speech database B-set [10]. This databaseconsists of 503 phonetically balanced sentences. We selectedspeech samples of one male speaker. We used Fujisaki model
parameters that had been manually annotated by an expert inthe ﬁ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[11], from which the Fujisaki model parameters were estimatedusing the present algorithm. The number of substates in the
HMM and the transition probability ϕi′ ,i were determined according to the manually annotated data of the ﬁrst 200 sentences. The parameter estimation algorithm was then testedon the remaining 303 sentences. We evaluated the accuracy ofthe 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 while
keeping the meaningfulness of the model parameters. log F0
RMSE was used to evaluate the reconstruction accuracy, whichmeasures the root mean squared error between an observed F0
contour and the estimated F0 contour. The detection rate wascalculated by performing matching between the estimated andground truth command sequences on a command-by-commandbasis by using a dynamic programming algorithm, where thetime difference between the estimated and ground truth commands shorter than 0.3 seconds was considered “matched” andthe local distance was set at zero (otherwise the local distance
was set at 1). 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 the sum 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 .
We chose Narusawa’s method [4] as a baseline method.
The present method obtained a detection rate of 69.5% while
the baseline method obtained 68.8%. This result conﬁrms that
our method was comparable to a state-of-the-art Fujisaki model
parameter extractor in terms of the detection rate. As for thelog F0 RMSE, on the other hand, the present method obtained
0.0611, while the baseline method obtained 0.1719. This result conﬁrms that our method was superior to the conventional
method in terms of the goodness-of-ﬁt property. Fig. 5 shows anexample of the estimated phrase and accent commands obtained
from a raw F0 contour.
To evaluate the pure behavior of the present parameter estimation algorithm, we also conducted a command estimation
experiment using a synthetic F0 contours. The synthetic F0

Figure 6: Estimation of note and expression commands [9].
contours were artiﬁcially created using the original Fujisaki
model with the abovementioned, manually annotated commandsequences. All other experimental conditions were the same asabove. The present method obtained a detection rate of 83.4%
while Narusawa’s method only obtained 72.6%. With this experiment, the present method was shown to be signiﬁcantly superior to the conventional method in terms of the detection rateof the command sequences.
Fig. 6 shows some examples of the decompositions of
singing voice F0 contours into note and expression components
by the method desribed in Sec. 5. As can be seen from these examples, pitch variations related to portamento and vibrato weresuccessfully separated from an observed F0 contour.

7. Conclusion
This paper introduced our ongoing work on generative modeling of F0 contours in speech and singing voice based on the
Fujisaki model. One important contribution of our work is that
the Fujisaki model has successfully been translated into an automatically trainable model. We believe that this will open thedoor to applying the present model to many speech applications such as speech synthesis, speaker identiﬁcation, speech
conversion, and dialogue systems, in which prosodic information plays a signiﬁcant role. Future work includes incorporating
the present model into the HMM-based speech synthesis system(HTS) [16] in such a way that the Fujisaki-model parameters
can be learned from a speech corpus in a uniﬁed manner.

8. Acknowledgements
We thank Prof. Keikichi Hirose (The University of Tokyo) for
kindly providing us with the manually annotated data associatedwith the ATR speech samples, Dr. Daisuke Saito (The University of Tokyo) and Dr. Daichi Mochihashi (The Institute of Statistical Mathematics) for fruitful discussions.

1829

9. 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. Le 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] K. Yoshizato, H. Kameoka, D. Saito, and S. Sagayama, “Hidden
Markov convolutive mixture model for pitch contour analysis of
speech,” in Proc. Interspeech 2012, 2012.
[8] T. Ishihara, H. Kameoka, K. Yoshizato, D. Saito, and S. Sagayama,“Probabilistic speech F0 contour model incorporating statistical
vocabulary model of phrase-accent command sequence,” submitted to Interspeech 2013.
[9] Y. Ohishi, H. Kameoka, D. Mochihashi, K. Kashino, “A stochastic
model of singing voice F0 contours for characterizing expressive
dynamic components,” in Proc. Interspeech 2012, 2012.
[10] A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri, H. Kuwabara,
and K. Shikano, “ATR japanese speech database as a tool of speech
recognition and synthesis,” Speech Communication, vol. 27, pp.
187–207, 1999.
[11] H. Kameoka, “Statistical speech spectrum model incorporating all-pole vocal tract model and F0 contour generating processmodel,” in Tech. Rep. IEICE, 2010, in Japanese.
[12] J. Sundberg, The science of the singing, Northern Illinois University Press, 1987.
[13] G. de Krom and G. Bloothooft, “Timing and accuracy of fundamental frequency changes in singing,” in Proc. ICPhS95.
[14] C. E. Seashore, “A musical ornament, the vibrato,” in Psychologyof Music, 1938, pp. 33–52.
[15] T. Saitou et al., “Acoustic and perceptual effects of vocal training
in amateur male singing,” in Proc. EUROSPEECH 2009.
[16] K. Tokuda, T. Yoshimura, T. Masuko, T. Kobayashi, T. Kitamura,“Speech parameter generation algorithms for HMM-based speechsynthesis,” in Proc. ICASSP, 2000, pp. 1315–1318.

1830