FAST ALGORITHM FOR STATISTICAL PHRASE/ACCENT COMMAND ESTIMATIONBASED ON GENERATIVE MODEL INCORPORATING SPECTRAL FEATURES
Ryotaro Sato1 , Hirokazu Kameoka2 , Kunio Kashino2
1

Graduate School of Information Science and Technology, The University of Tokyo
2
NTT Communication Science Laboratories, NTT Corporation
ABSTRACTAn important challenge in speech processing involves extracting
non-linguistic information from a fundamental frequency (F0 ) contour of speech. We propose a fast algorithm for estimating the model
parameters of the Fujisaki model, namely, the timings and magnitudes of the phrase and accent commands. Although a powerfulparameter estimation framework based on a stochastic counterpartof the Fujisaki model has recently been proposed, it still had room
for improvement in terms of both computational efficiency andparameter estimation accuracy. This paper describes our two contributions. First, we propose a hard expectation-maximization (EM)algorithm for parameter inference where the E step of the conventional EM algorithm is replaced with a point estimation procedure
to accelerate the estimation process. Second, to improve the parameter estimation accuracy, we add a generative process of a spectralfeature sequence to the generative model. This makes it possible
to use linguistic or phonological information as an additional clueto estimate the timings of the accent commands. The experiments
confirmed that the present algorithm was approximately 16 times
faster and estimated parameters about 3% more accurately than theconventional algorithm.
Index Terms— voice fundamental frequency contour, Fujisaki
model, prosodic information processing, EM algorithm
1. INTRODUCTION
Various speech recognition and synthesis systems have become
widespread in our daily lives, and the importance of non-linguistic
information in spoken communication has been highlighted in recent
years [1]. Voice fundamental frequency (F0 ) contours constitute oneof the most important physical entities correlated with socially essential information, such as the speaker’s identity, intention, attitude,and mood.
An F0 contour typically consists of relatively long-term frequency variations over the duration of a phrase and short-term
variations such as those found in accented syllables. The Fujisaki
model [2] describes a log F0 contour as the sum of these two components, which are called phrase and accent components, respectively. These components are modeled as the responses of criticallydamped second-order linear systems to impulse and step inputs,
called commands. The basic underlying parameters consist of thetiming and magnitudes of the impulse inputs (phrase commands)
and the onset/offset times and magnitudes of the step inputs (accentcommands).
The Fujisaki model allows us to approximate actual F0 contours
fairly well using a relatively small number of parameters when theparameters are appropriately chosen and has been applied to various
languages in multiple application areas. Although the parameters are

978-1-5090-4117-6/17/$31.00 ©2017 IEEE

5595

Fig. 1. HMM generating phrase/accent command functions.

generally difficult to estimate automatically from a raw F0 contour[3], a powerful parameter estimation framework based on a stochastic counterpart of the Fujisaki model has recently been proposed [4].
This method has been reported to yield reasonably good results by
taking advantages of powerful statistical inference and training techniques for the estimation. However, when we consider the various
applications that must be run in real time or faster under realistic
circumstances, its computational efficiency and accuracy need to befurther improved.
On the basis of the generative model described above, we propose an efficient algorithm for parameter inference where the E stepof the conventional EM algorithm is replaced with a point estimation procedure. We also propose using spectral information as anadditional clue to accurately estimate the parameters, assuming that
accents appearing in F0 contours are related to linguistic or phonological information. In the rest of this paper, Section 2 describes
the addressed problem, and Section 3 introduces the proposed algorithm. Section 4 reports the experimental results, and Section 5
concludes this paper.

2. GENERATIVE MODEL APPROACH TO
F0 PARAMETER ESTIMATION
2.1. Generative Model of F0 Command Functions
Several researchers have addressed the automatic parameter estimation of the Fujisaki model [3]. Our method is based on the generativemodel approach introduced by Kameoka et al. [4]. We briefly review
its framework.
The Fujisaki model assumes that an F0 contour y[t] on a loga-

ICASSP 2017

Fig. 2. State splitting to flexibly assign probabilities to state durations.
Fig. 3. Graphical representation of proposed model.
rithmic scale is given as the sum of three components:
y[k] = xp [k] + xa [k] + µb ,(1)

where xp [k] and xa [k] respectively are a phrase component and an
accent component at time frame k, and µb is a constant value. The
phrase and accent components are assumed to be the outputs of different second-order critically damped filters, Gp [k] and Ga [k], excited with pulse sequence up [k] (phrase commands) and rectangularpulse sequence ua [k] (accent commands), respectively:
xp [k] = Gp [k] ∗ up [k],(2)xa [k] = Ga [k] ∗ ua [k],(3)

• E step: Update P (sk = q|y, θ, µb , o) = P (sk = q|θ, o) for
each frame k and each state q using the Forward-Backward
algorithm.
• M step: Update o, θ, µb using the auxiliary function.
After convergence, MAP estimation for the state sequence s is
performed using the Viterbi algorithm to obtain the estimated command sequence oest = {[µp [k], µa [k]]}K−1
k=0 .
3. PROPOSED METHOD

where ∗ is convolution over time. The key idea of Kameoka’s approach [4] is that a phrase/accent command pair sequence is modeled
as an output sequence of the path-restricted hidden Markov model(HMM) with Gaussian state emission densities shown in Fig. 1.
Here the states p0 , an (n = 0, · · · , N − 1), and ri (i = 0, 1),
respectively, are states where only the phrase command is being
activated, i.e., µp [k] = E[up [k]] = C (p) [k], µa [k] = E[ua [k]] =0, only the accent command is being activated, i.e., µp [k] =(a)0, µa [k] = Cn , and neither command is being activated, i.e.,µp [k] = 0, µa [k] = 0. The path constraint shown in Fig. 1
restricts µp [k] = E[up [k]] to consisting of isolated deltas andµp [k] = E[µa [k]] to consisting of rectangular pulses. Hence,
estimating the state transition of the HMM directly amounts to estimating the Fujisaki-model parameters. Furthermore, we can directly
assign probabilities to the durations of self-transition by splitting
each state (except for p0 ) into a certain number of sub states, as
shown in Fig. 2. Given a state sequence s = (s0 , . . . , sK−1 ), thisHMM emits o[k] = [up [k], ua [k]]T according to the Gaussian state
emission density:o[k] ∼ N (o; ρ[k], Υ),(4)T

where ρ[k] = [µp [k], µa [k]]and Υ = diag[σp2 , σa2 ].

2.2. Conventional Parameter Estimation Algorithm
By using the conditional density P (y|o, µb ) reflecting the constraints (1)-(3) of the Fujisaki model, and the HMM likelihood
P (o|s, θ), the joint probability density function is written as
P (s, θ, µb , o, y) = P (y|o, µb )P (o|s, θ)P (s)P (µb )P (θ), (5){}(a) N −122where θ = {C (p) [k]}K−1
contains all the
k=0 , {Cn }n=0 , σp , σa
HMM parameters except for the transition probabilities, and s =K−1{sk }K−1
k=0 represents the state sequence. y = {y[k]}k=0 is the observed F0 contour. The statistical estimation is then summarized tomaximize P (o, θ, µb |y) for a given y. For parameter inference, we
previously proposed an algorithm based on an auxiliary function andthe EM algorithm [5] by treating s as a latent variable to be marginalized out and θ and µb as the model parameters to be estimated. The
procedure is as follows:

5596

We propose two improvements to the conventional method describedabove: the use of the hard EM algorithm and the introduction ofspectral features.
3.1. Modified EM Algorithm
In the method described in Sec. 2.2, our analysis revealed that
more than 90% of the computation time was spent by the ForwardBackward algorithm when computing the posterior state probabilities at each frame in the E step. Each output distribution of the
HMM is a Gaussian distribution, and so a logarithm and an exponential must be computed for every addition at each frame, and thisis computationally expensive. To sidestep these computations, we
propose using the joint probability density function P (s, θ, µb , o|y)instead of the marginal distribution P (θ, µb , o|y) as the objectivefunction. This results in a parameter estimation algorithm in which
the Forward-Backward procedure in the conventional algorithm
is replaced with a state decoding procedure using the Viterbi algorithm. Since the Viterbi algorithm is generally faster than the
Forward-Backward algorithm, we expect the entire computation
time to be reduced. This is called a hard EM algorithm.
The variables involved in this method are summarized at theright half of Fig. 3. The estimation process using the modified EM
algorithm can be written as follows:• E step (hard EM): Update s by maximizing P (o|s, θ)P (s)using the Viterbi algorithm.
• M step:
Update o, θ, µb (as with the conventional method).
We conjecture that this method can also improve the estimation accuracy. This can be explained in terms of the objectivefunctions considered in the conventional and proposed algorithms.
The conventional method can be interpreted as the optimization ofP (θ, µb , o|y), with s being integrated out. However, the parameters
we need to infer are µp [k] and µa [k] for each frame k, which are
uniquely determined by (θ, s). Therefore, the objective functionthat we should actually maximize is P (θ, µb , o, s|y) rather thanthat mentioned above. The conventional method [4] performs a
MAP estimation for s only at the final step, so the finally inferredparameters (µp [k], µa [k]) are not guaranteed to correspond to the
local maximum of P (θ, µb , o, s|y).

In contrast, our method maximizes P (θ, µb , o, s|y) directly.
This objective function can be factorized as Eq. (5) (ignoring a
constant factor) and so the modified E step can be interpreted as theconditional maximization of this function with respect to s while
keeping the other parameters fixed.
3.2. Introduction of Spectral Features
Several researchers have addressed the use of linguistic information
to automatically estimate F0 contour parameters [6, 7, 8]. The basic idea behind these approaches is the assumption that accent command functions are somewhat related to linguistic information. Morespecifically, the onsets of the accent commands tend to occur whenphoneme transitions occur. Since phoneme transitions are usuallyaccompanied by rapid changes in spectral feature values, we can expect the accuracy of accent command estimation to be improved if
we assume that the temporal variations of the spectral feature values
are likely to be large when an accent command is activated.
Fig. 3 shows the relationships between the variables in the proposed probabilistic model comprising both spectral variations andF0 commands.
Phoneme transition probability : z = {z[k]}K−1
k=0(0 ≤ z[k] ≤ 1)

(6)

Output sequence : v = {v[k]}K−1
k=0

(7)

{
P (z[k]|sk−1 , sk ) =f (z[k], φ)1(sk−1 = r1 , sk = an )(otherwise)

P (v[k]|z[k], Σ0 , Σ1 )
= (1 − z[k])N (v[k]|0, Σ0 ) + z[k]N (v[k]|0, Σ1 )

(8)

(9)

z[k] represents the probability that the kth frame corresponds
to the phoneme boundary. The prior of z[k] is given by the function f and the parameter φ. The distribution of the spectral feature vector v[k] is expressed by a two-component Gaussian mixture
distribution where z[k] is the weight. One component of this mixture is intended to express the distribution of the spectral features atphoneme boundaries and the other is intended to express that within
other segments. The norm of a ∆-spectral feature vector tends to belarge at a phoneme boundary, and our proposed probabilistic model
is designed to take this characteristic into account.
The F0 estimation process based on the EM algorithm is summarized as follows:• E step (hard EM): Update s by maximizing
P (z|s)P (o|s, θ)P (s) using the Viterbi algorithm.
• M step:
1. Update the variables of the prosody model s, θ, µb (thesame as the entire M step of the conventional method).
2. For each k ∈ [1, K − 1] that satisfies sk−1 = r1 andsk = ai , update z[k] to 1.
3. For each k ∈ [0, K − 1] where the kth frame is judged
to be unvoiced, update z[k] to 0.
4. Update Σ0 , Σ1 using the update function of the Gaussian Mixture Model P (v|z, Σ0 , Σ1 )(Here, we assume that Σ0 and Σ1 are diagonalized).
5. Update z by maximizing P (z|Σ0 , Σ1 , v).

5597

4. EXPERIMENTS
To evaluate the speed and accuracy of the methods, experiments
were conducted using real speech data, excerpted from the ATR
Japanese sentence database B-set [9]. We used 503 sentences spoken by one male speaker (MHT). For the ground truth data, we used
manually annotated Fujisaki model command function parameters.
F0 contours were estimated by using Nakatani’s method [10],and the baseline values µb were set for each sentence by using theminimum value of F0 in the voiced segments. The initial values of
the command functions were extracted by using Narusawa’s method[3].
We estimated the command functions using three methods: (1)the conventional method (C) [4], (2) the method introduced in Section 3.1 (P1), and (3) the method described in Section 3.2 (P2). For
the distribution of z in P2, we used f (z[k], a) = C exp(az[k]),and we empirically chose a = 10.0. For the spectral features,
we used ∆MFCC (to 12th order), ∆LPC (to 20th order), ∆power,∆LPC+∆power.
Throughout the experiments, we fixed α = 3.0 rad/s, β = 20.0
rad/s, N = 10, t0 = 8 ms, vp2 = 0.032 , va2 = 0.032 , vn2 = 1002 forunvoiced segments and vn2 = 0.032 for voiced segments. We divided
the 503 manually annotated sentences into a training set and a test
set; that is, the transition probabilities of the HMM were obtained in
the training step using the test set comprising the first 200 sentencesfrom the corpus, and then the evaluation test was performed using
the remaining 303 sentences.
We used two criteria as measures of estimation accuracy: (1) theroot mean squared error (RMSE) between the observed logF0 andthe synthesized value Gp [k] ∗ µp [k] + Ga [k] ∗ µa [k] + µb using theestimated µp and µa and the impulse responses Gp [k], Ga [k], overthe voiced segments, and (2) the deletion and insertion rates of thephrase and accent commands. The second criterion must be usedbecause the objective of the estimation is not only to minimize the
RMSE but also to obtain a compact description of the F0 contours.
The deletion/insertion rates are defined as follows. First, we matchthe estimated and ground truth command sequences on a commandby-command basis with a dynamic programming algorithm. By using a predefined time difference tolerance S, the estimated commands that found a match were judged “matched.” With the accentcommands, we compare the mean value of the onset time difference
and the offset time difference with S. Let N , Nest and Nmatchbe the number of commands in the ground truth data, the number of commands in the estimated result, and the number of commands judged to be matched, respectively. By using these values,
the deletion/insertion rates are defined as pdel = (N − Nmatch )/Nand pins = (Nest − Nmatch )/N . The numerators in this definition mean the number of deleted and inserted commands in the estimated results compared with the ground truth data. The relationshipbetween these rates and the detection rate pdet in [4] is written as
pdet = 1 − max(pdel , pins ). Note that the magnitudes of the commands are not introduced in the second criterion because the aim isto evaluate appropriateness of estimation in terms of the number ofcorrectly estimated commands. We performed the evaluation using
S = 0.1 sec throughout our experiments.
Fig. 4 shows the RMSE and insertion/deletion rates of command
functions estimated with the conventional and proposed methods described above. There were 20 iterations in the EM algorithm in thisexperiment.
With all of the proposed methods, the RMSE values are improved compared with the conventional method. The total deletion/insertion rate, shown in the right panel, stands for the mean

Table 1. Computation time for command function estimation.
MethodC [4]P1 (hard EM)
P2 (∆MFCC)

E step [s/iteration]
1.610
0.023
0.049

M step [s/iteration]
0.039
0.042
0.050

CPU: Core i7-6700K 4.0GHz, RAM:32GB
Windows 7 SP1, MATLAB R2016a
Length of sound data: 3.62 seconds

Fig. 4. RMSE and deletion/insertion rate.

Fig. 6. Example of processing results. Purple dotted line overlapping estimated value of command functions obtained with method
P2 represents z at each frame.
Fig. 5. Number of iterations vs. RMSE and insertion/deletion rates.
5. CONCLUSIONS
value of phrase/accent deletion and insertion rate weighted by the
numbers of commands in the ground truth data. The total rates
are also improved, and the technique using ∆MFCC achieved thebest rate among the methods tested here, reducing the total insertion/deletion error by 3 %.
The computation times are listed in Table 1, showing that the
proposed methods greatly increased the speed. For example, thecomputation time in E step is 70 times shorter in the hard EM method(P1) than in the conventional method (C), and the total computation
time is 16 times shorter for the proposed ∆MFCC method (P2) thanfor C. This means that the real time factor of the estimation process
is less than 1 with the proposed method (P2), if there are fewer than
36 iterations.
Fig. 5 shows the dependence of the RMSE and the insertion/deletion rates for phrase/accent commands on the number of
iterations. The RMSE decreases monotonically for all algorithms.
Nevertheless, the insertion/deletion rate obtains a local minimum
at around eight iterations. This phenomenon is presumed to occur
because our experiments were not conducted to fit the F0 contours
generated by the ground truth command functions but to fit theobserved F0 contours.
Fig. 6 shows an example of the observed F0 contour, the groundtruth of the commands, the estimated commands (P1 and P2 with
∆MFCC), and the ∆MFCC values. The dotted purple curve in thegraph second from the bottom visualizes the z estimated in P2. Thisshows that the phonological information successfully prevents theunnecessary switching of accent commands right after t = 3 sec.

5598

We described a way of improving the Fujisaki model parameter
estimation algorithm proposed in [4] in terms of computational
efficiency and parameter estimation accuracy. To accelerate thealgorithm, we proposed replacing the E step in the conventional algorithm with a state decoding procedure using the Viterbi algorithm.
To improve the estimation accuracy, we proposed adding a ∆spectral feature generative process to the original generative model,
which made it possible to estimate the F0 commands from F0 andspectral features. Experiments using real speech data showed thatthe proposed method achieved an approximately 70 times shortercomputational time in Estep and an approximately 16 times shorter
whole calculation time than the conventional method. The spectralfeatures were also shown to be effective. The detection rate was
the best when using ∆MFCC as the spectral features, and the number of wrong F0 command insertions or deletions can be reducedby 3 % compared with the case where no spectral features wereused. As future work, we plan a further study of the use of spectral or phonological information to improve the robustness of theestimation.
6. ACKNOWLEDGMENTS
We thank Prof. Keikichi Hirose for providing us with the manuallyannotated data associated with the ATR speech samples. We also
thank Dr. Tomohiro Nakatani (NTT), who kindly provided us withthe source code of the F0 extractor used in the experiment.

7. REFERENCES[1] K. Hirose, J. Tao (eds.): “Speech prosody in speech synthesis: Modeling and generation of prosody for high quality andflexible speech synthesis,” Springer (2015).
[2] H. Fujisaki, O. Fujimura (eds.): “A note on the physiologicaland physical basis for the phrase and accent components in thevoice fundamental frequency contour,” in Vocal Physiology:
Voice Production, Machanisms and Functions, Raven (1988).
[3] S. Narusawa, N. Minematsu, K. Hirose, H. Fujisaki: “A
method for automatic extraction of model parameters from fundamental frequency contours of speech,” in Proc. IEEE Int.
Conf. Acoust. Speech, Signal Proc. (ICASSP), pp. 509–512(2002).
[4] H. Kameoka, K. Yoshizato, T. Ishihara, K. Kadowaki, Y.
Ohishi, K. Kashino: “Generative modeling of voice fundamental frequency contours,” IEEE/ACM Trans. on Audio, Speechand Language Proc., vol.23, no.6, pp.1042–1053 (2015).
[5] A. P. Dempster, N. M. Laird, D. B. Rubin: “Maximum likelihood from incomplete data via the EM algorithm,” J. R. Statist.
soc. Ser. B, vol.39, no.1, pp.1–38 (1977).
[6] K. Hirose, Y. Furuyama, S. Narusawa, N. Minematsu, H. Fujisaki: “Use of linguistic information for automatic extractionof F0 contour generation process model parameters,” in Proc.
Eurospeech pp.141–144 (2003).
[7] H. M. Torres, J. A. Gurlekian, H. Mixdorff, H. Pfitzinger: “Linguistically motivated parameter estimation methods for a superpositional intonation model,” EURASIP Journal on Audio,
Speech, and Music Processing (2014).
[8] K. Hirose, Y. Furuyama, N. Minematsu: “Corpus-based extraction of F0 contour generation process model parameters,”
in Proc. Interspeech pp. 3257–3260 (2005).
[9] A. Kurematsu, K. Takeda, Y. Sagisaka, S. katagiri, H.
Kuwabara, K. Shikano: “ATR Japanese speech database asa tool of speech recognition and synthesis,” Speech Comm.,
vol.9, no.4, pp.357–363 (1990).
[10] T. Nakatani, S. Amano, T. Irino, K. Ishizuka, T. Kondo: “Amethod for fundamental frequency estimation and voicing decision: Application to infant utterances recorded in real acoustical environments,” Speech Comm., vol.50, no.3, pp. 203-?214(2008).

5599