ICASSP 2001 (The 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing) Proceedings, pp. V-3365-3368,
May 2001.

A PREDOMINANT-F0 ESTIMATION METHOD FOR CD RECORDINGS:
MAP ESTIMATION USING EM ALGORITHM FOR ADAPTIVE TONE MODELS
Masataka Goto
“Information and Human Activity”, PRESTO, Japan Science and Technology Corporation (JST). /Electrotechnical Laboratory (National Institute of Advanced Industrial Science and Technology).
1-1-4 Umezono, Tsukuba, Ibaraki 305-8568, JAPANgoto@etl.go.jp
nant-F0 Estimation method), that can detect the melody and basslines in complex mixtures containing simultaneous sounds of various instruments (even drums) [9]. PreFEst has the advantages thatit does not assume the number of sound sources, it does not locallytrace frequency components, and it does not rely on the existence of
the F0’s frequency component. It basically estimates the F0 of themost predominant harmonic structure in the input sound mixture;
it simultaneously takes into consideration all the possibilities of F0
and considers that the input mixture contains every possible harmonic structure with different weights (amplitude). It regards the
input frequency components as a weighted mixture of harmonicstructure tone models of all possible F0s and then finds the F0 of
the maximum-weight model corresponding to the most predominant harmonic structure.
PreFEst reported in our earlier paper [9] had three limitations.
First, although various kinds of harmonic structure appear at different F0s and even at the same F0, just a single harmonic-structure
tone model was prepared for each F0. Second, the shape of tonemodel was fixed as if one ideal tone model was always assumed:
the relative amplitude of each harmonic component was constant.
Because it did not exactly coincide with the harmonic structure
contained in the input, there was room for refining the tone-model
management in a more adaptive way. Third, even if prior knowledge about very rough F0 estimates of the melody and bass lines
is available, we were not able to incorporate it into the estimation.
Such prior rough estimates can be given in some practical applications where more precise F0 with less errors is required; for example, the analysis of expression in a recorded performance needs toestimate the actual F0 by using its rough estimate that can be given
by a score or by playing a MIDI instrument along to the originalrecorded performance.
In the following sections we describe how we extended ourprevious PreFEst so that it can overcome the above three limitations. We first give an overview of PreFEst described in [9] andthen describe the details of the three extensions that make it possible to overcome the above limitations. The main idea is to introduce multiple adaptive tone models and prior knowledge, and to
estimate the model parameters on the basis of the MAP (Maximum
A Posteriori Probability) estimation by using the ExpectationMaximization (EM) algorithm. Finally, we show experimental results of a real-time system based on the extended method.

ABSTRACTThis paper describes a predominant-F0 (fundamental frequency) estimation method called PreFEst, which can detect
melody and bass lines in monaural audio signals containing soundsof various instruments. While most previous methods premised
mixtures of a few sounds and had difficulty dealing with such complex signals, our method can estimate the F0 of the melody and basslines without assuming the number of sound sources in compactdisc recordings. In this paper we propose the following three extensions to our previous PreFEst to make it more adaptive and flexible:
introducing multiple harmonic-structure tone models, estimating
the shape of tone models, and introducing a prior distribution of its
shape and F0 estimates. These extensions were implemented by the
MAP (Maximum A Posteriori Probability) estimation by using the
Expectation-Maximization algorithm. Experimental results withcompact-disc recordings showed that our real-time system basedon the extended PreFEst achieved performance improvement.
1. INTRODUCTION
Our goal is to build a real-time system that can detect melody andbass lines in monaural complex real-world audio signals, such as
those sampled from commercially distributed compact discs. Thisdetection is an important initial step in computer emulation of human music understanding because the melody and bass lines arefundamental to the perception of Western music. In addition, thedetected melody and bass lines are useful in various practical applications, such as automatic music indexing for information retrieval(e.g., searching for a song by singing a melody), computer participation in live human performances, and analysis of recordings ofoutstanding performances.
Although this detection requires the estimation of the fundamental frequency (F0, perceived as pitch) of the melody and basslines, it has been considered very difficult to estimate the F0 incomplex audio signals sampled from compact discs. The main reasons are: in compact-disc recordings, the number of sound sources
cannot be assumed, the frequency components of one sound often overlap frequency components of simultaneous sounds, and theF0’s frequency component (the frequency component corresponding to the F0) is sometimes missing or very weak (missing fundamental). Most previous F0-estimation methods [1, 2, 3, 4], however, premised that the input contained just a single-pitch soundwith aperiodic noises. Although several methods for dealing withmultiple-pitch mixtures were proposed [5, 6, 7, 8], they dealt withat most three musical instruments or voices and had difficulty dealing with compact-disc recordings.
We therefore developed a method, called PreFEst (Predomi-

c0-7803-7041-4/01/$10.00 2001
IEEE

2. OVERVIEW OF PREFEST
PreFEst estimates the most predominant F0 in frequency-rangelimited sound mixtures. Since the melody line tends to have3365

the most predominant harmonic structure in middle- and highfrequency regions and the bass line tends to have the most predominant harmonic structure in a low-frequency region, we canestimate the F0s of the melody and bass lines by applying PreFEst
with appropriate frequency-range limitation [9, 10].
To use statistical methods, we represent the input bandpassfiltered frequency components as a probability density function(PDF), called observed PDF, which is estimated by using multirate
signal-processing techniques and instantaneous-frequency-related
measure [9, 10]. We then consider that the observed PDF has been
generated from a weighted-mixture model of tone models of all
possible F0s; a tone model is the PDF corresponding to a typical harmonic structure of every possible F0. Because the weights
of tone models represent the relative dominance of every possible
harmonic structure, we can regard those weights as the PDF of theF0 (F0’s PDF): the more dominant a tone model in the mixture,
the higher the probability of the F0 of its model. As explained inour earlier papers [9, 10], those weights (i.e., the F0’s PDF) can be
estimated by using the EM algorithm [11].
A simple way of determining the most predominant F0 is tofind the frequency that maximizes the F0’s PDF. This result is notstable, however, because peaks corresponding to the F0s of simultaneous sounds sometimes compete in the F0’s PDF for a moment
and are transiently selected, one after another, as the maximum. We
therefore consider the global temporal continuity of the F0 by using
a multiple-agent architecture in which agents track different temporal trajectories of the F0, and the final F0 output is determinedon the basis of the most dominant and stable F0 trajectory.

p(x,1 | F,2,µ(t)(F,2))

Tone Model (m=2)
p(x | F,2,µ(t)(F,2))(t)c (1|F,2)c(t)(2|F,2)

Tone Model (m=1)
p(x | F,1,µ(t)(F,1))c(t)(3|F,2)(t)c (1|F,1)c(t)(4|F,2)c(t)(2|F,1)c(t)(3|F,1)c(t)(4|F,1)

F

fundamental
frequency F+1200

F+1902
F+2400x [cent]x [cent]

Figure 1: Model parameters of multiple adaptive tone models.
p(x|F, m, µ(t) (F, m)) (Figure 1). The number of tone models is
Mi (1 ≤ m ≤ Mi ) where i denotes the melody line (i = m) or thebass line (i = b). The tone model that indicates where the harmonics of the F0, F , tend to occur is defined asHip(x|F, m, µ(t) (F, m)) =p(x, h|F, m, µ(t) (F, m)),

(2)

h=1

p(x, h|F, m, µ(t) (F, m))= c(t) (h|F, m) G(x; F + 1200 log2 h, Wi ),µ (F, m) = {c (h|F, m) | h = 1, ..., Hi },
(t)

(t)

(3)
(4)

(x−x0 )2
−
2σ 2

1
e
,
(5)
2πσ 2
where Hi is the number of harmonics considered, W2i is the variance of the Gaussian distribution G(x; x0 , σ), and c(t) (h|F, m) determines the relative amplitude of the h-th harmonic component(the shape of tone model) and satisfies
G(x; x0 , σ) = √3. THREE EXTENSIONS
To overcome the three limitations described in the Introduction, we
propose the following three extensions of the previous PreFEst.
[Extension 1] Introducing multiple tone models
We prepare multiple tone models for each F0 and consider their
mixture model.
[Extension 2] Estimating the shape of tone models
We consider, as model parameters, the relative amplitude ofeach harmonic component of all the tone models (called theshape of tone model) as well as their weights, and estimate themby using the EM algorithm.
[Extension 3] Introducing a prior distribution
To estimate the model parameters on the basis of their prior distribution, we use the MAP estimation instead of the maximum
likelihood estimation used in [9, 10]. Note that we can also take
into consideration a prior distribution of the shape of tone model
added to the model parameters by [Extension 2].
All of these extensions are dealt with in the process of estimating the F0’s PDF p(t)F 0 (F ) from the observed PDF (the PDF of the
bandpass-filtered frequency components), p(t)
Ψ (x). Here, t is thetime measured in units of frame-shift (10 msec), and x and F are
the log-scale frequency denoted in units of cents (a musical-intervalmeasurement). Frequency fHz in hertz is converted to frequency
fcent in cents as follows:
fHz
.
(1)
fcent = 1200 log2
3
440 × 2 12 −5

Hi

(t)

c (h|F, m) = 1.

(6)

h=1

We then consider that the observed PDF was generated fromthe following model p(x|θ(t) ) that is a weighted mixture of all possible tone models p(x|F, m, µ(t) (F, m)):(t)p(x|θ ) =

 Fhi Mi
Fliw(t) (F, m) p(x|F, m, µ(t) (F, m)) dF, (7)m=1θ (t) = {w (t) , µ(t) },(8)w (t) = {w (t) (F, m) | Fli ≤ F ≤ Fhi , m = 1, ..., Mi },(9)µ(t) = {µ(t) (F, m) | Fli ≤ F ≤ Fhi , m = 1, ..., Mi }, (10)where Fli and Fhi denote the lower and upper limits of the possible(allowable) F0 range and w (t) (F, m) is the weight of a tone modelp(x|F, m, µ(t) (F, m)) which satisfies

 Fhi Mi
Fliw(t) (F, m) dF = 1.

(11)

m=1

Because we cannot know a priori the number of sound sources,
it is important that we simultaneously take into consideration all
the possibilities of the F0 as expressed in Equation (7). If we
can estimate the model parameter θ (t) such that the observed PDF(t)p(t)Ψ (x) is likely to have been generated from the model p(x|θ ),(t)(t)
the weight w (F, m) can be interpreted as the F0’s PDF pF 0 (F )because w (t) (F, m) represents the relative dominance of the harmonic structure:

3.1. Weighted-mixture model of extended tone models
For [Extension 1] and [Extension 2], we prepare multiple tonemodels for each F0, F , and introduce the model parameterµ(t) (F, m) to the m-th tone model whose PDF is denoted asp(t)F 0 (F ) =

Mi

m=1

3366

w(t) (F, m) (Fli ≤ F ≤ Fhi ).

(12)

Q(θ (t) |θ(t) ) =3.2. Introducing a prior distribution

 ∞
−∞

(t)

For [Extension 3], we define a prior distribution p0i (θ ) of themodel parameter θ(t) as follows:(13)
p0i (θ(t) ) = p0i (w(t) ) p0i (µ(t) ),(t)
1 −βwi(t)
Dw (w(t) ;w (t) )0i
p0i (w ) =
e
,
(14)
Zw

 Fh Mi

1 − Fl ii
p0i (µ ) =
e
Zµ
(t)

m=1

(t)

(20)
EF,m,h [log p(x, F, m, h|θ (t) ) | x, θ(t) ] dx,
where Q(θ(t) |θ(t) ) is the conditional expectation of themean log-likelihood for the maximum likelihood estimation.
EF,m,h [a|b] denotes the conditional expectation of a with respect to the hidden variables F , m, and h with the probability
distribution determined by condition b.
2. (M-step)
Maximize QMAP (θ(t) |θ (t) ) as a function of θ(t) in order to obtainthe updated (improved) estimate θ(t) :θ(t) = argmax QMAP (θ(t) |θ(t) ).
(21)

(t)

βµi (F,m) Dµ (µ0i (F,m);µ(t) (F,m)) dF.

(15)Here, p0i (w(t) ) and p0i (µ(t) ) are unimodal distributions; p0i (w(t) )(t)(F, m) and p0i (µ(t) ) takes thetakes the maximum value at w0i(t)(t)(F, m) and µ(t)maximum value at µ0i (F, m), where w0i
0i (F, m)(t)(c0i (h|F, m)) are the most probable parameters. Zw and Zµ are(t)(t)
the normalization factors, and βwi
and βµi(F, m) are the parameters determining how much emphasis is put on the maximum value;(t)(t)= 0 and βµi(F, m) = 0 represent the noninformative prior
βwi
distribution (uniform distribution). In Equations (14) and (15),(t)(t); w(t) ) and Dµ (µ(t)
Dw (w0i
0i (F, m); µ (F, m)) are defined as thefollowing Kullback-Leibler’s information:(t); w(t) ) =Dw (w0i Fhi 
Mi
Fli

(t)
w0i
(F, m) log

m=1

θ (t)

In the E-step, Q(θ (t) |θ(t) ) is expressed as
Q(θ (t) |θ(t) ) =

=

h=1

(22)

= log(w (t) (F, m) p(x, h|F, m, µ(t) (F, m))).
From Equation (13) the log prior distribution is given bylog p0i (θ(t) ) = − log Zw Zµ

(23)

(17)

 Fhi Mi (t) (t)βwi
w0i (F, m) log

Fli(t)(F, m)
+ βµi

m=1
Hi

(t)

c0i (h|F, m) log

h=1

(t)
(F, m)
w0i
w(t) (F, m)c(t)0i (h|F, m)c(t) (h|F, m)
dF. (24)

Regarding the M-step, Equation (21) is a conditional problemof variation, where the conditions are given by Equations (6) and(11). This problem can be solved by using the following EulerLagrange differential equations with Lagrange multipliers λw andλµ :The problem to be solved is to estimate the model parameter θ (t) on
the basis of the prior distribution p0i (θ(t) ) when we observe p(t)
Ψ (x).
The MAP (Maximum A Posteriori Probability) estimator of θ (t) is
obtained by maximizing∞−∞

m=1 h=1p(F, m, h|x, θ ) log p(x, F, m, h|θ (t) )dF dx,
where the complete-data log-likelihood is given bylog p(x, F, m, h|θ(t) )−3.3. MAP estimation using EM algorithm

(t)
(t)
p(t)
Ψ (x) (log p(x|θ ) + log p0i (θ )) dx.

Fli

p(t)
Ψ (x)

(t)

(t)
Dµ (µ(t)0i (F, m); µ (F, m))c(t) (h|F, m).
c0i (h|F, m) log 0ic(t) (h|F, m)

 ∞  Fhi 
Mi Hi

−∞

(t)
(F, m)
w0i
dF, (16)w(t) (F, m)

Hi

(t)

p(t)
Ψ (x)

(18)

∂
∂w (t)

Because this maximization problem is too difficult to be solved analytically, we use the EM algorithm to estimate θ (t) . While the EM
algorithm is usually used for computing maximum likelihood estimates from incomplete observed data, it can also be used for computing MAP estimates as described in [11]. In the maximum likelihood estimation, the EM algorithm iteratively applies two steps,
the expectation step (E-step) computing the conditional expectation of the mean log-likelihood and the maximization step (M-step)
maximizing its expectation. On the other hand, in the MAP estimation, it iteratively applies the E-step computing the sum of theconditional expectation and the log prior distribution and the Mstep maximizing it. With respect to θ (t) , each iteration updates the‘old’ estimate θ(t) = {w(t) , µ(t) } to obtain the ‘new’ improved
estimate θ(t) = {w(t) , µ(t) }.
By introducing hidden (unobservable) variables F , m, andh, which respectively describe which F0, which tone model, andwhich harmonic component were responsible for generating eachobserved frequency component at x, we can specify the two steps
as follows:
1. (E-step)Compute the following QMAP (θ(t) |θ(t) ) for the MAP estimation:(19)
QMAP (θ(t) |θ(t) ) = Q(θ (t) |θ(t) ) + log p0i (θ(t) ), ∞ 
Hi−∞ h=1

(t)p(t))Ψ (x) p(F, m, h|x, θ(log w(t) (F, m) + log p(x, h|F, m, µ(t) (F, m))) dx
(t)
(t)
w0i
(F, m) log− βwi− λw (w
∂
∂c(t)

 ∞
−∞

(t)

(t)
(F, m)
w0i
w(t) (F, m)

(F, m) − Mi (Fh1i −Fli ) )

= 0,

(25)

(t)
p(t)
) (log w(t) (F, m)
Ψ (x) p(F, m, h|x, θ+ log c(t) (h|F, m) + log G(x; F + 1200 log2 h, Wi )) dx(t)(F, m) c(t)− βµi0i (h|F, m) logc(t)0i (h|F, m)c(t) (h|F, m)− λµ (c(t) (h|F, m) − H1i )From these equations we get
 ∞
1
(t)
p(t)
) dx
w(t) (F, m) =
Ψ (x) p(F, m|x, θ
λw
−∞

= 0.



(t) (t)w0i (F, m) ,
+ βwi

3367

(26)

(27)

c(t) (h|F, m) =

1
λµ

 ∞
−∞

The system was tested on excerpts of 10 songs in popular, jazz,and orchestral genres. The input monaural audio signals — eachcontaining a single-tone melody and the sounds of several instruments — were sampled from compact discs. We evaluated the detection rates by comparing the estimated F0s with the correct F0s
that were hand-labeled by using the F0 editor program we previously developed [9].
In our experiment the system correctly detected the melody
and bass lines for most of each audio sample; the average detection rate for the melody line was 88.4% and that for the bass line
was 79.9%. The results of comparing the extended PreFEst (withthe two adaptive tone models) with the previous PreFEst without
any extension showed that the detection rates for three songs weregreatly improved (at most 29.2% improvement).

(t)
p(t)
) dx
Ψ (x)p(F, m, h|x, θ(t)(F, m)c(t)+ βµi0i (h|F, m) .

(28)

In these equations, λw and λµ are determined from Equations (6)and (11) as(t)
λw = 1 + βwi
,
(29)

 ∞

λµ =
−∞

(t)
(t)
p(t)
) dx + βµi(F, m).
Ψ (x)p(F, m|x, θ(30)

From the Bayes’ theorem, p(F, m, h|x, θ (t) ) and p(F, m|x, θ(t) )are given byw(t) (F, m) p(x, h|F, m, µ(t) (F, m)), (31)
p(F, m, h|x, θ(t) ) =p(x|θ (t) )5. CONCLUSIONw(t) (F, m) p(x|F, m, µ(t) (F, m)). (32)
p(x|θ (t) )Finally we obtain the ‘new’ parameter estimates w(t) (F, m) andc(t) (h|F, m):p(F, m|x, θ (t) ) =w(t) (F, m) =(t)(t) (t)wML(F, m) + βwi
w0i (F, m)(t)
1 + βwi,We have described how to improve a method called PreFEst that
detects the melody and bass lines in complex real-world audio signals. The MAP estimation executed by using the EM algorithm enabled three extensions: introducing multiple tone models, estimating the shape of tone models, and introducing a prior distributionof its shape and F0 estimates. Experimental results showed that areal-time system using the extended PreFEst improved in performance and is robust enough to estimate the F0s of the melody andbass lines in compact-disc recordings.
Although the three extensions we made have great potential,
we have not fully exploited them. In the future, for example, a
lot of tone models could be prepared by analyzing various kindsof harmonic structure appearing in a music (sound) database. We
also plan to report experimental results that have shown the effectiveness of using a prior distribution of F0 estimates.

(33)c(t) (h|F, m) =
(t)
(t)
(t)
wML
(F, m) c(t)ML (h|F, m) + βµi (F, m)c0i (h|F, m)
(t)
(t)
wML
(F, m) + βµi(F, m), (34)(t)(F, m) and c(t)where wML
ML (h|F, m) are the following maximum
likelihood estimates when the noninformative prior distribution
(t)
(t)
(βwi
= 0 and βµi(F, m) = 0) is given:

 ∞

(t)
wML
(F, m) =−∞ACKNOWLEDGMENTS: I thank Shotaro Akaho and Hideki
Asoh for their valuable discussions.

p(t)
Ψ (x)

w(t) (F, m) p(x|F, m, µ(t) (F, m)) Fhi Mi

ν=1

Fli

c(t)
ML (h|F, m) =
w

(t)

1
(t)
wML
(F, m)

ν=1

6. REFERENCES

dx, (35)

 ∞−∞[1] L. R. Rabiner et al., “A comparative performance study of several
pitch detection algorithms,” IEEE Trans. on ASSP, vol. ASSP-24,
no. 5, pp. 399–418, 1976.
[2] A. Nehorai and B. Porat, “Adaptive comb filtering for harmonic
signal enhancement,” IEEE Trans. on ASSP, vol. ASSP-34, no. 5,
pp. 1124–1138, 1986.
[3] T. Abe et al., “Robust pitch estimation with harmonics enhancement
in noisy environments based on instantaneous frequency,” in ICSLP
96, pp. 1277–1280, 1996.
[4] H. Kawahara et al., “Fixed point analysis of frequency to instantaneous frequency mapping for accurate estimation of F0 and periodicity,” in Eurospeech 99, pp. 2781–2784, 1999.
[5] C. Chafe and D. Jaffe, “Source separation and note identification inpolyphonic music,” in ICASSP 86, pp. 1289–1292, 1986.
[6] H. Katayose and S. Inokuchi, “The kansei music system,” Computer
Music Journal, vol. 13, no. 4, pp. 72–77, 1989.
[7] G. J. Brown and M. Cooke, “Perceptual grouping of musical sounds:
A computational model,” JNMR, vol. 23, pp. 107–132, 1994.
[8] K. Kashino and H. Murase, “Music recognition using note transition
context,” in ICASSP 98, pp. 3593–3596, 1998.
[9] M. Goto, “A robust predominant-f0 estimation method for real-time
detection of melody and bass lines in CD recordings,” in ICASSP
2000, pp. II–757–760, 2000.
[10] M. Goto, “A real-time music scene description system: Detecting
melody and bass lines in audio signals,” in Working Notes of the
IJCAI-99 Workshop on CASA, pp. 31–40, 1999.
[11] A. P. Dempster et al., “Maximum likelihood from incomplete data
via the EM algorithm,” J. Roy. Stat. Soc. B, vol. 39, no. 1, pp. 1–38,
1977.

p(t)
Ψ (x)

(F, m) p(x, h|F, m, µ(t) (F, m)) Fhi Mi
Fliw(t) (η, ν) p(x|η, ν, µ(t) (F, ν)) dηw(t) (η, ν) p(x|η, ν, µ(t) (F, ν)) dηdx. (36)

After the above iterative computation, the F0’s PDF p(t)F 0 (F )estimated by considering the prior distribution can be obtained
from w(t) (F, m) according to Equation (12). We can also obtainc(t) (h|F, m), which is the relative amplitude of each harmoniccomponent of all the multiple tone models p(x|F, m, µ(t) (F, m)).
All three extensions are thus fulfilled.
4. EXPERIMENTAL RESULTS
The extended PreFEst has been implemented in a real-time systemthat takes a musical audio signal as input and outputs the detected
melody and bass lines in several forms, such as audio signals forauralization and computer graphics for visualization [9, 10]. The
current implementation uses the following parameter values with
two adaptive tone models: Fhm = 8400 cent, Flm = 3600 cent, Mm =2, Hm = 16, Wm = 17 cent, Fhb = 4800 cent, Flb = 1000 cent, Mb = 2,
Hb = 6, and Wb = 17 cent. For the prior distribution of the shape of
tone models, we use c(t)0i (h|F, m) = αi,m gm,h G(h; 1, Ui ), wherem is 1 or 2, αi,m is a normalization factor, gm,h is 2/3 (m = 2 andh is even) or 1 (otherwise), Um = 5.5, and Ub = 2.7.

3368