INTERSPEECH 2010

Statistical Modeling of F0 Dynamics in Singing Voices
Based on Gaussian Processes with Multiple Oscillation Bases
Yasunori Ohishi, Hirokazu Kameoka, Daichi Mochihashi, Hidehisa Nagano, Kunio Kashino
NTT Communication Science Laboratories, NTT Corporation{ohishi,kameoka,nagano,kunio}@cs.brl.ntt.co.jp, daichi@cslab.kecl.ntt.co.jp

10.21437/Interspeech.2010-254

F0 [cent]

Abstract
We present a novel statistical model for dynamics of various
singing behaviors, such as vibrato and overshoot, in a fundamental frequency (F0) contour. These dynamics are the important cues for perceiving individuality of a singer, and can be auseful measure for various applications, such as singing skill
evaluation and singing voice synthesis. While most previous
studies have modeled the dynamics using a second-order linear
system, the automatic and accurate estimation of model parameters has yet to be accomplished. In this paper, we first develop a complete stochastic representation of the second-ordersystem with Gaussian processes from parametric discretization,and propose a complete, efficient scheme for parameter estimation using the Expectation-Maximization (EM) algorithm. Experimental results show that the proposed method can decompose an F0 contour into a musical component and a dynamics
component. Finally, we discuss estimating singing styles fromthe model parameters for each singer.
Index Terms: Singing voices, Fundamental frequency (F0),Gaussian Processes, EM algorithm, Singing voice synthesis

46
System output:
Step signal:
Impulse response:

[sec]

8

2

F0 sequence:

Residue:

Second-orderlinear system

Gaussian
white noise
0

0

Time

0

Time

0

Time

Time

Figure 1: Decomposition of F0 contour into three components
based on second-order linear system
On the other hand, we have previously proposed a stochastic framework that learns these parameters and the step-wisesignal simultaneously from an observed F0 contour [13]. As
seen in Fig. 1, we divided the F0 contour into segments anddecomposed the F0 contour for each segment into three components, i.e., the step signal f (t), the impulse response h(t) of thesecond-order system, and the residue (t). f (t) indicates therelative pitch that the singer attempts to sing. h(t) represents
the temporal attack of the musical note and overshoot. System
output y(t) denotes the convolution of f (t) and h(t). (t) is
the residual component including vibrato. Therefore, we approximated the second-order differential equation given by the
inverse Laplace transform of Eq. (1) by the difference equation:

1. Introduction
The goal of this work is to characterize both musical-note information and the dynamics of various singing behaviors, suchas vibrato and overshoot, in a sung melodic contour, i.e., an F0
sequence. The importance of dynamics for perceiving individuality of a singer is reported in [1, 2] based on psycho-acoustic
experiments, and this means that the dynamics can be a usefulmeasure for identifying of singing styles [3], singing skill evaluation [4], and the synthesis of more natural and various singing
voices [5, 6, 7, 8]. On the other hand, musical-note information
is useful for various applications such as Query-By-Hummingand the automatic clustering of songs [9, 10].
Most previous studies have used a deterministic, secondorder linear system to represent the F0 dynamics of singing
voices [1, 11]. The transfer function of the system is described
as
Ω2
,
(1)
H(s) = 2
s + 2ζΩs + Ω2
where ζ and Ω denote the damping ratio and the damped natural frequency, respectively. In the F0 control model for singing
voices, the dynamics are represented not only by Critical damping (ζ = 1) used in the Fujisaki model [12] that controls the F0
contours of speaking voices, but also by Over-damping (ζ > 1),
Under-damping (0<ζ<1), and Steady oscillation (ζ=0). Thisis because the F0 fluctuations in singing voices are larger andmore rapid than those of speaking voices. In [1], using theF0 contour generated by the convolution of a step-wise signal
which corresponds to the musical-note sequence and the impulse response of Eq. (1), natural and expressive singing voices
were synthesized. However, ζ and Ω were controlled manuallybased on psycho-acoustic experiments.

Copyright © 2010 ISCA

6000
5800
5600
5400
5200
5000
0

(a := 1/Ω2 , b := 2ζ/Ω),(aA + bB + C)y = f(2)

where f = [f1 , f2 , . . . , fN ] and y = [y1 , y2 , . . . , yN ] denote discrete signals sampled with a sampling period Δ. N is
the signal length, and aA + bB + C represents the inverse impulse response. We define A, B, and C as follows:
T

2 1

3

Δ

T

2

1
2Δ

1
6 0
6 −2
7
0
Δ
6 1
6Δ
7
−2
1
6
6 1
7
Δ
Δ
A := 6 Δ
7 , B := 6− 2Δ
6
6
7
..
..
..
4
4
5
.
.
.
−2
1
1
0
0
Δ
Δ
Δ
3
2
0
0
7
61 0
7
60 1
0
7.
C := 6
7
6
.
.
.
5
4
..
..
..

0

0

1

3
1
2Δ

0
..
.

7
7
7
7,
7
5

0
1
2Δ

..
.
1
− 2Δ

..

.
0

1
2Δ

(3)

0

Then, we estimated the model parameters iteratively based on
the maximum-likelihood approach, so that the residual error between both sides of Eq. (2) were minimized. However, the estimation accuracy was poor in parts of the transition between
musical notes. Since various dynamics, such as the attack andvibrato, are mixed in each segment, we consider that over-fitting
occurs with the conventional method which estimates the model
parameters of the second-order linear system directly.
In this paper, we propose a complete stochastic representation of the second-order linear system with Gaussian processes

2598

26 - 30 September 2010, Makuhari, Chiba, Japan

(GPs) [14]. Therefore, we parametrically approximate the impulse response of the system using multiple oscillation basesand explicitly introduce the residual component › as a randomvariable. Then, we develop a efficient scheme for parameter
estimation using the EM algorithm. Our experimental results
show that the proposed method can appropriately decompose
the F0 contour into the target musical note and the dynamics.
Furthermore, we also discuss estimating singing styles from themodel parameters for each singer.

It should be noted that this model provides us with a particular example of GPs. A key point about GPs is that the joint
distribution over N elements of y is specified completely by
the second-order statistics, namely the mean and covariance. In
vanilla GPs, since the mean is mainly zero, the covariance matrix, which is usually denoted by the Gram matrix K, is represented by combining multiple kernel matrices linearly, referred
to as the Multiple Kernel Learning (MKL) [15]. This technique
can exploit the temporal correlations between observations instead of treating the observations as i.i.d. The GPs with the
MKL have received increased attention in the machine-learning
community. However, in our case, since the mean includes
the variable Ψ−1 , we represent Ψ by the linear combinationof multiple oscillation bases. In this respect, we refer to thisrepresentation as GPs with the multiple basis learning.

2. Impulse response approximation withmultiple oscillation bases
We introduce a discrete-time representation of Eq. (1) to replace
the difference approximation Eq. (2). The impulse response of
Eq. (1) has four cases:8√√´−ζΩt `
ζ 2 −1Ωt− ζ 2 −1ΩtΩe>√e, (ζ > 1)
−
e
>
>
2 −1
>
2
ζ
> −ζΩt `
p
´
<
Ωe
2
√
sin( 1 − ζ Ωt ,(0 < ζ < 1)h(t) =
1−ζ 2
>
2
−Ωt
>
>Ω te
,
(ζ = 1)>>:
Ω sin(Ωt),(ζ = 0)3.2. Likelihood function and prior probability
We introduce the residual component › (Gaussian white noise,› ∼ N (0, βI N )) and assume that the observed F0 sequenceo = [o1 , o2 , . . . , oN ]T is given by output y with › so thato = y + ›, where β is a hyperparameter representing the variance of the noise. Since we suppose that y and › are mutuallyindependent, from the definition of GPs, the likelihood functionof Θ := {Ψ, u, β} is given as follows:

Discretizing these impulse responses at Δ, the input-output
relation is described by y = Φf . For example, when ζ = 1, Φ
is described as a lower triangular matrix
2Ω2 Δe−ΩΔ
6 2Ω2 Δe−2ΩΔ
6
Φ=6
..
4
.
2
N Ω Δe−N ΩΔ

30Ω2 Δe−ΩΔ
..
.
...

7
7
7.
5

..

.
2Ω2 Δe−2ΩΔ

n 1
o
1
T −1
exp
−
Σ
(o
−
—)
,
(o
−
—)
P (o|Θ) =2(2π)N/2 |Σ|1/2”“(6)— = Ψ−1 u, Σ = αΨ−1 (Ψ−1 )T + βI N .
We assume that the prior distributions for Ψ, u, and β are
independent, which yields P (Θ) = P (Ψ)P (u)P (β), and that
P (u) and P (β) are uniform distributions. For P (Ψ), we also
assume the independence of parameters, which yields

Ω2 Δe−ΩΔ

However, since h(t) has multiple cases, we compose matrix Φ
as follows:Φ−1  w1 Υ(1) + w2 Υ(2) + . . . + wI Υ(I) ,j(4)
P (Ψ) =

where we determine the values of ζ and Ω manually and preliminarily calculate I oscillation bases {Φ(1) , Φ(2) , . . . , Φ(I) }.
For convenience, we define the impulse response Υ(i) :=(Φ(i) )−1 of the inverse filter and approximate Φ−1 by combining these matrices linearly, where we assume that weight parameters w := {w1 , w2 , . . . , wI } are sparse. This means thatΦ−1 are expressed by only a few oscillation bases. Therefore,
the input-output relation of the system is obtained by(w1 Υ(1) + w2 Υ(2) + . . . + wI Υ(I) )y = f .

P (a)P (b),P (w),(Diﬀerence approximation),(Impulse response approximation)and that P (a) and P (b) are uniform distributions. P (w) corresponds to the sparsity cost described in Section 2, for which anatural choice is a generalized Gaussian prior
P (w) =

(5)

I
Y

p
p
λp
exp−λ |wi |2Γ(1/p)i=1(7)Note that the transfer function Eq. (1) is converted into the form
Ψy = f as seen from Eqs. (2) and (5).

where p and λ are the parameters that determine the shape of the
distribution. When 0 < p < 2, P (w) becomes super-Gaussianand promotes sparsity if the norm of w is bounded.

3. Statistical modeling of second-orderlinear system based on Gaussian Processes

4. Parameter estimation algorithm basedon Expectation-Maximization algorithm

We statistically model the input-output relations of the systemrepresented by Eqs. (2) and (5) based on GPs.

Given observed F0 sequence o, we want to determine the estimate of Θ that maximizes the posterior density P (Θ|o) ∝P (o|Θ)P (Θ). However, it is difficult to obtain an optimum
solution for the maximum a posterior (MAP) estimate of Θ analytically. This is because (1) o is given by y with ›, (2) the
objective function is nonlinear with respect to w. To cope withthese problems, (1) we partition o into y and › using the EM
algorithm, (2) Applying the auxiliary function method [16] to
the M-step, we design an auxiliary function of the Q-function.

3.1. Modeling of input step signal
We assume input f is a step signal (see Fig. 1). We statisticallymodel f as a random variable generated from the multivariate Gaussian distribution N (u, αI N ), where we define u =[u1 , . . . , uN ]T = u[1, 1, . . . , 1]T = u1. Scalars u and α are
the relative pitch and a hyperparameter representing the variance of the distribution, respectively. α is set experimentally.
I N denotes the N × N identity matrix. Since y is a linear combination of Gaussian distributed variables
given by the elements´`
of f , y is itself Gaussian y ∼ N Ψ−1 u, αΨ−1 (Ψ−1 )T .

4.1. Definition of complete data
When applying the EM algorithm to the current MAP estimation problem, the first step is to define the “complete data”. We2599

denote the complete data as x := [y T , ›T ]T . Taking the conditional expectation of the log-likelihood given x and Θ = Θand then adding log P (Θ), we obtain the Q-function as follows:
“
”
h
c 1
Q(Θ, Θ ) =log |Λ−1 | − tr Λ−1 E[xxT |o; Θ ]2
i
T −1
+ 2m Λ E[x|o; Θ ] − mT Λ−1 m + log P (Θ), (8)„» –» −1 –
–«
»1 T
Ψ Ψ
0
y
Ψ u
−1
α
x :=
, m :=,, Λ :=
1
›
I
0
0
β N
where the many-to-one relationship between the complete data
x and the incomplete data o is described as o = Hx with
H := [I N I N ]. E[x|o; Θ ] and E[xxT |o; Θ ] are describedbyE[x|o; Θ ] = m + ΛH T (HΛH T )−1 (o − Hm), (9)
(10)

u=

(12)

I
N X
(i)
” X
wi Υn,nwi Υ(i)γi,n log
,
≥
n,n
γi,n
n=1 i=1i=1(13)(0 < p ≤ 1),(14)

with auxiliary variables w̄ := {w̄1 , w̄2 , . . . , w̄I } and ‚ :={γ1,1 , . . . , γI,N }. Hence, we define the auxiliary function
f + (w, u, β, w̄, ‚)
I
N X
(i)
X
N
wi Υn,n1γi,n log+ uT Ψx̄ylog αβ +2γαi,nn=1 i=11
1
1 T
tr(ΨT ΨRy ) −tr(R ) −
u u
2α
2β
2α
I “
”
X
p|w̄i |p−1 wi + |w̄i |p − p|w̄i |p ,
− λp

1 ` ´
tr R .
N(17)

The effectiveness of using the proposed method to decompose
F0 contours into musical-note components and dynamics components is evaluated experimentally. We compare the parameter estimation accuracy of three methods, i.e., the conventional
method based on the maximum-likelihood approach [13] (Conventional), the proposed method using the difference approximation (Method 1), and the proposed method using the impulse
response approximation (Method 2). a and b were initially set
at a = 100 and b = 20 calculated by ζ = 1.0 and Ω = 0.1.
Using ζ values that varied from 0.01 to 2 at 0.02 intervals andΩ values that varied from 0.05 to 3 at 0.005 intervals, we calculated Υ (I = 5100). α, β, and wi were initially set at α = 2,
β = 100, and wi = 1/I. u was initially set at the median valueof the observed signal o. λ and p were set at λ = 10000 and
p = 0.8. We determined these initial parameters empirically.
For the first experiment we evaluated the step signal and theimpulse response decomposed from the F0 contour to determine
whether the proposed methods were affected by local minimum
problem. First, we synthesized one hundred F0 contours artificially using randomly-determined ζ, Ω, and u values, where N ,α, and β were set at 300, 2, and 100, respectively. Then, we
estimated model parameters for each F0 contour. Finally, we
calculated the root mean square error (RMSE) between the step
signals, which we obtained by using the original u and the estimated parameter u, respectively. We also calculated the RMSE
between the impulse responses, which we obtained by using theoriginal ζ and Ω values and the model parameters, respectively.
Tab. 1 shows the average values of the RMSEs between thesesignals over all segments. The best value was obtained with
Method 2. We consider that Method 2 solves the local minimum problem more effectively.

(i)

:= −

β=

5. Experiment

where Υn,n is the nth diagonal component of the matrix Υ(i) .
We utilize the auxiliary function method [16] to maximize
f (w, u, β). To define an auxiliary function of Eq. (12), we
use two inequalities:

n=1|wi |p ≤ p|w̄i |p−1 wi + |w̄i |p − p|w̄i |p

1 T
1 Ψx̄y ,N

If the convergence criterion of Eq. (12) is not satisfied, then letΘ ← Θ and return to the E-step. We can similarly derive theM-step for the difference approximation method Eq. (2) basedon the auxiliary function method.

N
I
“X
”
X
N
log αβ +log
wi Υ(i)
n,n
2
n=1
i=1

I
“X

(16)

M-step: Update w, u, and β according to Eqs. (16) and (17)

I

log

= 0.

E-step: Evaluate E[x|o; Θ ] and E[xxT |o; Θ ] using the current parameter values Θ and update w̄ and ‚.

1
1
1+ uT Ψx̄y −tr(ΨT ΨRy ) −tr(R )
α
2α
2β

N
X

wi

Initial step: Initialize Θ (w, u, and β).

We can derive the M-step update formulae for all the model parameters by employing the definitions of the complete data andthe prior. In this subsection we derive the M-step update formulae for the impulse response approximation method describedin Section 2 owing to space limitations, where Ψ correspondsto Φ−1 of Eq. (4). Collecting terms that depend on Θ from Eq.
(8), we obtain the objective function as follows:X
1 T|wi |p ,u u − λp
2α
i=1

n=1

The proposed algorithm is summarized as follows:

4.2. M-step update formulae−NXγi ,nSolving this nonlinear simultaneous equation for i = 1, . . . , I
based on the coordinate descent method, we obtain w. Update
rules for u and β are the followings:

For convenience, we segment E[x|o; Θ ] and E[xxT |o; Θ ]into small sections as follows:
–
» –
»
∗
x̄y
Ry
, E[xxT |o; Θ ] =. (11)
E[x|o; Θ ] =
x̄
∗
R

f (w, u, β) := −I
1 X “ T (i) T (i ) ”
1
tr Ry ΥΥwi − uT Υ(i ) x̄yα i=1α+ λp p|w̄i |p−2 wi −E[xxT |o; Θ ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|o; Θ ]E[x|o; Θ ]T .

where f (w, u, β) ≥ f + (w, u, β, w̄, ‚) is satisfied. We have
P
(i)
(i )
equality if w̄i = wi and γi,n = w̄i Υn,n / i w̄i Υn,n , hence
Eq. (15) satisfies the definition of the auxiliary function method.
By iteratively updating E[x|o; Θ ], E[xxT |o; Θ ], Θ, w̄, and‚, f (w, u, β) will converge to a stationary point. As the update
rules for w̄ and ‚ are shown above, we only need to derive theupdate rule for Θ. Differentiating Eq. (15) partially w.r.t. wi
and setting to zero, we obtain

−

(15)

i=1

2600

Table 1: Root mean square error (RMSE) between signals

Conventional
Method 1
Method 2

21.7
1.94
0.67

o and —
[cent]
44.7
42.9
34.7

Impulseresponses
5.29×10−2
65.2×10−4
2.98×10−4Classical singer(Female)0.18

ΩStep signals

0.20

Shorterrise time

0.16

Amateur(Female)0.14
Longerrise time

0.12

Classical singer(Male)Pop singer(Female)Amateur(Male)
0.46

0.50

Pop singer(Male)

0.54

0.58

0.62

ζ
Over-damping
Figure 3: Estimated ζ and Ω for each singer
Under-damping

60

5000
4600
4200

O
0

2

4

6

8

10

Time [sec]

40
20
0

Figure 2: Comparing observed signal o with signal — generated
by model parameters of Method 2

Classical singers-400 -200 0 200 400
u [cent]

Amateurs

60
Frequency

Classical singer(Female)5400

Frequency

F0 [cent]

5800

40
20
0

-400 -200 0 200 400
u [cent]

Figure 4: Histograms of u for classical and amateur singers
and derived a model parameter estimation algorithm based on
the EM method. Our experiments showed that the proposed
method can solve the local minimum problem and decompose
an F0 contour into a musical note component and a dynamics
component. Furthermore, we examined the differences between
singing styles based on the estimated parameters. In the future,
we plan to divide the F0 contour into segments using Hidden
Markov Model (HMM), and evaluate its ability to detect particular singing behaviors automatically such as vibrato and overshoot and singing voice synthesis. We also want to employ the
proposed method for automatic speech recognition and gesturerecognition.

For the second experiment we evaluated the convergence
performance of the proposed methods using real F0 contours.
We used singing signals provided by six singers. We recorded
signals for each gender in the professional classical, professional pop, and amateur categories. Each subject sang songs
with Japanese lyrics (two patterns) and hummed, in both cases
without musical accompaniment. The songs were “Twinkle,
Twinkle, Little Star” and “Ode to Joy”. In all, 36 singing signals were recorded. The F0 contour was estimated every 5 ms
(Δ = 5ms) using YIN [17] and represented in cents, so that
one-tempered semitone corresponded to 100 cents. Using thelinear interpolation, we smoothed F0 in voiceless consonants.
We then divided the F0 contours into 1323 segments manually
as seen in Fig. 1 and subtracted the starting F0 value over eachsegment so that the starting F0 value of the segment was zero.
Finally, we estimated the model parameters for each segment.
The average values of the RMSE between an observed signal o and the signal — generated by model parameters are shownin Tab. 1. We obtained the smallest RMSE for Method 2.
The temporal attack of the musical note was successfully represented by the model parameters (see Fig. 2). Hence, the RMSE
corresponds to a residual component such as vibrato. Since
Method 2 represents the dynamics by a linear combination ofvarious oscillation bases and regularizes weighting parametersw, we consider that over-fitting for the model is controlled andthe RMSE is the smallest among three methods. Furthermore,
the periodic components such as vibrato are included in theresidue ›. Instead of modeling › by Gaussian white noise, weplan to model › using the periodic kernel similarly with GPs.
Fig. 3 shows ζ and Ω for each singer. Estimating w foreach segment using Method 2, we plotted the average values ofζ and Ω corresponding to the maximum value of the elements ofw over all segments. If ζ and Ω are small, this means that theoscillation phenomenon is Under-damping such as overshootand the rise time of the musical note is longer. Since amateurs
have poor singing skills, ζ and Ω are small. Fig. 4 shows histograms of parameter u for classical and amateur singers. We
confirmed that the histogram has peaks at the integral multiplesof the semitone for classical and pop singers. However, the histogram peaks for amateurs are blurred. This means that it isdifficult for amateurs to sing at the correct pitch. In the future,
we plan to analyze singing styles based on estimated ζ, Ω, andu values using a large singing database.

7. References[1] T. Saitou et al., “Speech-To-Singing Synthesis: Converting speaking voices
to singing voices by controlling acoustic features unique to singing voices,”in WASSPA 2007, pp. 215–218.
[2] ——, “Acoustic and perceptual effects of vocal training in amateur malesinging,” in EUROSPEECH 2009, pp. 832–835.
[3] T. Kako et al., “Automatic identification for singing style based on sung
melodic contour characterized in phase plane,” in Proc. ISMIR 2009, pp.
393–397.
[4] T. Nakano et al., “An automatic singing skill evaluation method for unknown melodies using pitch interval accuracy and vibrato features,” in Proc.
ICSLP 2006, pp. 1706–1709.
[5] J. Sundberg, “The KTH synthesis of singing,” Advances in Cognitive Psychology. Special issue on Music Performance, vol. 2, pp. 131–143, 2006.
[6] J. Bonada and A. Loscos, “Sample-based singing voice synthesizer by spectral concatenation,” in Proc. SMAC 2003.
[7] T. Nakano et al., “VocaListener: A Singing-to-Singing synthesis system
based on iterative parameter estimation,” in Proc. SMC 2009, pp. 343–348.
[8] S. Fukayama et al., “Orpheus: Automatic composition system considering
prosody of Japanese lyrics,” in Proc. ICEC 2009, pp. 309–310.
[9] R. B. Dannenberg, W. P. Birmingham, et al., “A comparative evaluationof search techniques for Query-by-Humming using the MUSART testbed,”
JASIST, vol. 58, no. 5, pp. 687–701, 2007.
[10] P. Proutskova and M. Casey, “You call That singing? ensemble classification
for multi-cultural collections of music recordings,” in Proc. ISMIR 2009, pp.
759–764.
[11] N. Minematsu et al., “Prosodic modeling of Nagauta singing and its evaluation,” in Proc. SpeechProsody 2004, pp. 487–490.
[12] H. Fujisaki, “A note on the physiological and physical basis for the phrase
and accent components in the voice fundamental frequency contour,” inVocal Physiology: Voice Production, Mechanisms and Functions, (O. Fujimura, ed.). Raven Press, 1988, pp. 347–355.
[13] Y. Ohishi et al., “Parameter estimation method of F0 control model forsinging voices,” in Proc. ICSLP 2008, pp. 139–142.
[14] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine
Learning. MIT Press, Cambridge, Mass, USA, 2006.
[15] F. Bach et al., “Multiple kernel learning, conic duality, and the smo algorithm,” in Proc. ICML 2004, pp. 6–13.
[16] H. Kameoka et al., “Complex NMF: A new sparse representation for acoustic signals,” in Proc. ICASSP 2009, pp. 3437–3440.
[17] A. de Cheveigné and H. Kawahara, “YIN, a fundamental frequency estimator for speech and music,” JASA, vol. 111, no. 4, pp. 1917–1930, 2002.

6. Concluding remarks
This paper proposed a statistical representation of F0 dynamics
in singing voices based on GPs with multiple oscillation bases,

2601