1042

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 23, NO. 6, JUNE 2015

Generative Modeling of Voice Fundamental
Frequency Contours
Hirokazu Kameoka, Kota Yoshizato, Tatsuma Ishihara, Kento Kadowaki, Yasunori Ohishi, and Kunio Kashino

Abstract—This paper introduces a generative model of voice
fundamental frequency ( ) contours that allows us to extract
conprosodic features from raw speech data. The present
tour model is formulated by translating the Fujisaki model, awell-founded mathematical model representing the control mechanism of vocal fold vibration, into a probabilistic model described
as a discrete-time stochastic process. There are two motivations
behind this formulation. One is to derive a general parameter
estimation framework for the Fujisaki model that allows the introduction of powerful statistical methods. The other is to construct
an automatically trainable version of the Fujisaki model that
we can incorporate into statistical-model-based text-to-speechsynthesizers in such a way that the Fujisaki-model parameters
can be learned from a speech corpus in a uniﬁed manner. It could
also be useful for other speech applications such as emotion recognition, speaker identiﬁcation, speech conversion and dialoguesystems, in which prosodic information plays a signiﬁcant role. We
quantitatively evaluated the performance of the proposed Fujisaki
model parameter extractor using real speech data. Experimental
results revealed that our method was superior to a state-of-the-art
Fujisaki model parameter extractor.
Index Terms—Expectation-maximization algorithm, Fujisaki
model, prosody, voice fundamental frequency contour.

I. INTRODUCTIONPROSODY assists the listener to interpret an utterance bygrouping words into larger information units and drawing
attention to speciﬁc words. It also plays an important role inconveying various types of non-linguistic information such as
the identity, intention, attitude and mood of the speaker. Since
the voice fundamental frequency ( ) contour is an important

Manuscript received October 22, 2013; revised January 09, 2015; accepted
March 13, 2015. Date of publication April 01, 2015; date of current version
April 15, 2015. This work was supported by the JSPS KAKENHI under Grants
26730100 and 26280060. Earlier versions of this work were presented at SAPA2010 [1], Speech Prosody 2012 [2], Interspeech 2012 [3], and Interspeech 2013[4] as workshop/conference papers. The associate editor coordinating the reviewof this manuscript and approving it for publication was Dr. Nobutaka Ono.
H. Kameoka is with the Graduate School of Information Science and Technology, The University of Tokyo, 113-8656, Japan, and also with NTT Communication Science Laboratories, NTT Corporation, 243-0198, Japan (e-mail:kameoka@hil.t.u-tokyo.ac.jp; kameoka.hirokazu@lab.ntt.co.jp).
K. Yoshizato, T. Ishihara, and K. Kadowaki are with the Graduate Schoolof Information Science and Technology, The University of Tokyo, 113-8656,
Japan (e-mail: yoshizato@hil.t.u-tokyo.ac.jp; ishihara@hil.t.u-tokyo.ac.jp;
kadowaki@hil.t.u-tokyo.ac.jp).
Y. Ohishi and K. Kashino are with the NTT Communication Science Laboratories, NTT Corporation, 243-0198, Japan (e-mail: ohishi.yasunori@lab.ntt.
co.jp; kashino.kunio@lab.ntt.co.jp).
Color versions of one or more of the ﬁgures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identiﬁer 10.1109/TASLP.2015.2418576

acoustic correlate of many prosodic constructs, modeling andanalyzing
contours is 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. It is also important to notethat
contours indicate intonation in pitch accent languages.
Ancontour consists of long-term pitch variations overthe duration of prosodic units and short-term pitch variations inaccented syllables. The former usually contribute to phrasing
while the latter contribute to 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 [5], [6] is a well-founded mathematical model that describes ancontour as the sum ofthese two contributions. This model approximates actual
contours of speech fairly well when the model parameters areappropriately chosen, and its validity has been demonstrated
for many typologically diverse languages [6]–[13]. Since
prosodic features in speech are predominantly characterized by
the levels and timings of the phrase and accent components,
one important challenge is to solve the inverse problem of
estimating the Fujisaki-model parameters automatically from araw
contour.
However, this problem has proved difﬁcult to solve. Several
techniques have already been developed [6], [14]–[18], but so
far with limited success due to the difﬁculty of ﬁnding optimal
parameters under the constraints imposed in the Fujisaki model.
While the Fujisaki model describes a deterministic process forgenerating voicecontours, this paper proposes formulating astochastic counterpart of the Fujisaki model. As will be shownin the subsequent sections, this makes it possible to use powerfulstatistical inference techniques for estimating the underlying parameters of the Fujisaki model. Another important motivation
for this formulation is to construct an automatically trainableversion of the Fujisaki model that we can smoothly incorporate
into text-to-speech synthesis systems or speech conversion systems so as to guarantee the naturalness of computer-generated
speech.
The rest of this paper is organized as follows.
Section II brieﬂy reviews the original Fujisaki model.
Section III describes a discrete-time version of the Fujisaki
model, on which basis Section IV formulates a generativemodel of voicecontours. Section V presents two iterativealgorithms, which locally maximize the posterior density of
the Fujisaki model parameters given an observed
contour.
Section VI presents experimental evaluations of the present
method in terms of its ability as a Fujisaki model parameter
extractor.

2329-9290 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

KAMEOKA et al.: GENERATIVE MODELING OF VOICE FUNDAMENTAL FREQUENCY CONTOURS

1043

where is the sampling period of the discrete-time representation. By undertaking this transform, the transfer function of theinverse system
can be written in the -domain as
(8)
where

Fig. 1. Original Fujisaki model [5].

II. ORIGINAL FUJISAKI MODEL
The Fujisaki model [5], shown in Fig. 1, assumes that ancontour on a logarithmic scale,, where is time, is the superposition of three components: a phrase component, anaccent component, and a base value :(1)
The phrase and accent components are considered to be associated with mutually independent types of movement of the thyroid cartilage with different degrees of freedom and muscular reaction times. The phrase component
consists of the largescale pitch variations over the duration of the prosodic units,
and the accent component
consists of the smaller-scale
pitch variations in accented syllables. These two components
are modeled as the outputs of second-order critically damped
ﬁlters; one being excited with a command function
consisting of Dirac deltas (phrase commands), and the other withconsisting of rectangular pulses (accent commands):
(2)
(3)
where

denotes convolution over time and(4)(5)

is a constant value related to the lower bound of the, below which no regular vocal fold vibration
speaker’s
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 aparticular speaker. It has been shown thatrad/s andrad/s can be used as default values [6], [15].
III. DISCRETIZED FUJISAKI MODEL
In this section, we apply a backward difference -to- transform to the phrase and accent control mechanisms described ascontinuous-time linear systems in order to obtain a discrete-timeversion of the Fujisaki model. The reason for the discretization
will be made apparent later. The transfer function (the Laplacetransform of the impulse response) of the phrase control mechanism is given in the -domain as(6)
The backward difference transform approximates the time differential operator by the backward difference operator in the-domain such that
(7)

(9)
(10)
(11)
(12)
and
to denote the discrete-time version
Let us useof the phrase command function and phrase component, respectively, where is the discrete-time index.
can thus beregarded as the output of a constrained all-pole system whosecharacteristics are governed by a single parameter (or )(13)
In the same way, the relationship between the accent commandand the accent component
is described asfunction
(14)
where
(15)
(16)
(17)
(18)
Altogether, the discrete-time version of the Fujisaki model canbe expressed as the superposition of the three components:.
IV. GENERATIVE MODEL OF VOICECONTOURS
Here, we model the probabilistic generative process of aspeechcontour based on the discrete-time version of theFujisaki model.
A. Modeling Phrase and Accent Command Pair
We ﬁrst describe the process for generating the phrase andaccent command functions,and. In the original Fujisaki model, they must 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 accent command. This means that a phrase command will not occur
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,
andare reciprocallyconstrained and so they should not simply be modeled separately. One challenge as regards the estimation of the Fujisaki
model parameters has been how to deal with the optimization
problem under these constraints. As a convenient way ofincorporating these requirements into the command functions,1044

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 23, NO. 6, JUNE 2015

Fig. 3. A duration-explicit representation of the hidden states. Splitting state
into substates
,
,
, and
allows us to parametrize the duration
corresponds to the probability ofof each hidden state. For example,
with 4 consecutive times.
staying at stateFig. 2. Command function modeling with HMM.

we propose modeling theMarkov model (HMM).
We denote theandit is normally distributed:

andpair bypair using a hidden
and assume that
(19)

where

Eq. (19) can be viewed as an HMM in which the output distribution of each state is a Gaussian distribution and the mean vector
evolves over time as a result of the state transition. The
consists of the means of the phrase and acmean vector
cent command functions,and. The use of an HMM
allows us to incorporate assumptions 1–3 intoandby constraining the path of the state transitions as illustrated inFig. 2.
The present HMM consists ofdistinct states, , ,and.
is the number of possible values thatthe magnitude of each accent command can take. It can thus be
understood as the resolution of magnitude “quantization”: thelarger it becomes, the more ﬁnely the model is able to express
the accent command function. In state ,andareboth restricted to zero. In state ,can take a non-zerovalue,, whereasis still restricted to zero. In state, no self-transitions are allowed. In state ,andbecome zero again. This path constraint restricts
to con,sisting of isolated deltas. State leads to statesin each of which
can take a different non-zero value,
whereasis forced to be zero. Direct state transitions fromstate to state () without passing through state areto consisting of rectnot allowed. This constraint restricts
angular pulses. It should also be noted that this HMM ensuresthat no more than one command will be active at each point intime. The use of the HMM described above for modeling thecommand functions has been our primary reason for translating
the Fujisaki model into its discrete-time counterpart.
The state segments correspond to the timings and durationsof phrase and accent commands. If the statistical distributionsof the state durations can be trained a priori, they can be usefulfor estimating the timings of phrase and accent commands.
While an ordinary HMM assumes the state durations to begeometrically distributed, it would be more convenient if wewere allowed to assume arbitrary distributions. To allow arbitrary distributions, we propose splitting each state into a certain
number of substates such that they all have exactly the sameemission densities. Fig. 3 shows an example of the splittingof state. The number of substates is set at a sufﬁciently
large value and the transition probability from substate
to substateis set at 1 for. This state splitting
allows us to assume arbitrary distributions over the durations
for which the process stays in statethrough the settingsof the transition probability. The transition probability fromsubstate
to substatecorresponds to theprobability of the present HMM generating a rectangular pulsethat has a particular duration. In the same way, we split states
andto parameterize the probability of the spacing between
phrase and accent commands. Note that this is equivalent to
the explicit-duration HMM proposed by Ferguson [19]. Alternatively, the use of a hidden semi-Markov model [20], [21]would also be appropriate for the same purpose. Henceforth,
we use the notation
,
,
and
. Let
be the logarithm of thetransition probability from state and . To sum up, the present
HMM is deﬁned as follows:
Output sequence:
Set of states:
State sequence:
Output distribution:

Transition probability:

The free parameters to be estimated in our command function
model consist of the state sequence,, and the mean and,variance of the output distribution of each state,,,. Hereafter, we use to denoteand to denote the rest of the parameters:

KAMEOKA et al.: GENERATIVE MODELING OF VOICE FUNDAMENTAL FREQUENCY CONTOURS

The generating process for the phrase and accent components
is summarized as follows: First, the state sequenceis generated according to a Markov chain. Given the state
sequence, the mean sequences
andare determined by

1045

where
(30)
(31)
(32)
(33)

(20)
and
the present HMM generates theandpair according to Eq. (19). From Eq. (13) and Eq. (14),andare then fed through different all-pole systems
and.
to generate the phrase and accent components,
Given

B. Uncertainty of

By using the linear equation given by Eqs. (13) and (14), thevectors consisting of the phrase and accent components
(34)
(35)
and

can be written in terms ofObservations

contours and voiced/unIn the following, we assume thatvoiced (V/UV) segments are obtained in advance by using apitch extractor and a V/UV detector. For real speech,
valuesshould not always be considered reliable. For example,
estimates obtained with a pitch extractor in unvoiced regions would
be totally unreliable. When performing parameter inference, weshould trust only reliable observations and neglect unreliable
ones. To incorporate the degree of uncertainty ofobservacontour
as ations, we consider modeling an observed
superposition of the “ideal”
contour, i.e.,,
and a normally distributed noise component

,
(36)
(37)

where

..

.

..

.

..

.

..

.

..

(38)
.

(21)
represents the degree of uncertainty of theobwhereservation at time , which is assumed to be given. For example,
one simple way would be to set
at a small value near 0
for voiced regions and a sufﬁciently large value for unvoiced
regions. By denoting

..

(39)
.

Hence, it follows from Eqs. (28), (29), (36) and (37) that
(40)

(22)
the entire

contour is given by(41)
We refer toas the base component and let

be

(23)

(42)
It follows from Eqs. (21) and (22) that

C. Likelihood Function and Prior Probabilities
In this subsection, we derive the probability density functionof an observed
contour,, based on theprobabilistic modeling of the command functions and the reliability modeling presented in the previous subsections. From
Eq. (19),(24)

is normally distributed
(43)

where
(44)
(45)
Let us deﬁne a vector consisting of observed

as

(25)
If we let

and

be

(46)
Hence,(26)(27)
we can write Eqs. (24) and (25) in vector notation:
(28)
(29)

(47)
We can conclude from Eqs. (47), (40), (41) and (43) that
normally distributed such that

is

(48)

1046

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 23, NO. 6, JUNE 2015

(49)

It can be shown that an iterative procedure consisting of maximizing
with respect to (the maximization step) and(the expectation step) locally maximizes
substituting into
the posterior density. The expectation step computes
andaccording to Eqs. (54) and (55) bysubstituting the current parameter estimate into .
Now, if we partition
into threeblocks andinto nine
blocks such that

As for the prior probability of , we assume that the phrase
control parameter , accent control parameter and state sequence
are independent of each other. Recall that we
is a ﬁrst-order Markov chain.
assumed in IV-A that
We further assume that all other parameters are uniformly distributed. Thus,(56)
where stands for blocks that we can ignore hereafter, the auxiliary function can be rewritten in a more convenient form:

where, , , ,. Overall, the likelihood functionof the Fujisaki model parameters given can be written as

(50)
(51)

V. PARAMETER OPTIMIZATION ALGORITHM(57)

A. Expectation-Maximization (EM) approach
Here we describe an iterative algorithm, which locally maximizes the posterior density ofgiven ,. By regarding the set consisting of
the phrase, accent and base components,,
as the complete data, this problem can be viewed as an incomplete data problem, which can be dealt with using the
Expectation-Maximization (EM) algorithm [22], [23].
The log-likelihood of given the complete data is given as(52)

The update formula for each parameter in the maximization step
can be derived using Eq. (57).
1) State sequence: Leaving only the terms inthat depend on, we have(58). Here the notation
where
is used to denote the -th element of a vector. The state
sequence
maximizing
can be solved
efﬁciently using the Viterbi algorithm as follows. We ﬁrst
set
at
(59)

where denotes equality up to a term independent of . Taking
the conditional expectation of Eq. (52) with respect to givenand, and then adding
to both sides, we obtainan auxiliary function

for all the hidden states . If we consider state
starting state, we shall set
at

to be the(60)

We can compute
via

for

recursively(53)
Because the relationship between the incomplete data and thecomplete data can be written as
where
,
and
are given explicitly as
(54)
(55)

(61)
The most likely transition for each state should be registered at each recursion, so that the most likely state sequence can be traced
at the end of the recursion with
, where. Substituting

KAMEOKA et al.: GENERATIVE MODELING OF VOICE FUNDAMENTAL FREQUENCY CONTOURS

the updated state sequence
into Eq. (20), we ﬁnallyobtain the updatedand .
2) Magnitude of phrase command:ismaximized with respect towhen6) Baseline value
when

1047

:

is maximized with respect to

(69)

(62)
3) Magnitude of accent command
mized with respect towhen:

is maxi-(63)7) Variances of state emission densities
maximized with respect to
and

, :
when

is
(70)
(71)

4) Phrase control parameter : Let us assume a Gaussianprior distribution over such that(64)
Leaving only the terms in
havethat depend on , we

(65)
Now, let

..

.

..

.

..

.

..

.

..

..

.

..

.

..

.

..

.

(66)
.

then from Eqs. (9)–(11),can be written as
(67)

(or
) with respect
The partial derivative of
to is a quartic function, equal up to a constant factor to

(68)
, can
and its roots, namely the stationary points ofbe solved algebraically, from which we can ﬁnd the optimal.
5) Accent control parameter: Let us again assume a Gaussian prior distribution oversuch that. As the derivation follows in exactly
the same manner as above, we shall omit it.

To summarize, we obtain the following iterative algorithm
that guarantees monotonic convergence to a local maximum of:
the posterior densitya) (E-step) Update , , ,,andvia Eqs. (54)and (55).
b) (M-step) Increase
w.r.t.
through the following updates:a) Update by using the Viterbi algorithm.
b) Update via Eqs. (62), (63), (70) and (71).
c) Update (and ) by solving the root of Eq. (68)d) Update
via Eq. (69).
Return to 1) until convergence.
The complexity of this algorithm is, whereisthe number of hidden states.
B. Parameter Inference Under Non-negativity Constraints
It has been shown that phrase and accent commands must
be non-negative in many non-tonal languages such as Japanese,
English, German and Spanish [6]–[9]. In V-A, we treatedandas latent variables (i.e., parameters to be marginalizedout), and did not explicitly take the non-negativity constraints
onand
into consideration. While the method described in V-A can be generally applied even for such languages
as Scandinavian, Portuguese and Chinese in which phrase andaccent commands can be negative [10]–[13], this subsection focuses on parameter estimation under the non-negativity constraint. To impose the non-negativity constraint explicitly, it isconvenient to treat
andas model parameters insteadof latent variables.
Throughout this subsection, let us assume for simplicity thatand are constants. Now, we ﬁrst describe an expanded version of the generative process of :
(72)
(73)
(74)
Note that it can be readily veriﬁed that marginalizing out
andreduces Eqs. (72)–(74) to Eq. (48). Hereafter,
we use to denote the set consisting ofand . Insteadof , we consider treating
as a latent variable. Namely,
we are concerned with maximizing the posterior density. We can obtain the joint probability density
by marginalizing
(75)

1048

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 23, NO. 6, JUNE 2015

with respect to . We notice that the only difference from thejoint probability densitythat we wanted to maximize
in V-A is that is replaced with . For convenience of notation,
we use to denote the set consisting of , and . Hence,.
Here we describe an iterative algorithm that searches for the
maximum a posteriori estimates of
by locally maximizing
given , using the generalized EM algorithm. An auxiliary function for the posterior density of interest can be written
as

where

are auxiliary variables satisfying(78)

We can verify that the equality in this inequality holds whenandare given by
(79)
(80)
We can use this inequality to construct a lower bound function
for
:
(81)
where

(76)
where

(82)
and

Note thatandstand for the discrete-time versions ofand, respectively. It can be shown that an iterativealgorithm that consists of computing
where(via the Forward-Backward
algorithm), increasing
with respect to , and thensubstituting into
locally maximizes the posterior.
Here, we must ensure that increasing
with respect tois performed subject to non-negativity. This can be accomplished by invoking the idea described in [24] as follows.
By using Jensen’s inequality we obtain.
Note that we have used to denoteWe can use this lower bound function to derive an update equation for each element of that guarantees a certain increase of.
As mentioned above, the maximization of this lower boundfunction with respect to can be achieved when it is given
by Eqs. (79) and (80). With
ﬁxed, the maximization ofwith respect to other parameters can be achieved
analytically as follows. First,
is maximized withrespect to when

(83)

(84)

where
(85)
Next,

(77)

is maximized with respect to

when
(86)

KAMEOKA et al.: GENERATIVE MODELING OF VOICE FUNDAMENTAL FREQUENCY CONTOURS

Finally,andis maximized with respect to

,

1049

,

when
(87)
(88)

(89)

(90)
The above equations can be obtained by setting the partial
with respect to,,,derivatives ofandat zero, respectively. Here, it is important
to note that when Eq. (83) (or Eq. (84)) becomes negative,(or) is the maximizer ofsubject tonon-negativity, since is a quadratic (strictly convex) functionof(or). Thus, under the non-negativity constraint,
andshall be written as
the update rules of
(91)
(92)
It can be shown thatis non-decreasing under the updates of and with the above update equations since(93)
To summarize, we obtain the following algorithm that guarantees the convergence to a local maximum of the posterior density of interest and the non-negativity ofand :1) Update using the Forward-Backward algorithm.
w.r.t. through the following updates:
2) Increasea) Update via Eqs. (79) and (80).
b) Update via Eqs. (91) and (92).
c) Update via Eqs. (87)–(90).
d) Update
via Eq. (86).
Return to 1) until convergence.
As with the algorithm in V-A, the complexity of this algorithm is. After convergence, we search for the optimal state sequence from the output sequence
by using the Viterbi algorithm.
VI. EVALUATION OF PARAMETER ESTIMATION ACCURACY
A. Parameter Estimation Using Real Speech Data
To evaluate the parameter estimation accuracy of the algorithms proposed in V-A and V-B, we conducted an experimentusing real speech data, excerpted from the ATR Japanese speech
database B-set [25]. This database consists of 503 phoneticallybalanced sentences. We selected speech samples of one malespeaker (MHT). The ground truth data of the Fujisaki model
parameters had been manually annotated by an expert in the
speech prosody ﬁeld. In these ground truth data, the baseline
values were all set at
Hz. We chose the Fujisaki modelFig. 4. Overview of the experiment in VI-A.

parameter extractor developed by Narusawa [15] as a baseline
method for comparison.
Fig. 4 shows the experimental scheme of the evaluation.
contours were extracted using a method we had previously developed [26], from which the Fujisaki model parameters were
estimated using the present algorithm. V/UV segments wereobtained by simple energy thresholding. The constant parameters were ﬁxed at,ms,
rad s,
rad s,,,
for unvoicedfor voiced regions.
was set at theregions andminimum
value in the voiced regions. The initial values
of were set at the values obtained with Narusawa’s method[15]. The EM algorithm was then run for 20 iterations. The
number of substates in the HMM and the transition probability
were determined according to the manually annotated dataof the ﬁrst 200 sentences. The parameter estimation algorithm
was then tested on the remaining 303 sentences.
We evaluated the accuracy of the parameter estimation basedon the following two criteria: (1) the detection rate of the phrase
and accent commands, and (2) the root mean squared error(RMSE) between an observed log
contour,, and an estiover the voicedmated model,
regions. The aim of this experiment was to conﬁrm whether
the present method is able to achieve accurate model ﬁtting
while ensuring that the estimated parameters are linguisticallyreasonable. The log
RMSE indicates how well the estimatedFujisaki model ﬁts an observed
contour. The detection rateof the phrase and accent commands indicates how linguisticallyreasonable the estimated parameters are. The detection rate of
the phrase and accent commands was calculated by matching
the estimated and ground truth command sequences on acommand-by-command basis using a dynamic programmingalgorithm. If the time difference between the estimated andground truth phrase commands was shorter thanseconds,
the estimated phrase command was considered “matched” andthe local distance was set at zero. Otherwise the local distance
was set at 1. As for the accent commands, we took the averageof the time difference between the onsets of the estimated andground truth accent commands and the time difference between
their offsets. In the same way, when the average time difference
was shorter than seconds, the estimated accent command wasconsidered matched. The magnitudes of the phrase and accentcommands were not taken into account in our evaluation. Thisis because the magnitude estimation was very sensitive to the1050

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 23, NO. 6, JUNE 2015

Fig. 5. Example of command detection (1 of 4). (1) An observed
contour
contours (dotted line) alongin voiced regions (solid line) and the estimatedwith (2) the estimated phrase and accent commands. (a) Narusawa's method[15]. (b) Algorithm proposed in V-A. (c) Algorithm proposed in V-B.

Fig. 6. Example of command detection (2 of 4). (1) An observed
contour
contours (dotted line) alongin voiced regions (solid line) and the estimatedwith (2) the estimated phrase and accent commands. (a) Narusawa's method[15]. (b) Algorithm proposed in V-A. (c) Algorithm proposed in V-B.

baselinevalue, which was set differently in the present
method and in the manual annotation. Let
be the numberof commands in the ground truth command sequences andbe the sum of
for all 303 sentences. We deﬁned thedetection rate as

results show, “P2” yielded the highest model ﬁtting accuracy.
Figs. 5–8 show some examples of observed
contours andcontours obtained with the conventional andthe estimatedpresent methods, from which we can conﬁrm that the present
methods (especially “P2”) were able to ﬁt the Fujisaki model toobserved
contours fairly well.
It should be noted that the detection rate tended to drop signiﬁcantly when
. Considering the fact that averagesyllable durations are typically about 0.2 [s], deviations longer
than 0.1 [s] from the true positions are not negligible. This implies that the proposed methods still have plenty of room for
improvement.

(94)
where
,
and
are the total numbers of insertion, substitution and deletion errors, respectively.
Table I shows the detection rate results for phrase and accentcommands with
. “C,” “P1” and “P2” refer
to Narusawa’s method [15], and the algorithms proposed in V-A
and V-B, respectively. The results showed that the proposed
algorithms “P1” and “P2” were superior to Narusawa’s method,and “P2” was slightly superior to “P1” in terms of detection
rate. The left graph of Fig. 9 shows the log
RMSEs. As the

B. Parameter Estimation Using Synthetic

Contours

To evaluate the pure behavior of the present parameter estimation algorithms, we also conducted a command estimation
contours. The synthetic
experiment using synthetic

KAMEOKA et al.: GENERATIVE MODELING OF VOICE FUNDAMENTAL FREQUENCY CONTOURS

Fig. 7. Example of command detection (3 of 4). (1) An observed
contour
contours (dotted line) alongin voiced regions (solid line) and the estimatedwith (2) the estimated phrase and accent commands. (a) Narusawa's method[15]. (b) Algorithm proposed in V-A. (c) Algorithm proposed in V-B.

contours were artiﬁcially created using the original Fujisaki
model with the abovementioned, manually annotated commandsequences. All other experimental conditions were the same asabove. Fig. 10 provides an overview of this experiment.
Table II shows detection rate results for command sequenceswith different settings. The log
RMSEs are shown in theright graph of Fig. 9. As the results show, the proposed algorithms were again signiﬁcantly superior to the conventional
method in terms of both the detection rate of the command sequences and the model ﬁtting accuracy. However, as regards
the accent command detection rate, the proposed methods were
outperformed by the conventional method. The conventional
method uses the fact that the maxima and minima of the ﬁrst
derivative of thecontour of the Fujisaki model correspondto the onsets and offsets of accent commands with a constantdelay ofif the contributions of the phrase components canbe disregarded. One reason why the conventional method was

1051

contourFig. 8. Example of command detection (4 of 4). (1) An observed
contours (dotted line) alongin voiced regions (solid line) and the estimatedwith (2) the estimated phrase and accent commands. (a) Narusawa's method[15]. (b) Algorithm proposed in V-A. (c) Algorithm proposed in V-B.

Fig. 9. The root mean squared errors (RMSEs) between observed log
tours and estimated models.

con-able to accurately detect accent commands might be that this experiment used synthetic
contours created using the Fujisaki
model as the test data, which agrees particularly well with theabove assumption.

1052

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 23, NO. 6, JUNE 2015

TABLE I
DETECTION RATE (%) OF PHRASE AND ACCENT COMMANDS WITH DIFFERENTTABLE II
DETECTION RATE (%) OF PHRASE AND ACCENT COMMANDS WITH DIFFERENTSETTINGS EVALUATED ON REAL SPEECH DATASETTINGS EVALUATED ON SYNTHETICCONTOURS

ated with the ATR speech samples. We also thank Prof. Shigeki
Sagayama (Meiji University) and Dr. Daisuke Saito (The University of Tokyo) for fruitful discussions.
REFERENCES

Fig. 10. Overview of the experiment described in VI-B.

VII. CONCLUSION
This paper proposed introducing a generative model of voicecontours for estimating prosodic features from raw speech
contour model by transdata. We formulated the present
lating 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. There were two motivations behind this formulation.
One was to derive a general parameter estimation framework for
the Fujisaki model that allows the introduction of powerful algorithms such as the Viterbi algorithm, forward-backward algorithm and EM algorithm. The other was to construct an automatically trainable version of the Fujisaki model that we can naturally incorporate into statistical speech synthesis and conversion
frameworks. We quantitatively evaluated the performance of the
proposed Fujisaki model parameter extractor using real speechdata. Experimental results revealed that our method was superior to a state-of-the-art Fujisaki model parameter extractor. The
application of the present
contour model to prosody generation for text-to-speech synthesis is one of our ongoing projects.
A preliminary study is presented in [27].
ACKNOWLEDGMENT
We thank Prof. Keikichi Hirose (The University of Tokyo)
for kindly providing us with the manually annotated data associ-[1] H. Kameoka, J. Le Roux, and Y. Ohishi, “A statistical model of speech
contours,” in Proc. SCA Workshop Statist. Percept. Audition(SAPA’10), 2010, pp. 43–48.
[2] 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, pp.
175–178.
[3] K. Yoshizato, H. Kameoka, D. Saito, and S. Sagayama, “Hidden
Markov convolutive mixture model for pitch contour analysis of
speech,” in Proc. 13th Annu. Conf. Int. Speech Commun. Association(Interspeech’12), 2012.
[4] H. Kameoka, K. Yoshizato, T. Ishihara, Y. Ohishi, K. Kashino, andcontours,” in Proc.
S. Sagayama, “Generative modeling of speech
14th Annu. Conf. Int. Speech Commun. Assoc. (Interspeech’13), 2013,
pp. 1826–1830.
[5] H. Fujisaki, , O. Fujimura, Ed., “A note on the physiological and physical basis for the phrase and accent components in the voice fundamental frequency contour,” in Vocal Physiology: Voice Production,
Mechanisms and Functions. New York, NY, USA: Raven, 1988.
[6] 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.
[7] H. Fujisaki and S. Ohno, “Analysis and modeling of fundamental frequency contours of English utterances,” in Proc. 4th Eur. Conf. SpeechCommun. Technol. (EUROSPEECH’95), 1995, vol. 2, pp. 985–988.
[8] H. Mixdorff and H. Fujisaki, “Analysis of voice fundamental frequency
contours of German utterances using a quantitative model,” in Proc.
3rd Int. Conf. Spoken Lang. Process. (ICSLP’94), 1994, vol. 4, pp.
2231–2234.
[9] H. Fujisaki, S. Ohno, K. Nakamura, M. Guirao, and J. Gurlekian,“Analysis of accent and intonation in Spanish based on the command-response model,” in Proc. 3rd Int. Conf. Spoken Lang. Process.
(ICSLP’94), 1994, vol. 1, pp. 355–358.
[10] H. Fujisaki, M. Ljungqvist, and H. Murata, “Analysis and modelingof word accent and sentence intonation in Swedish,” in Proc. IEEE
Int. Conf. Acoust., Speech, Signal Process. (ICASSP’93), 1993, pp.
211–214.
[11] H. Fujisaki, S. Narusawa, S. Ohno, and D. Freitas, “Analysis and modcontours of Portuguese utterances based on the commandeling of
response model,” in Proc. 8th Eur. Conf. Speech Commun. Technol.
(EUROSPEECH’03), 2003, vol. 3, pp. 2317–2320.
[12] C. Wang, H. Fujisaki, R. Tomana, and S. Ohno, “Analysis of fundamental frequency contours of standard Chinese in terms of a command-response model and its application to synthesis by rule of intonation,” in Proc. 6th Int. Conf. Spoken Lang. Process. (ICSLP’00),
2000, vol. 3, pp. 326–329.

KAMEOKA et al.: GENERATIVE MODELING OF VOICE FUNDAMENTAL FREQUENCY CONTOURS[13] H. Fujisaki, W. Gu, and K. Hirose, “The command-response model forcontours of Cantonese utterances,” in Proc. 7th
the generation ofInt. Conf. Signal Process. (ICSP’04), 2004, vol. 1, pp. 655–658.
[14] H. Mixdorff, “A novel approach to the fully automatic extraction of
Fujisaki model parameters,” in Proc. IEEE Int. Conf. Acoust., Speech,
Signal Process. (ICASSP’00), 2000, vol. 3, pp. 1281–1284.
[15] S. Narusawa, N. Minematsu, K. Hirose, and H. Fujisaki, “A method forautomatic extraction of model parameters from fundamental frequency
contours of speech,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal
Process. (ICASSP’02), 2002, pp. 509–512.
[16] P. S. Rossi, F. Palmieri, and F. Cutugno, “A method for automatic extraction of Fujisaki-model parameters,” in Proc. Speech Prosody’02,
2002, pp. 615–618.
model for[17] P. S. Rossi, F. Palmieri, and F. Cutugno, “Inversion ofnatural-sounding speech synthesis,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Process. (ICASSP’03), 2003, pp. 520–523.
[18] H. R. Pﬁtzinger, H. Mixdorff, and J. Schwarz, “Comparison of Fustylizers,” in Proc. 10th Annu. Conf.
jisaki-model extractors andInt. Speech Commun. Assoc. (Interspeech’09), 2009, pp. 2455–2458.
[19] J. D. Ferguson, “Variable duration models for speech,” in Proc. Symp.
Applicat. Hidden Markov Models to Text and Speech, 1980, pp.
143–179.
[20] M. Russell and R. Moore, “Explicit modeling of state occupancy inhidden Markov models for automatic speech recognition,” in Proc.
IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP’85), 1985,
pp. 5–8.
[21] S. Levinson, “Continuously variable duration hidden Markov models
for speech analysis,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal
Process. (ICASSP’86), 1986, pp. 1241–1244.
[22] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihoodfrom incomplete data via the EM algorithm,” J. R. Statist. Soc. Ser. B,
vol. 39, no. 1, pp. 1–38, 1977.
[23] M. Feder and E. Weinstein, “Parameter estimation of superimposed
signals using the EM algorithm,” IEEE Trans. Acoust., Speech, Signal
Process., vol. 36, no. 4, pp. 477–489, Apr. 1988.
[24] H. Kameoka, T. Nakatani, and T. Yoshioka, “Robust speech dereverberation based on non-negativity and sparse nature of speech spectrograms,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process.
(ICASSP’09), 2009, pp. 45–48.
[25] A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri, H. Kuwabara, andK. Shikano, “ATR Japanese speech database as a tool of speech recognition and synthesis,” Speech Commun., vol. 27, pp. 187–207, 1999.
[26] H. Kameoka, “Statistical approach to multipitch analysis,” Ph.D. dissertation, The Univ. of Tokyo, Tokyo, Japan, 2007.
[27] K. Kadowaki, T. Ishihara, N. Hojo, and H. Kameoka, “Speech prosody
generation for text-to-speech synthesis based on generative model ofcontours,” in Proc. 15th Annu. Conf. Int. Speech Commun. Assoc.
(Interspeech’14), 2014, pp. 2322–2326.

Hirokazu Kameoka received B.E., M.S. and Ph.D.
degrees all from the University of Tokyo, Japan, in
2002, 2004 and 2007, respectively. He is currently
a Distinguished Researcher at NTT Communication Science Laboratories, Nippon Telegraph andTelephone Corporation, and is also an Adjunct
Associate Professor at the University of Tokyo.
His research interests include audio, speech, andmusic signal processing and machine learning. He
has served as a publicity chair of ISMIR 2009 anda program committee member of ISMIR 2014 and
ISMIR 2015. Since 2015, he has been an associate editor of the IEEE/ACM
TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING. He is amember of IEEE, the Institute of Electronics, Information and Communication
Engineers (IEICE), the Information Processing Society of Japan (IPSJ) andthe Acoustical Society of Japan (ASJ). He received 13 awards over the past 10
years, including the IEEE Signal Processing Society 2008 SPS Young Author
Best Paper Award. He is the author or co-author of about 90 articles in journal
papers and peer-reviewed conference proceedings.

1053

Kota Yoshizato received the B.E. and M.S. degrees
from the University of Tokyo, Japan, in 2011 and2013, respectively. Since then, he has been a Software Engineer at DeNA Corp.

Tatsuma Ishihara received the B.E. and M.S. degrees from the University of Tokyo, Japan, in 2012
and 2014, respectively. Since then, he has been a Research Engineer with the Toshiba Corporate Research& Development Center.

Kento Kadowaki received the B.E. degree from
Osaka University, Japan, in 2013 and the M.S.
degree from the University of Tokyo, Japan, in 2015.
He is currently with Nomura Research Institute Co.,
Ltd., Japan.

Yasunori Ohishi received his Ph.D. degree in information science in 2009 from Nagoya University,
Japan. In 2009, he joined NTT Communication Science Laboratories. His research interests generally
concern statistical acoustic signal processing, music
information retrieval, and unsupervised machinelearning with audio applications. In 2009, he received the 36th Awaya Kiyoshi Science Promotion
Award from Acoustical Society of Japan (ASJ).

Kunio Kashino received B.E., M.S., and Ph.D.
degrees all from the University of Tokyo, Japan,
in 1990, 1992, and 1995, respectively. He is currently a Senior Distinguished Researcher at Nippon
Telegraph and Telephone Corporation, Executive
Manager of Media Information Laboratory, NTT
Communication Science Laboratories, and a Visiting
Professor at National Institute of Informatics. His
research interests include media information processing and audio information processing. He has
been involved in the research projects on computational auditory scene analysis since 1990 and robust media search since 1996.
He is a senior member of the IEEE and the Institute of Electronics, Information
and Communication Engineers (IEICE), and a member of the Information
Processing Society of Japan (IPSJ), the Acoustical Society of Japan (ASJ),
and the Japanese Society for Artiﬁcial Intelligence (JSAI). He received morethan 10 awards over the past 20 years, including the IEEE TRANSACTIONS ON
MULTIMEDIA Prize Paper Award.