A Statistical Model of Speech F0 Contours
Hirokazu Kameoka, Jonathan Le Roux, Yasunori Ohishi1
1

NTT Communication Science Laboratories, NTT Corporation, Japan{kameoka,leroux,ohishi}@cs.brl.ntt.co.jp

Abstract
This paper proposes a statistical model of speech fundamental frequency (F0 ) contours, based on the formulation of thediscrete-time stochastic process version of the Fujisaki model,
which is known as a well-founded mathematical model representing the control mechanism of vocal fold vibration. There
are two important motivations for this statistical formulation.
One is to derive a general parameter estimation framework for
the Fujisaki model, allowing for the introduction of powerfulstatistical methods, and the other is to introduce a measure ofspeech naturalness in terms of an F0 contour through a probability distribution assumption, that can be incorporated into
many statistical speech processing problems such as speech
analysis, synthesis, separation, denoising and dereverberation.
Index Terms: speech F0 contour, statistical model

1. Introduction
The fundamental frequency (F0 ) contour in normal speech contains various types of information including emotional and othernon-linguistic information such as the speaker’s identity, moodand level of attention, and plays as important a role in our daily
speech communication as formants, through which we encodea phonemic sequence to convey linguistic information to thelistener(s). An F0 contour is a realization of the vocal fold oscillation with slowly varying frequencies, whose dynamics aregoverned by a combination of different factors, in particular thelength and elasticity of vocal folds, laryngeal muscle tension,and subglottal air pressure. All possible F0 contours produced
by a particular speech apparatus should thus be characterized
and constrained by the presumably small number of parameters governing the control mechanism of vocal fold vibration.
Therefore, how well the F0 contour of a certain sound matches
the mechanical constraint is an important factor that determines
how likely it is that the sound originates from a speech utterance. Accordingly, modeling the dynamics of the F0 contour ofspeech can be potentially very beneﬁcial for any speech applications that could be improved by taking account of the naturalness in terms of the F0 contour.
The Fujisaki model [1–3] is a well-founded mathematical model consisting of a set of physiologically and physicallymeaningful parameters, which describes the process of generating F0 contours by vocal folds in a reasonably simpliﬁedform. This model is known to approximate actual F0 contoursof speech surprisingly well when the model parameters are chosen appropriately, and its validity has been shown for many, typologically diverse languages. For this reason, and thanks to theintuitive association of the model parameters with the mechanical factors in the control mechanism of phonation, the Fujisaki
model has been widely used with notable success to design F0

contours for synthesizing natural speech. On the other hand,several techniques have been proposed for solving the inverse
problem of estimating the Fujisaki model parameters from raw
F0 contour observations [3–5], for the purpose of incorporating
the extracted parameters in automatic speech/emotion recognition systems in some way to improve their performance, but so
far with limited success due to the analytical complexity of theFujisaki model.
In this paper, we formulate a discrete-time stochastic process version of the Fujisaki model. Our motivation behind thisstatistical formulation is twofold. Firstly, by making the bestuse of statistical techniques, we expect to be able to derive apowerful framework for the estimation of the Fujisaki model
parameters, which has conventionally been considered a difﬁcult task. Secondly, it should enable us to represent the notionof speech naturalness in terms of the F0 contour through a probability distribution assumption. This will allow us to smoothly
incorporate our probabilistic model as an additional speech naturalness measure into many statistical speech processing problems such as speech separation, denoising and dereverberation.

2. Original Fujisaki model
The Fujisaki model [1–3], shown in Fig. 1, assumes that an F0
contour on a logarithmic scale, y(t), where t is time, is the superposition of two contributions associated with mutually independent types of movement of the thyroid cartilage with different degrees of freedom and muscular reaction times, referred
to as the phrase component, yp (t), and the accent component,ya (t), respectively [2]. The phrase component consists of the
major-scale pitch variations over the duration of the prosodic
units, which are characterized by a fast rise followed by a slower
fall. The accent component consists of the smaller-scale pitch
variations in accented syllables. These two components aremodeled as the outputs of second-order critically-damped ﬁlters, one being excited with Dirac deltas (phrase commands),
and the other with rectangular pulses (accent commands). The
linear systems producing yp (t) and ya (t), namely phrase control and accent control mechanisms, are characterized by(
α2 te−αt (t ≥ 0)
Gp (t) =,(1)
0(t < 0)(1 − (1 + βt)e−βt (t ≥ 0)
Sa (t) =,(2)0(t < 0)where Gp (t) is the impulse response of the phrase control
mechanism and Sa (t) is the step response of the accent control
mechanism. α and β are natural angular frequencies of the twosecond-order systems, which are known to be almost constantphrase command
functionwherephrasecontrolmechanism

accent command
functionaccentcontrolmechanism

Figure 1: Block diagram of Fujisaki model
within an utterance as well as across utterances for a particularspeaker. Although certain differences exist across speakers, it
has been shown that α = 3rad/s and β = 20rad/s can be used as
default values. The two components yp (t) and ya (t) are further
added by a constant value yb related to the lower bound for the
speaker’s F0 , below which no regular vocal fold vibration canbe maintained. The log F0 contour, y(t), is thus expressed asXy(t) = yb +Ap,i Gp (t − T0,i )
i

+

X

˘
¯
Aa,j Sa (t − T1,j ) − Sa (t − T2,j ) ,(3)jwhere T0,i and Ap,i denote the onset time and amplitude of thei-th phrase command, and T1,j , T2,j and Aa,j denote the onset
time, offset time and amplitude of the j-th accent command.
Although the two functions Gp (t) and Sa (t) look different,
it should be made clear that the impulse response of the accent
control mechanism, Ga (t), has the same form as the phrase control mechanism, Gp (t), namely(
β 2 te−βt (t ≥ 0)
Ga (t) =,(4)0(t < 0)since Ga (t) is deﬁned by the time differentiation of the step
response Sa (t). The accent component can thus be seen as the
convolution of Ga (t) and a command function composed of aset of stepwise functions.

3. Discretized Fujisaki model
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. 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 + α)2(5)The backward difference transform approximates the time differential operator s by the backward difference operator in thez-domain such thats≃

1 − z −1,t0(6)where t0 is the sampling period of the discrete-time representation. By undertaking this transform, the transfer function of theinverse system Gp−1 (s) can be written in the z-domain as
Hp−1 (z) = a2 z −2 + a1 z −1 + a0 ,(7)

a2 = (ψ − 1)2 ,a1 = −2ψ(ψ − 1),

(8)
(9)

a0 = ψ 2 ,(10)

and ψ = 1+1/(αt0 ). Here we use up [k] to denote the discretetime version of the phrase command function where k indicatesthe discrete-time index. The discrete-time version of the phrase
component, yp [k], can thus be regarded as the output of a constrained all-pole system whose characteristics are governed bya single parameter ψ (or α), such thatup [k] = a0 yp [k] + a1 yp [k − 1] + a2 yp [k − 2].

(11)

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

(13)
(14)

b0 = φ2 ,(15)

2where φ = 1 + 1/(βt0 ). For the baseline F0 , we use yb [k]to denote the discrete-time version of yb (t). Altogether, thediscrete-time version of the Fujisaki model can be expressed as
the superposition of the three components:
y[k] = yp [k] + ya [k] + yb [k].

(16)

4. Statistical formulation
4.1. Modeling phrase and accent command pair with HMM
We describe here how to model the command functions up [k]and ua [k]. To explain the phrase and accend commands in aphysiologically and linguistically meaningful way, it is important that they satisfy the following requirements:
1. Phrase commands are a set of impulses and accent commands are a set of step-wise functions.
2. A phrase command occurs at the start of an utterance orafter the offset of an accent command in the precedingphrase, and is followed by the onset of the next accentcommand. This means that a phrase command will notoccur while an accent command is being activated.
3. The onset of an accent command is followed by its offset.
This means that neighboring accent commands will notoverlap each other.
According to assumption 2, up [k] and ua [k] are reciprocallyconstrained and so they should not simply be modeled separately. One challenge in the estimation of the Fujisaki model
parameters has been how to deal with the optimization problem under these constraints. As a convenient way of incorporating these assumptions into the command functions, we propose modeling the up [k] and ua [k] pair using a hidden Markov
model (HMM). The use of an HMM for modeling the command functions is our primary reason for considering modelingthe discrete-time version of the Fujisaki model.
Let us arrange up [k] and ua [k] into a vector o[k]. We assume that o[k] is a random vector driven by an additive white
Gaussian noise (ϵp [k], ϵa [k])T such that
"
# "
# "
#
up [k]µp [k]
ϵp [k]
o[k] :==+.
(17)ua [k]µa [k]
ϵa [k]

Let ϵp [k] ∼ N (0, σp2 ) and ϵa [k] ∼ N (0, σa2 ) be mutually independent, theno[k] ∼ N (‌[k], Υ)
#
"
"
σp2
µp [k], Υ :=‌[k] :=µa [k]
0

0
σa2

#

p1

a0

p0

a1

(18)
.

(19)

a2

Equation (18) can be viewed as an HMM in which the output
distribution of each state is a Gaussian distribution. The mean
vector ‌[k] is thus considered to evolve in time as a result of thestate transition. This way of thinking allows us to incorporate
assumptions 1–3 into µp [k] and µa [k] by simply constraining
the path of the state transitions, as illustrated in Fig. 2.
The present HMM consists of N +3 distinct states, p0 , p1
and a0 , · · · , aN . In state p0 , µp [k] and µa [k] are both constrained to be zero. Given that the model is in state p0 , it is only
allowed either to stay in that state or to move to state p1 . In
state p1 , µp [k] can take a non-zero value as a function of time,
Ap [k], whereas µa [k] is still restricted to zero. It is important tonote that in state p1 , no self-transitions are allowed and only thetransition to state a0 is possible. In state a0 , µp [k] and µa [k] areagain both assumed to be zero. Hence, this speciﬁc path constraint restricts µp [k] to consisting of isolated deltas. State a0
leads to states a1 , · · · , aN , in each of which µa [k] can take a(n)different non-zero value Aa assumed to be constant in time,
whereas µp [k] is forced to be zero. Direct state transitions fromstate an to state an′ (n ̸= n′ , 1 ≤ n ≤ N , 1 ≤ n′ ≤ N ) without passing through state a0 are not allowed. This constraint
restricts µa [k] to consisting of rectangular pulses. It should also
be noted that the present HMM ensures that no more than one
command will be active at each point in time. To sum up, theproposed HMM is deﬁned as follows:
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]|sk = i) = N (ci [k], Υ)
8` ´T
0, 0
(i = p0 )>>><`A [k], 0´T (i = p )p1
ci [k] = ` ´T
>
0,
0
(i
=
a
0)
>
>
:`
(n) ´T0, Aa(i = an )Transition probability: ϕi′ ,i := log P (sk = i|sk−1 = i′ )For simplicity, we treat the transition probabilities ϕi′ ,i asconstant parameters in this paper, so that the free parameters tobe determined in our command function model consist of themagnitude of the phrase command, Ap [k], the state sequence,(n)sk , the magnitude of the accent command, {Aa }N
n=1 , and the22variance of the output distribution, σp , σa . Hereafter we use θu
to denote all these parameters:
˘
(n) N
2
2¯
θu := {Ap [k], s[k]}K(20)
k=1 , {Aa }n=1 , σp , σa .
Once the state sequence {sk }K
k=1 is speciﬁed, the mean seK
quences, {µp [k]}K
k=1 and {µa [k]}k=1 , namely the phrase andaccent command functions, are determined simultaneously by"#µp [k]= csk [k].
(21)
µa [k]

Profile ofa3

p0 p1 a0

a1

a0 a2

p0 p1 a0

state

Figure 2: Command function modeling with HMM.
4.2. Likelihood function and prior probabilities
In this subsection, we derive the probability density function of
the F0 contour, y[1], · · · , y[K], based on the statistical modeling of the command functions presented in the previous subsection. From Eqs. (18) and (19),up [k]|θu ∼ N (µp [k], σp2 ),(22)

ua [k]|θu ∼ N (µa [k], σa2 ).

(23)

Since up [k] and ua [k] are both assumed to be driven by white
noise, it can be written in vector notation as
up |θu ∼ N (—p , Σp ), Σp = σp2 I,(24)

ua |θu ∼ N (—a , Σa ), Σa = σa2 I,(25)

where up:=(up [1], · · · , up [K])T , ua:=(ua [1], · · · , ua [K])T ,—p := (µp [1], · · · , µp [K])T , —a :=(µa [1], · · · , µa [K])T . By
using the linear equation given in Section 3, the phrase component y p :=(yp [1], · · · , yp [K])T and the accent componenty a :=(ya [1], · · · , ya [K])T can be written in terms of up andua , respectively, such thatup = Ay p ,(26)

ua = By a ,

(27)

where
2
3
2
3
a0
O
b0
O
6a1 a0
7
6b1 b07
6
7
6
7
6a2 a1 a0
7
6b2 b1 a07
A :=67, B :=6
7. (28)6
7
6 .. .. ..
7.. .. ..
4
5
4
. . .
. . . 5Oa 2 a1 a0Ob 2 b1 b0
Hence, it follows from Eqs. (24) and (25) that`´
y p |θu , α ∼ N A−1 —p , A−1 Σp (A−1 )T ,`´
y a |θu , β ∼ N B −1 —a , B −1 Σa (B −1 )T .

(29)
(30)

As for the base component yb [k], we assume that it isalso a random variable driven by additive white Gaussian noiseϵb [k] ∼ N (0, σb2 ) such thatyb [k] = µb + ϵb [k],(31)y b |µb ∼ N (µb 1, Σb ), Σb = σb2 I.

(32)

and henceWe use θb := {µb , σb2 } to denote the parameters related to thebase component. Let ϵξ [j] and ϵξ′ [j ′ ] be mutually independentwhen (ξ, j) ̸= (ξ ′ , j ′ ), then y p , y a and y b are assumed to bestatistically independent. Care must be taken that this does notmean that —p and —a are independent. It therefore follows from
Eqs. (29), (30) and (32) that their sum y = y p + y a + y b will
also be normally distributed such that`y|Θ ∼ N A−1 —p + B −1 —a + µb 1,´
A−1 Σp (A−1 )T + B −1 Σa (B −1 )T + Σb ,

(33)

ȷ
ﬀ
|Σ−1 |1/2
1
T −1
exp
−
(y
−
—)
Σ
(y
−
—)
,
2
(2π)K/2−1y ← — + ΣM T (M ΣM T )−1 (y obs − M —),(36)

and in the maximization step, the model parameters are updated
using the data hypothetically completed in the expectation stepΘ ← argmax log P (y|Θ)P (Θ),

(37)

Θ

where Θ := {θu , ψ, φ, θb }. Overall, the likelihood function of
the Fujisaki model parameters Θ given y can be written as
P (y|Θ) =the complete data) are updated according to−1where P (y|Θ) and P (Θ) are deﬁned as Eqs. (34) and (35),respectively. Though trivial, when M = I for example, whichmeans there is no missing region, Eq. (36) becomes y ← y obs ,implying that we can treat the raw F0 data as is as the complete
data.
5.2. Proposed model as an F0 contour prior

5. Practical problem setting

Given some statistical model of speech features (e.g., waveform, spectrum or others) which takes F0 values as free parameters, let us consider the problem of estimating the unknown
parameters of the given model by incorporating the present statistical F0 contour model as a prior distribution over the F0 parameters. In this case the Fujisaki model parameters play therole of hyperparameters in the entire system under analysis.
Suppose that we are given a sequence of observed speechfeatures on the whole sample period, D = {dk }K
k=1 , whichis assumed to be generated according to a statistical model
P (D|y, Ξ), where y denotes the F0 parameters and Ξ contains all other parameters independent of y. By incorporating
Eq. (34) as a prior distribution over y, and assuming D is notdirectly dependent on the hyperparameter Θ, the Maximum A
Posteriori (MAP) estimation problem can be formalized as

In this section, we present some practical problems into which
the proposed statistical F0 contour model can be incorporated.

{ŷ, Θ̂, Ξ̂} = argmax log P (D|y, Ξ)P (y|Θ)P (Θ)P (Ξ).

—=A

−1

Σ=A

—p + B —a + µb 1,(34)` T ´−1` T ´−1−1
Σp A+ B Σa B+ Σb .

As for the prior probability of Θ, we assume that the parameters are independent of each other, the parameters other than
the state sequence, {s[k]}K
k=1 , phrase control parameter, ψ, andthe accent control parameter, φ, are uniformly distributed, and{s[k]}K
k=1 is a ﬁrst-order Markov chain:
P (Θ) ∝ P (ψ)P (φ)P (s1 )KYP (sk |sk−1 ).

(35)

k=2

y,Θ,Ξ

5.1. Parameter estimation from a raw F0 contour
The ﬁrst problem involves estimating the Fujisaki model parameters from a raw F0 contour. Here the raw F0 contour refers
to F0 data assumed to have been extracted using some F0 detection method from a speech signal of interest. Looking backat Eq. (34), it can be seen that we have thus far implicitly assumed that we are given a set of F0 observations on the wholesample period. However, F0 s should only exist in the voiced
region, which means that generally the observable data should
be “incomplete” in the sense that the F0 data in the unvoiced
regions are missing. Hence, when talking about the estimationof the Fujisaki model parameters from a raw F0 contour, it isgenerally necessary to deal with such incomplete data. In a statistical framework, this can simply be viewed as a missing data
imputation problem, which can be effectively dealt with using
the Expectation-Maximization (EM) algorithm [7].
Let us deﬁne y ∈ RK as the “complete” data, which consist of the observed raw F0 data and the missing data. If we
could estimate which region was the missing region in advance,
for example by using a voicing determination algorithm suchas the one incorporated in YIN [8], the relationship between′
the observed raw F0 data, y obs ∈ RK such that K ′ ≤ K,
and the complete data, y ∈ RK , could be written explicitly as
y obs = M y, where M is a K ′ -by-K binary matrix that hasexactly one entry 1 in each row and 0’s elsewhere. Owing tospace limitations, we omit all the mathematical details and only
provide the procedure needed in practice: in the expectation
step, the complete data y (strictly speaking, the expectation of

Although it may depend on how P (D|y, Ξ) is deﬁned, in manycases we need to employ an iterative method in which each iteration comprises the following conditional maximization steps:Ξ ← argmax log P (D|y, Ξ)P (Ξ),

(38)

Ξ

y ← argmax log P (D|y, Ξ)P (y|Θ),(39)Θ ← argmax log P (y|Θ)P (Θ).

(40)

y

Θ

Equation (39) updates the F0 contour estimate such that the statistical speech model best explains the observed speech features
under the prior distribution deﬁned by the Fujisaki model parameters obtained at the previous iteration. Eq. (40) then updates the Fujisaki model parameters such that the present statistical F0 contour model best explains the updated F0 contour
estimate y.

6. Parameter optimization process
We can notice that the algorithms in either of the two examples
shown in the previous section consist of performing the same
optimization step, that is, Eqs. (37) and (40). Here we describe
an iterative algorithm, which locally maximizes the posterior
density of Θ given y, P (Θ|y) ∝ P (y|Θ)P (Θ). By regarding a set consisting of the phrase, accent and base components,
T
T T
x := (y T
p , y a , y b ) , as the complete data, this problem can beviewed as an incomplete data problem, which can be dealt withagain using the EM algorithm. In this case, the log-likelihoodof Θ given the complete data is given as
1
1
log |Λ−1 | − (x − m)T Λ−1 (x − m),2
2
2 3
2 −1 3
yp
A —p
6 7
7
6
x := 4y a 5 , m := 4B −1 —a 5 ,(41)clog P (x|Θ) =yb
2 T −1A Σp A
6−1
Λ := 4
O

µb 1
O

B

O

T

O

3

7
O 5.

Σ−1
a B

1
δ1 (i) = − (o[1] − ci [1])T Υ−1 (o[1] − ci [1]).
2

for all the hidden states i. We can then compute δk (i) for k =2, · · · , K recursively via–
Υ−1 (o[k] − ci [k]) + ϕi′,i .

Taking the conditional expectation of Eq. (41) with respect to x
given y and Θ = Θ′ , and then adding log P (Θ) to both sides,
we obtain the Q function
1hQ(Θ, Θ ) =log |Λ−1 | − tr(Λ−1 E[xxT |y; Θ′ ])
2
i
′

c

+ 2mT Λ−1 E[x|y; Θ′ ] − mT Λ−1 m + log P (Θ). (42)

Because the relationship between the incomplete data and thecomplete data can be written as y = Hx where H :=[I, I, I], E[x|y; Θ] and E[xxT |y; Θ] are given explicitly as
E[x|y; Θ] = m + ΛH T (HΛH T )−1 (y − Hm), (43)
E[xxT |y; Θ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|y; Θ]E[x|y; Θ]T . (44)The expectation step computes E[x|y; Θ′ ] and E[xxT |y; Θ′ ]according to Eqs. (43) and (44) by substituting the current parameter estimate into Θ′ .
Now if we let E[x|y; Θ′ ] be partitioned into three K × 1
blocks and E[xxT |y; Θ′ ] into nine K × K blocks such that
2 3
3
2
x̄p
Rp ∗ ∗
6 7
7
6
E[x|y; Θ′ ] = 4x̄a 5, E[xxT |y; Θ′ ] = 4 ∗ Ra ∗ 5 , (45)∗x̄bwhere ∗ stands for blocks that we can ignore hereafter, then the
Q function can be rewritten in a more convenient form:
h
c 1
T −1
−1
log |AT Σ−1
Q(Θ, Θ′ ) =p A| + log |B Σa B| + log |Σb |2
T −1− tr(AT Σ−1
p ARp ) + 2—p Σp Ax̄p
T −1
− tr(B T Σ−1
a BRa ) + 2—a Σa B x̄a

−1
T −1
T −1
− —T
p Σp —p − —a Σa —a − —b Σb —b+ log P (Θ).

2) Magnitude of phrase command Ap [k]: Q(Θ, Θ) is maximized with respect to Ap [k] whenAp [k] = [Ax̄p ]k (k ∈ Tp1 ), Tp1 = {k|sk = p1 }.

(50)3) Magnitude of accent command Aa : Q(Θ, Θ′ ) is maxi(n)mized with respect to Aa when
(n)

A(n)
=
a

X
1
[B x̄a ]k , Tan = {k|sk = an }.
|Tan |

(51)

k∈Tan

4) Phrase control parameter ψ: Let us assume a Gaussianprior distribution over ψ such thatψ ∼ N (µψ , 1/νψ2 ).

(52)

Leaving only the terms in Q(Θ, Θ′ ) that depend on ψ, we have´
1 ` T −1
tr A Σp ARp
2
1 2
−1
+ —T
νψ (ψ − µψ )2 .
p Σp Ax̄p −2I2 (ψ) = log |A| −(53)

Now, leti(46)

The update formula for each parameter in the maximization step
can be derived using Eq. (46).
1) State sequence s1 , · · · , sK : Leaving only the terms inQ(Θ, Θ′ ) that depend on s := {sk }K
k=1 , we haveK

I1 (s) := −(49)

The most likely transition for each state should be registeredat each recursion Ψk (i) = argmaxi′ [δk−1 (i′ ) + ϕi′,i ], so thatthe most likely state sequence can be traced at the end of therecursion with sk−1 = Ψk (sk ) (k = K, · · · , 2), where sK =argmaxi δK (i). Substituting the updated state sequence {sk }into Eq. (21), we ﬁnally obtain the updated —p and —a .

∗ Rb

T −1
− tr(Σ−1b Rb ) + 2—b Σb x̄b

(48)

»
1
δk (i) = max
δk−1 (i′ ) − (o[k] − ci [k])T
i′
2

Σ−1
b

O

where o[k] := ([Ax̄p ]k , [B x̄a ]k )T . Here the notation [·]k isused to denote the k-th element of a vector. The state sequence{sk }K
k=1 maximizing I1 (s) can be solved efﬁciently using theViterbi algorithm as follows. We ﬁrst set δ1 (i) at

1X(o[k] − csk [k])T Υ−1 (o[k] − csk [k])
2

3
3
2
0
O
1
O
7
7
6 2 0
6−2 1
7
7
6
6
7,
7, U 1 := 6−2 2 0
1 −2 1U 2 := 6
7
7
6
6
.
.
.
.
.
.
.. .. .. 5.. .. .. 5
4
4
O
−2 2 0O
1 −2 1
3
2
0
O
7
60 0
7
6
7,
1 0 0
(54)
U 0 := 67
6 . . .
.. .. .. 5
4
O
1 0 0
2

k=1

+ log P (s1 ) +

K
X
k=2

then from Eqs. (8)–(10), A can be written aslog P (sk |sk−1 ),(47)

A = U 2 ψ2 + U 1 ψ + U 0 .

(55)

T
T
2 2
2
+ {tr((2U T
2 U 0 + U 1 U 1 )Rp ) − 2—p U 2 x̄p + σp νψ }ψ
2 2
2T+ {tr(U T
1 U 0 Rp ) − —p U 1 x̄p − 2σp νψ µψ }ψ − 2Kσp ,and its roots, namely the stationary points of Q(Θ, Θ′ ), can besolved algebraically, from which we can ﬁnd the optimal ψ.

F0 (Hz)Magnitude

4T
3
2tr(U T2 U 2 Rp )ψ + 3tr(U 2 U 1 Rp )ψ

Magnitude

The partial derivative of I2 (ψ) (or Q(Θ, Θ′ ) itself) with respectto ψ is a quartic function, equal up to a constant factor to

10

2

6KOGa U
0.5

06KOGa U

0.5

0

0

1

2

3

4

5
Time (s)

6

7

8

9

Figure 3: Test data: F0 contour in solid line and phrase component in dotted

7) Variances of additive noises σp2 , σa2 and σb2 : Q(Θ, Θ′ ) is
maximized with respect to σp2 , σa2 and σb2 when
` `
´
´
T
σp2 = tr AT ARp − 2—T
p Ax̄p + —p —p /K,
` `
´
´
T
σa2 = tr B T BRa − 2—T
a B x̄a + —a —a /K,`´
σb2 = tr(Rb ) − 2µb 1T x̄b /K + µ2b .

(56)
(57)
(58)

7. Experiment
Before applying our model to practical problems such as those
introduced in Section 5, we focus solely on testing the behavior of the optimization algorithm presented in Section 6. Thisis important because the devised algorithm is simply a local
search algorithm, which is not guaranteed to solve the global
optimum, and how seriously the local optima can affect the result of the parameter estimation is not yet understood. For thisreason, we only present the result of a numerical analysis that
we ran on artiﬁcial F0 contour data created using the original(continuous-time) Fujisaki model. There are two points that wewant to investigate using this experiment. One is the possibilityof avoiding undesirable solutions that may be caused by local
optima. The other is the inﬂuence of the approximation errors
originating from the discretization.
The test data, depicted in Fig. 3, were created with the following settings. We set the total duration at 10s, the sampling
period at 10ms, α = 3, β = 20, yb = 4, T0,1 = 1, T0,2 = 3,
T0,3 = 4.5, T0,4 = 6.5, Ap,1 = 0.7, Ap,2 = 0.6, Ap,3 = 0.3,
Ap,4 = 0.4, T1,1 = 1.5, T2,1 = 1.8, Aa,1 = 0.5, T1,2 = 2.1,
T2,2=2.5, Aa,2=0.4, T1,3=3.4, T2,3=3.8, Aa,3=0.6, T1,4=5.0,
T2,4=5.5, Aa,4=0.3, T1,5=6.8, T2,5=7.2, Aa,5=0.6, T1,6=7.5,
T2,6 =7.9, and Aa,6 =0.3. The conditions for the present algorithm were as follows. The iteration was run for 10 iterations.
N was set at 10, so that the total number of hidden states was
13. The state transition probabilities were set respectively at
ϕp0 ,p0 = log(0.999), ϕp0 ,p1 = log(0.001), ϕp1 ,a0 = log(1.0),
ϕa0 ,a0 =log(0.999), ϕan ,a0 =log(0.001), ϕa0 ,an =log(0.0001),
ϕan ,an =log(0.899), ϕan ,p0 =log(0.1), with 1 ≤ n ≤ 10.
The parameter estimates of —, —p and —a obtained using
the present algorithm are shown in Fig. 4. By comparing Fig. 4
with Fig. 3, one can conﬁrm that the parameter estimates are
extremely close to the correct values. Similar results were obtained when using initial parameters that were randomized differently. This fact indicates that we can solve the optimization
problem given by Eq. (37) or Eq. (40) quite satisfactorily withthe present algorithm, which strongly motivates us to deal with

Magnitude
Magnitude6) Mean of base component µb : Q(Θ, Θ′ ) is maximized withrespect to µb when µb = 1T x̄b /K.

F0 (Hz)line (top), phrase commands (middle), and accent commands (bottom).

5) Accent control parameter φ: Let us again assume a Gaussian prior distribution over φ such that φ ∼ N (µφ , 1/νφ2 ). As
the derivation follows in exactly the same manner as above, weshall omit it.

10

6KOGa U
0.5

06KOGa U

0.5

0

0

1

2

3

4

5
Time (s)

6

7

8

9

Figure 4: Parameter estimates: µ in solid line and A−1 µp+µb in dotted line(top), µp (middle), and µa (bottom).
the practical problems introduced in Section 5. Another important conclusion drawn from the results was that the approximation error due to the discretization did not appear to be a crucial
matter.

8. Conclusion
This paper proposed a statistical model of speech F0 contours,
based on the formulation of a discrete-time stochastic process
version of the Fujisaki model. Some examples of practical problems into which the proposed model could be incorporated werementioned. A parameter estimation framework for the proposedmodel based on the EM approach was derived, and the behaviorof the devised algorithm was tested on artiﬁcially created data.

9. References[1] H. Fujisaki, S. Nagashima, “A model for synthesis of pitch contours of connected speech,” Annual Report of Engineering Research Institute, University of Tokyo, Vol. 28, pp. 53–60, 1969.
[2] H. Fujisaki, “A note on the physiological and physical basis forthe phrase and accent components in the voice fundamental frequency contour,” In Vocal Physiology: Voice Production, Mechanisms and Functions, (O. Fujimura, ed.) Raven Press, pp. 347–
355, 1988.
[3] H. Fujisaki, 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.
[4] H. Mixdorf, “A novel approach to the fully automatic extractionof Fujisaki model parameters,” Proc. Intl. Conf. Acoust., Speech,
Signal Process. (ICASSP), Vol. 3, pp. 1281–1284, 2000.
[5] S. Narusawa, N. Minematsu, K. Hirose, H. Fujisaki, “A method
for automatic extraction of model parameters from fundamental
frequency contours of speech,” Proc. Intl. Conf. Acoust., Speech,
Signal Process. (ICASSP), Vol. 1, pp. 509–512, 2002.
[6] Y. Ohishi, H. Kameoka, K. Kashino, K. Takeda, “Parameter estimation method of F0 control model for singing voices,” Proc. Intl.
Conf. Spoken Language Process. (ICSLP), pp. 139–142, 2008.
[7] A.P. Dempster, N.M. Laird, D.B. Rubin, “Maximum likelihoodfrom incomplete aata via the EM algorithm,” J. of Royal StatisticalSociety Series B, Vol. 39, pp. 1–38, 1977.
[8] A. de Cheveigné, H. Kawahara, “YIN, a fundamental frequency
estimator for speech and music,” J. Acoust. Soc. Am., Vol. 111,
No. 4, pp. 1917–1930, 2002.