STATISTICAL APPROACH TO FUJISAKI-MODEL PARAMETERESTIMATION FROM SPEECH SIGNALS AND ITS QUANTITATIVEEVALUATION
Kota Yoshizato† , Hirokazu Kameoka†‡ , Daisuke Saito† , Shigeki Sagayama† ,†Graduate School of Information Science and Technology, The University of Tokyo, Japan‡NTT Communication Science Laboratories, NTT Corporation, Japan{yoshizato,kameoka,dsaito,sagayama}@hil.t.u-tokyo.ac.jp

Abstract
We have previously proposed a statistical model of speech F0
contours, which is based on the discrete-time version of the Fujisaki model. One advantage of this model is that it allows us
to introduce statistical methods to learn the Fujisaki-model parameters from speech F0 contours. This paper proposes several
modiﬁcations to our previous model and parameter inference
algorithm, and quantitatively evaluates the performance of our
modiﬁed parameter inference algorithm.
Index Terms: Speech F0 contours, statistical model, Fujisaki
model, hidden Markov model, EM algorithm

1. Introduction
The fundamental frequency (F0 ) contours in speech containvarious types of non-linguistic information such as the speaker’s
identity, emotion and level of attention. They also indicate
intonation in pitch accent languages. Modeling F0 contoursof speech utterances is therefore potentially useful for many
speech applications, in particular speech synthesis.
Thanks to the increasing availability of speech databases,
speech synthesis systems based on statistical models such as
hidden Markov models (HMMs) have attracted particular attention in recent years. Since unnatural F0 contours result in a synthesis that sounds “emotionless” to human listeners, one of the
primary challenges in speech synthesis technology is to create anatural-sounding F0 contour for the utterance as a whole. The
weakness of the current statistical text-to-speech systems is thatthey do not satisfactorily represent the macroscopic variationsof F0 s in natural speech.
The Fujisaki model [1] is a well-founded mathematicalmodel, which describes the process by which the whole F0
contour of a speech utterance is generated. The notable feature of the Fujisaki model is that it consists of physiologicallyand physically meaningful parameters (called phrase and accentcommands) and is able to ﬁt F0 contours of real speech wellwhen they are chosen appropriately. Thus, one way of enablingstatistical speech synthesizers to generate natural sounding F0
contours would be to incorporate the Fujisaki model into thestatistical model so that its parameters can be learned from aspeech corpus in a uniﬁed manner. However, the Fujisaki model
has an analytically complex form, making it difﬁcult to incorporate it into statistical speech synthesis systems as is.
To this end, we have recently introduced a statistical modelof the discrete-time version of the Fujisaki model [2]. Thismodel makes the best use of powerful statistical methods suchas the HMM and the Expectation-Maximization (EM) algo-rithm to estimate the Fujisaki model parameters. The aim of
this paper is to (1) make several modiﬁcations to our previous
model and parameter estimation algorithm, and (2) conﬁrm theperformance of our model through quantitative evaluations of
its parameter estimation accuracy.
The rest of this paper is organized as follows. Section 2
brieﬂy reviews the original Fujisaki model. Section 3 reviews
our previously introduced model, i.e., a discrete counterpart of
the Fujisaki model and its statistical model formulation. Section 4 presents an algorithm for estimating the Fujisaki model
parameters from observed real speech F0 contour data. Section
5 shows results of a quantitative evaluation obtained by conducting an experiment using real speech data excerpted fromthe ATR speech database. Section 6 concludes this paper.

2. Original Fujisaki Model
The Fujisaki model [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 yp (t), an accent component ya (t), and a base component yb . The phrase component
consists of the major-scale pitch variations over the duration of
the prosodic units, and the accent component consists of thesmaller-scale 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 (phrase commands), andthe other with ua (t) consisting of rectangular pulses (accentcommands). The baseline component is a constant value related
to the lower bound of the speaker’s F0 , below which no regular vocal fold vibration can be maintained. The log F0 contour,
y(t), is thus expressed asy(t) =yp (t) + ya (t) + yb ,

(1)

where
yp (t) =Gp (t) ∗ up (t),(α2 te−αt
Gp (t) =0ya (t) =Ga (t) ∗ ua (t),(β 2 te−βt
Ga (t) =0(2)(t ≥ 0),(t < 0)(3)(4)(t ≥ 0).
(t < 0)(5)∗ denotes convolution over time. α and β are natural angularfrequencies of the two second-order systems, which are known

generated from an HMM with the speciﬁc topology illustrated
in Figure 1. The output distribution of each state is a Gaussian
distribution
o[k] ∼N (‌[k], Υ[k]) ,
»
–
» 2
µ [k]
υ [k]
‌[k] = p, Υ[k] = pµa [k]
0

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 ofstaying at state an with 4 consecutive times.

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. Statistical Speech F0 Contour Model
This section describes a statistical model based on a discretizedversion of the Fujisaki model.
3.1. Discretization
To obtain a discrete-time version of the Fujisaki model, we apply a backward difference s-to-z transform to the phrase andaccent control mechanisms. With this discretization, the relationships between the ﬁlter inputs and outputs are given as
up [k] =a0 yp [k] + a1 yp [k − 1] + a2 yp [k − 2],(6)

a2 = (ψ − 1)2 , a1 = −2ψ(ψ − 1), a0 = ψ 2 ,ua [k] =b0 ya [k] + b1 ya [k − 1] + b2 ya [k − 2],

(7)
(8)

b2 = (φ − 1)2 , b1 = −2φ(φ − 1), b0 = φ2 ,(9)where k is the discrete time index, yp [k], up [k], ya [k] and ua [k]are the discrete-time versions of the phrase component, phrasecommand function, accent component and accent command
function, respectively, ψ = 1 + 1/(αt0 ), φ = 1 + 1/(βt0 ),and t0 is the sampling period of the discrete-time representation.
3.2. Statistical formulation
In the original Fujisaki model, the phrase commands and accent commands are assumed to consist of Dirac deltas andrectangular pulses, respectively. 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–0,υa2 [k](10)(11)

where the mean vector ‌[k] and variance matrix Υ[k] are considered to evolve in time as a result of the state transition.
To parameterize the durations of the state transitions, eachstate is split into a certain number of substates such that they
all have exactly the same emission densities. Figure 2 shows anexample of the splitting of state an . The number of substates
is set at a sufﬁciently large value and the transition probability from substate an,l to substate an,l+1 is set at 1 for l ̸= 0.
This state splitting allows us to ﬂexibly control the durations
for which the process stays in state an through the settings ofthe transition probability. The transition probability from substate an,0 to substate an,l (l ≥ 1) corresponds to the probability of the present HMM generating a rectangular pulse thathas 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 thenotation p0 = {p0,0 , p0,1 , . . .}, a0 = {a0,0 , a0,1 , . . .}, andan = {an,0 , an,1 , . . .}. The state splitting described above is
one of the modiﬁcations we have made to our previous model[2]. The present HMM is now 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]|s8`k = ´i)T = N (‌[k], Υ[k])> 0, 0(i ∈ p0 , a0 )<`´T‌[k] = Ap [k], 0 (i = p1 )>:`(n) ´TAa» 0,– (i ∈ an )2σp,i0
Υ[k] =20
σa,iTransition probability: ϕi′ ,i = log P (sk = i|sk−1 = i′ )We further assume that the baseline`component´ is normally distributed with mean µb , yb [k] ∼ N µb , υb2 .
For real speech F0 contours, we must take account of
the uncertainty in the observed F0 data, since observed data
should not always be considered reliable. For example, F0 estimates in unvoiced regions would be unreliable. Thus, we consider y[k] to be the superposition of a latent component corresponding to` the true
´ F0 and a noise component such thatyn [k] ∼ N 0, υn2 [k] . The variance υn2 [k] corresponds to thedegree of uncertainty of the observed F0 at time k. Therefore,y[k] = yp [k] + ya [k] + yb [k] + yn [k].
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 ,y p = (yp [1], . . . , yp [K])T , y a = (ya [1], . . . , ya [K])T ,y b = (yb [1], . . . , yb [K])T , y n = (yn [1], . . . , yn [K])T ,y = (y[1], · · · , y[K])T .
Then, we can write up and ua as
up = Ay p ,(12)ua = By a ,(13)be partitioned into four K × 1 blocks such that E[x|y; Θ] =(x̄Tp , x̄Ta , x̄Tb , x̄Tn )T .
The M step update formulas are given as follows.

where
2
3
2
3
b0
O
a0
O
6a1 a0
7
6b1 b07
7
6b2 b1 a0
7
6
1 a0
A = 6a2 a
,
B
=
7
6
7.
4 ... ... ... 5
4 ... ... ... 5Oa 2 a1 a0
O
b2 b1 b0

(14)1) State sequence: Leaving only the terms in Q(Θ, Θ′ ) that
depend on s := {sk }K
k=1 , we haveK

2
2
For simplicity, we treat ϕi′ ,i , µb σp,i
, σa,i, υb2 , υn2 [k], α, β as(n) NKconstants and Θ = {{Ap [k], s[k]}k=1 , {Aa }n=1 } as the freeparameters to be estimated. To sum up, the likelihood functionof the Fujisaki model parameters Θ given y is given asȷﬀ|Σ−1 |1/2
1
T −1
P (y|Θ) =
exp
−
(y
−
—)
Σ
(y
−
—)
,
2
(2π)K/2— = A−1 —p + B −1 —a + µb 1,(15)` T ´−1` T ´−1−1−1
Σ = A Σp A
+ B Σa B+ Σb + Σn ,2 2
υa [1]..
Σa = 4O
2 2
υn [1]..
Σn = 4OO.

υa2 [K]O.

υn2 [K]

3
3

5 . (16)

We present here an iterative algorithm that searches for the unknown parameters Θ by locally maximizing P (y|Θ) given y.
By regarding x = (y Tp , y Ta , y Tb , y Tn )T as the complete data thisproblem can be viewed as an incomplete data problem, which
can be dealt with using the EM algorithm. The log-likelihood
function of Θ given x is written as
1
1
log |Λ−1 | − (x − m)T Λ−1 (x − m),
22 3
22 −1 3
yp
A —p
6y 7
6B −1 — 7
a7
a7
6
x=6
,
m
=
(17)
4y b 54 µb 1 5 ,yn0
2 T −1
3
A Σp A
O
O
O
T −1
6
O
B Σa BOO 77.
Λ−1 = 6
−1
4
O
O
Σb
O 5
O
O
O
Σ−1
n
c

log P (x|Θ) =The auxiliary function is thus given byc1hlog |Λ−1 | − tr(Λ−1 E[xxT |y; Θ′ ])2i+ 2mT Λ−1 E[x|y; Θ′ ] − mT Λ−1 m + log P (Θ). (18)E[x|y; Θ] and E[xxT |y; Θ] are given explicitly as
E[x|y; Θ] = m + ΛH T (HΛH T )−1 (y − Hm), (19)E[xxT |y; Θ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|y; Θ]E[x|y; Θ]T ,k=1+ log P (s1 ) +K
Xlog P (sk |sk−1 ),

(21)

k=2

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 I(s) can be solved efﬁciently using theViterbi algorithm.

Ap [k] = [Ax̄p ]k (k ∈ Tp1 ), Tp1 = {k|sk = p1 },
X
1
A(n)
=
[B x̄a ]k , Tan = {k|sk = an }.
a|Tan |

5,

4. Parameter Optimization Process

Q(Θ, Θ′ ) =1X(o[k] − ‌[k])T Υ[k]−1 (o[k] − ‌[k])22) Magnitudes of phrase and accent commands: Q(Θ, Θ′ )(n)
is maximized with respect to Ap [k] and Aa when

where
3
2 2
υp [1]
O
..
5,
Σp =4
.
2
O
υ [K]
2 2
3p
υb
O
Σb =4 . . . 5 ,Oυb2I(s) := −(20)

by using the relationship y = Hx, where H = [I, I, I, I].
These are the values to be updated at the E step. Let E[x|y; Θ]

(22)
(23)

k∈Tan

The M step described above involves the computation ofAx̄p and B x̄a . This amounts to computing the ﬁlter inputsūp and ūa from the estimates of the phrase and accent components, x̄p and x̄a , using the “inverse ﬁltering” matrices A
and B. Note that as ūp and ūa correspond to the estimatesof command functions, they should be non-negative. However, the computations of Ax̄p and B x̄a allow each elementof the resulting vector to have a negative value. To effectively
avoid allowing such unwanted estimates, we consider ﬁndingūp and ūa that minimize ∥A−1 ūp −x̄p ∥22 and ∥B −1 ūa −x̄a ∥22
subject to non-negativity, instead of simply computing Ax̄pand
P B x̄a . In elementwise2 notation,
P they can be written as2
k |Gp [k] ∗ ūp [k] − x̄p [k]| andk |Ga [k] ∗ ūa [k] − x̄a [k]| ,where Gp [k] and Ga [k] are the discrete-time versions of the impulse responses of the phrase and accent control mechanisms.
Fortunately, this non-negative deconvolution problem can besolved efﬁciently by employing the method described in [3].

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 modeland the HMM-based text-to-speech synthesis system [4] so that
the Fujisaki-model parameters as well 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 “non-statistical” methods such as [5]. 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 [5]. Thus, we quantitatively
evaluated the parameter estimation accuracy of the present algorithm using real speech data, excerpted from the ATR Japanese
speech database B-set [6]. 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 speechprosody as the ground truth data, where the baseline component
was set constant at log 60 Hz.

Figure 3: An example of command sequence matching.
Table 1: Accuracy rate and Error rate (S=0.3s).
All Commands
Init
Estimated

A
0.688
0.697

EI
0.088
0.127

ED
0.224
0.177

Phrase Commands
Init
Estimated

A
0.647
0.680

EI
0.109
0.207

ED
0.244
0.112

Accent Commands
Init
Estimated

A
0.711
0.708

EI
0.076
0.083

ED
0.213
0.207

rate EI as (NEsum − NMsum )/NAsum , the deletion error rate
ED as (NAsum − NMsum )/NAsum , and the accuracy rate A as
1 − EI − ED . Table 1 shows the result of our quantitative evaluation with S = 0.3 s. The top, middle, and bottom tables show
the accuracy and error rates of the phrase and accent commands,
the phrase commands alone, and the accent commands alone,respectively. The “Init” row shows the accuracy and error ratesof the initial command sequence (which was obtained with themethod described in [5]), and the “Estimated” row shows thatof the estimated command sequence after the EM iterations. As
the result shows, our method performed slightly better than thestate-of-the-art technique [5].

6. Conclusion
In this paper, we made several modiﬁcations to our previously
reported model of speech F0 contours and parameter estimation algorithm. We evaluated the parameter estimation accuracy of the present method using real speech data. Future work
will include incorporating the present model into the HMMbased speech synthesis system (HTS) [4] in such a way that theFujisaki-model parameters can be learned from a speech corpus
in a uniﬁed manner. We are also currently exploring a reasonable model for the F0 contours of a singing voice [8].

7. Acknowledgement
F0 contours were extracted using the method described in[7], from which the Fujisaki model parameters were estimatedusing the present algorithm. The constant parameters were
ﬁxed respectively at N = 20, t0 = 8 ms, α = 3.0 rad/s,
β = 20.0 rad/s, υp2 [k] = 0.22 , υa2 [k] = 0.022 , υb2 = 0.0012 ,υn2 [k] = 1015 for unvoiced regions and υn2 [k] = 10−15 forvoiced regions. µb was set at the minimum log F0 value in thevoiced regions. The initial values of Θ were set at the values
obtained with the method described in [5]. The EM algorithm
was then run for 20 iterations. The number of substates in thepresent 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 parameter estimation accuracy in the following manner: We performed matching between the estimatedand ground truth command sequences as illustrated in Figure
3 on a command-by-command basis using dynamic programming. If the time difference between the estimated and groundtruth 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 wasset at 1. As for the accent commands, we took the average of thetime difference between the onsets of the estimated and groundtruth 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 was shorter
than S seconds, the estimated accent command was consideredmatched. The magnitudes of the phrase and accent commands
were not taken into account in our evaluation. This is because
the magnitude estimation was very 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 commands in the estimated and ground truth command sequences,
NM be the number of the matched commands between the twosequences, NEsum , NAsum , and NMsum be the sum of NE ,NA , NM for all 303 sentences. We deﬁned the insertion error

We thank Prof. Keikichi Hirose of the University of Tokyo, whokindly provided us with the manually annotated ground truth
data associated with the ATR speech samples.

8. References[1] H. Fujisaki, In Vocal Physiology: Voice Production, Mechanisms
and Functions, Raven Press, 1988.
[2] H. Kameoka, J. L. Roux, and Y. Ohishi, “A statistical model ofspeech F0 contours,” in Proc. ISCA Tutorial and Research Workshop on Statistical And Perceptual Audition (SAPA), 2010, pp. 43–
48.
[3] H. Kameoka, T. Nakatani, and T. Yoshioka, “Robust speech dereverberation based on non-negativity and sparse nature of speech
spectrograms,” in In Proc. 2008 IEEE Intl. Conf. Acoust., Speech,
Signal Process. (ICASSP 2009), 2009, pp. 45–48.
[4] K. Tokuda, T. Yoshimura, T. Masuko, T. Kobayashi, and T. Kitamura, “Speech parameter generation algorithms for HMM-based
speech synthesis,” in Proc. 2000 IEEE Intl. Conf. Acoust., Speech,
Signal Process. (ICASSP 2000), 2000, pp. 1315–1318.
[5] 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. 2002 IEEE Intl. Conf.
Acoust., Speech, Signal Process. (ICASSP 2002), 2002, pp. 509–
512.
[6] 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.
[7] H. Kameoka, “Statistical speech spectrum model incorporating allpole vocal tract model and F0 contour generating process model,”in Technical report of the Institute of Electronics, Information andCommunication Engineers (IEICE), 2010, in Japanese.
[8] Y. Ohishi, H. Kameoka, D. Mochihashi, H. Nagano, and
K. Kashino, “Statistical modeling of F0 dynamics in singing voices
based on Gaussian processes with multiple oscillation bases,” inProc. 11th Annual Conference of the International Speech Communication Association (INTERSPEECH 2010), 2010, pp. 2598–
2601.