Proceedings of the 18th International Congress on Acoustics (ICA 2004), pp.II-1085-1088, April 2004. (Invited Paper)

A Predominant-F0 Estimation Method for Polyphonic Musical Audio Signals
Masataka Goto
National Institute of Advanced Industrial Science and Technology (AIST)
m.goto@aist.go.jp

Abstract
In this paper I introduce a method, called PreFEst, for estimating the fundamental frequency (F0) of simultaneous
sounds in monaural polyphonic audio signals. Most previous F0-estimation methods have had difficulty dealing
with such complex audio signals because these methods
were designed to deal with mixtures of only a few sounds.
Without assuming the number of sound sources, PreFEst
can estimate the relative dominance of every possible harmonic structure in the input mixture. It treats the mixture
as if it contains all possible harmonic structures with different weights, and estimates their weights by MAP estimation. PreFEst can obtain the melody and bass lines byregarding the most predominant F0 in middle- and highfrequency regions as the melody line and the one in a lowfrequency region as the bass line. Experimental resultswith compact-disc recordings showed that a real-time system implementing this method was able to detect melody
and bass lines about 80% of the time these existed.

mates the F0 of the most predominant harmonic structure— the most predominant F0 corresponding to the melody
or bass line — within an intentionally limited frequency
range of the input mixture. It simultaneously takes into
consideration all possibilities for the F0 and treats the input mixture as if it contains all possible harmonic structures with different weights (amplitudes). It regards a
probability density function (PDF) of the input frequency
components as a weighted mixture of harmonic-structure
tone models (represented by PDFs) of all possible F0s andsimultaneously estimates both their weights corresponding to the relative dominance of every possible harmonic
structure and the shape of the tone models by MAP (Maximum A Posteriori Probability) estimation considering
their prior distribution. It then considers the maximumweight model as the most predominant harmonic structure and obtains its F0. The method also considers theF0’s temporal continuity by using a multiple-agent architecture.

2. Predominant-F0 estimation method:
PreFEst

1. Introduction

PreFEst-front-end

Instantaneous frequency calculation

Extracting frequency components
Limiting frequency regions

BPF for melody line

BPF for bass line

Frequency components
Forming F0's probability density function
F0's probability density function
Agent

Salience
detector

Int

Agent

era

cti

on

Agent

Most predominant F0 trajectory
Melody line

Bass line

Figure 1: Overview of PreFEst.

PreFEst-core PreFEst-back-end

II - 1085

Audio signals

Detecting melody line

Tu2.H.4

Figure 1 shows an overview of PreFEst. PreFEst consistsof three components, the PreFEst-front-end for frequency
analysis, the PreFEst-core to estimate the predominant
F0, and the PreFEst-back-end to evaluate the temporalcontinuity of the F0. Since the melody line tends to havethe most predominant harmonic structure in middle- andDetecting bass line

The estimation of the fundamental frequency (F0) of simultaneous sounds in monaural polyphonic mixtures isimportant to build a computational model that can understand audio signals in a human-like fashion. Moreover,
it is useful in various practical applications, such as automatic indexing for information retrieval and intelligent
editing of audio signals. It has, however, been considered difficult to estimate the F0 in such real-world audio signals; this is because the number of sound sourcesin them generally cannot be assumed, because the frequency components of one sound often overlap the frequency components of simultaneous sounds, and because
the F0’s frequency component (the frequency component
corresponding to the F0) is sometimes very weak or missing (missing fundamental). Most previous F0 estimation
methods have been premised upon the input audio signal
containing just a single-pitch sound with aperiodic noise.
Although several methods for dealing with multiple-pitch
mixtures have been proposed [1, 2, 3, 4], these require thatthe number of simultaneous sounds be assumed and havedifficulty estimating the F0 in complex audio signals likethose sampled from music compact discs.
This paper introduces a Predominant-F0 Estimation
method, called PreFEst, that I developed during 1999 and2001 [5, 6, 7, 8]. PreFEst can estimate the fundamental frequency (F0) of melody and bass lines in monaural
audio signals containing simultaneous sounds of various
musical instruments. Unlike previous methods, PreFEstdoes not assume the number of sound sources, locallytrace frequency components, or even rely on the existenceof the F0’s frequency component. PreFEst basically esti-

high-frequency regions and the bass line tends to have themost predominant harmonic structure in a low-frequency
region, the F0s of the melody and bass lines can be estimated by applying the PreFEst-core with appropriate
frequency-range limitation.

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

2.2. PreFEst-core: Estimating the F0’s probability
density function
For each set of filtered frequency components represented
as an observed PDF p(t)
Ψ (x), the PreFEst-core forms aprobability density function of the F0, called the F0’s
PDF, p(t)F 0 (F ), where F is the log-scale frequency in cents.
We consider each observed PDF to have been generated
from a weighted-mixture model of the tone models of all
the possible F0s; a tone model is the PDF corresponding
to a typical harmonic structure and indicates where theharmonics of the F0 tend to occur. Because the weights
of tone models represent the relative dominance of every possible harmonic structure, these weights can be regarded as the F0’s PDF: the more dominant a tone modelis in the mixture, the higher the probability of the F0 of
its model.
2.2.1. Weighted-mixture model of adaptive tone models
To deal with diversity of the harmonic structure, thePreFEst-core can use several types of harmonic-structure
tone models. The PDF of the m-th tone model for each F0
F is denoted by p(x|F, m, µ(t) (F, m)) (Figure 2), where
the model parameter µ(t) (F, m) represents the shape ofthe tone model. The number of tone models is Mi (1 ≤m ≤ Mi ) where i denotes the melody line (i = m) or thebass line (i = b). Each tone model is defined byHip(x|F, m, µ(t) (F, m)) =p(x, h|F, m, µ(t) (F, m)), (1)h=1(t)= c(t) (h|F, m) G(x; F + 1200 log2 h, Wi ),µ (F, m) = {c(t) (h|F, m) | h = 1, ..., Hi },(x−x0 )21e− 2σ2 ,G(x; x0 , σ) = √2πσ 2c(t)(3|F,2)(t)c (1|F,1)c(t)(4|F,2)c(t)(2|F,1)
F

fundamental
frequency F+1200c(t)(3|F,1)c(t)(4|F,1)
F+1902
F+2400x [cent]x [cent]

Figure 2: Model parameters of multiple adaptive tone models.

where Hi is the number of harmonics considered, Wi
is the standard deviation σ 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 thetone model) and satisfiesHic(t) (h|F, m) = 1.
(5)
h=1

In short, this tone model places a weighted Gaussian distribution at the position of each harmonic component.
We then consider the observed PDF p(t)
Ψ (x) to havebeen generated from the following model p(x|θ(t) ),which is a weighted mixture of all possible tone modelsp(x|F, m, µ(t) (F, m)): Fhi
Mi(t)p(x|θ ) =w (t) (F, m) p(x|F, m, µ(t) (F, m)) dF,
Fli m=1(6)(7)

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

m=1

Because we cannot know a priori the number of soundsources, it is important that we simultaneously take
into consideration all F0 possibilities as expressed in
Equation (6). If we can estimate the model parameter θ (t)such that the observed PDF p(t)
Ψ (x) is likely to have been
generated from the model p(x|θ(t) ), the weight w (t) (F, m)
can be interpreted as the F0’s PDF p(t)F 0 (F ):Miw (t) (F, m) (Fli ≤ F ≤ Fhi ).
(11)
p(t)
F 0 (F ) =p(x, h|F, m, µ (F, m))(t)c (1|F,2)c(t)(2|F,2)

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

2.1. PreFEst-front-end: Forming the observed probability density functions
The PreFEst-front-end first uses an STFT-based multirate filter bank in order to obtain adequate time-frequency
resolution under the real-time constraint. It then extracts frequency components by using an instantaneousfrequency-related measure [5, 8] and obtains two setsof bandpass-filtered frequency components, one for the
melody line (261.6 - 4186 Hz) and the other for the bassline (32.7 - 261.6 Hz). To enable the application of statistical methods, each set of the bandpass-filtered components is represented as a probability density function(PDF), called an observed PDF, p(t)
Ψ (x), where t is thetime measured in units of frame-shifts (10 ms), and x is
the log-scale frequency denoted in units of cents 1 .

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

m=1

(2)
(3)
(4)

1 In this paper I define that frequency f
Hz in hertz is converted to frequency fcent in cents so that there are 100 cents to a tempered semitone3
and 1200 to an octave: fcent = 1200 log2 (fHz / (440 × 2 12 −5 )).

2.2.2. Introducing a prior distribution
To use prior knowledge about F0 estimates and the tonemodel shapes, we define a prior distribution p0i (θ(t) ) of
θ(t) as follows:(12)
p0i (θ(t) ) = p0i (w(t) ) p0i (µ(t) ),1 −β(t) Dw (w(t) ;w(t) )
(t)
0i
e wi
,
(13)
p0i (w ) =Zw

II - 1086

 Fh  i (t)βµi (F,m) Dµ (µ(t)(F,m);µ(t) (F,m)) dF
1 − Fl i Mm=10iip0i (µ ) =
e
.
Zµ
(14)
Here p0i (w(t) ) and p0i (µ(t) ) are unimodal distributions:(t)p0i (w(t) ) takes its maximum value at w0i(F, m) and(t)(t)
p0i (µ ) takes its maximum value at µ0i (F, m), where
(t)
(t)
w0i
(F, m) and µ(t)0i (F, m) (c0i (h|F, m)) are the mostprobable parameters. Zw and Zµ are normalization fac(t)(t)
and βµi(F, m) are parameters determining
tors, and βwi
how much emphasis is put on the maximum value. The
prior distribution is not informative (i.e., it is uniform)(t)(t)
and βµi(F, m) are 0, corresponding to the casewhen βwiwhen no prior knowledge is available. In Equations (13)(t)(t)
; w(t) ) and Dµ (µ(t)
and (14), Dw (w0i
0i (F, m); µ (F, m))
are the following Kullback-Leibler information:
 Fhi
Miw (t) (F, m)(t)(t)
Dw (w0idF,(15); w(t) ) =w0i(F, m) log 0iw (t) (F, m)
Fli
(t)

m=1

(t)
Dµ (µ(t)0i (F, m); µ (F, m))

=

Hi

c(t)
0i (h|F, m)
.(16)
c(t)
(h|F,
m)
log
0i
c(t) (h|F, m)
h=12.2.3. MAP estimation using the EM algorithm
The problem to be solved is to estimate the model parameter θ(t) , taking into account the prior distribution p0i (θ(t) ),(t)iswhen we observe p(t)
Ψ (x). The MAP estimator of θ
obtainedbymaximizing
 ∞
(t)
(t)
p(t)
(17)
Ψ (x) (log p(x|θ ) + log p0i (θ )) dx.
−∞

Because this maximization problem is too difficult tosolve analytically, we use the Expectation-Maximization(EM) algorithm [9], which is an algorithm iteratively applying two steps — the expectation step (E-step) and the
maximization step (M-step) — to compute MAP estimatesfrom incomplete observed data (i.e., from p(t)
Ψ (x)). With
respect to θ(t) , each iteration updates the old estimate
θ(t) = {w(t) , µ(t) } to obtain the new (improved) estimate
θ(t) = {w (t) , µ(t) }. For each frame t, w(t) is initializedwith the final estimate w(t−1) after iterations at the previous frame t − 1; µ(t) is initialized with the most probableparameter µ(t)0i in the current implementation.
By introducing the hidden (unobservable) variables
F , m, and h, which, respectively, describe which F0,
which tone model, and which harmonic component were
responsible for generating each observed 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:
QMAP (θ(t) |θ (t) ) = Q(θ(t) |θ(t) ) + log p0i (θ(t) ),
(18)
 ∞
p(t)
Q(θ(t) |θ (t) ) =
Ψ (x)
−∞

EF,m,h [log p(x, F, m, h|θ(t) ) | x, θ(t) ] dx, (19)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 expec-

tation of a with respect to the hidden variables F , m,and h, with the probability distribution determined bycondition b.
2. (M-step)
Maximize QMAP (θ(t) |θ (t) ) as a function of θ (t) to obtain the updated (improved) estimate θ (t) :θ(t) = argmax QMAP (θ(t) |θ(t) ).
(20)
θ(t)
(t) (t)In the E-step, Q(θ |θ ) is expressed as ∞  Fhi Mi Hi(t) (t)p(t)
Q(θ |θ ) =Ψ (x)−∞ Flim=1 h=1p(F, m, h|x, θ(t) ) log p(x, F, m, h|θ(t) )dF dx, (21)where the complete-data log-likelihood is given bylog p(x, F, m, h|θ (t) )= log(w(t) (F, m) p(x, h|F, m, µ(t) (F, m))).

(22)

Regarding the M-step, Equation (20) is a conditional
problem of variation, where the conditions are given by
Equations (5) and (10). This problem can be solved by using Euler-Lagrange differential equations with Lagrange
multipliers [7, 8] and we obtain the following new parameter estimates:(t) (t)w (t) (F, m) + βwi
w0i (F, m)w(t) (F, m) = ML
,
(23)
(t)
1 + βwic(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),(24)(t)(F, m) and c(t)where wML
ML (h|F, m) are, when the nonin(t)(t)= 0 and βµi(F, m) = 0)
formative prior distribution (βwi
is given, the following
maximum
likelihood
estimates:

(t)
wML
(F, m) =

∞

−∞

p(t)
Ψ (x)

w (t) (F, m) p(x|F, m, µ(t) (F, m))dx, (25)
 Fhi Mi
(t) (η, ν) p(x|η, ν, µ(t) (F, ν)) dη
ν=1 w
Fli
 ∞
1
(t)
cML (h|F, m) =
p(t)
Ψ (x)
(t)
wML (F, m) −∞w(t) (F, m) p(x, h|F, m, µ(t) (F, m))dx. (26)
 Fhi Mi
(t) (η, ν) p(x|η, ν, µ(t) (F, ν)) dη
ν=1 w
Fli
For an intuitive explanation of Equation (25), we callw (t) (F,m) p(x|F,m,µ(t) (F,m)) Fhi the decomposition
Mi
(t)
(t)
Fli

ν=1

w

(η,ν) p(x|η,ν,µ(F,ν)) dηfilter. For the integrand on the right side of Equation (25),
we can consider that, because of this filter, the value ofp(t)Ψ (x) at frequency x is decomposed into (is distributed
among) all possible tone models p(x|F, m, µ(t) (F, m))(Fli ≤ F ≤ Fhi , 1 ≤ m ≤ Mi ) in proportion to the numerator of the decomposition filter at x. The higher theweight w (t) (F, m), the larger the decomposed value given
to the corresponding tone model. Note that the value ofp(t)Ψ (x) at different x is also decomposed according to adifferent ratio in proportion to the numerator of the decomposition filter at that x. Finally, the updated weightII - 1087

(t)wML(F, m) is obtained by integrating all the decomposedvalues given to the corresponding m-th tone model for the
F0 F .
We think that this decomposition behavior is the advantage of PreFEst in comparison to previous comb-filterbased or autocorrelation-based methods [1, 2, 3]. Thisis because these previous methods cannot easily support
the decomposition of an overlapping frequency component (overtone) shared by several simultaneous tones andtend to have difficulty distinguishing sounds with overlapping overtones. In addition, PreFEst can simultane(t)ously estimate all the weights wML(F, m) (for all the rangeof F ) so that these weights can be optimally balanced:
it does not determine the weight at F after determining
the weight at another F . We think this simultaneous estimation of all the weights is an advantage of PreFEst
compared to previous recursive-subtraction-based methods [1, 4] where components of the most dominant harmonic structure identified are subtracted from a mixture
and then this is recursively done again starting from theresidue of the previous subtraction. In these methods,
once inappropriate identification or subtraction occurs,
the following recursions starting from the wrong residue
become unreliable.
After the above iterative computation of Equations(23) and (24), the F0’s PDF p(t)F 0 (F ) can be obtained fromw (t) (F, m) according to Equation (11). We can also obtain the tone-model shape c(t) (h|F, m), which is the relative amplitude of each harmonic component of all typesof tone models p(x|F, m, µ(t) (F, m)).

and the bass lines provided by a bass guitar or a contrabass: the average detection rate was 88.4% for the melody
line and 79.9% for the bass line.

4. Conclusion
I have introduced the PreFEst method that estimates themost predominant F0 in a monaural sound mixture without assuming the number of sound sources. Although
PreFEst has great potential, I have not fully exploited
it. In the future, for example, many different tone models could be prepared by analyzing or learning various
kinds of harmonic structure that appear in music andmultiple peaks in the F0’s PDF, each corresponding toa different sound source, could be tracked simultaneously by using a sound source discrimination method.
While I dealt with only harmonic-structure tone models inthis paper, PreFEst can be applied to any weighted mixture of arbitrary tone models (even if their components
are inharmonic) by simply replacing Equation (2) withp(x, h|F, m, µ(t) (F, m)) = c(t) (h|F, m)parbit (x; F, h, m),where parbit (x; F, h, m) is an arbitrary PDF (h is merely
the component number in this case). Even with this general tone model, in theory the F0’s PDF can be estimatedby using the same Equations (23) and (24). Both any
harmonic- and inharmonic-structure tone models can also
be used together. Moreover, PreFEst can be applied to
non-music audio signals. In fact, Masuda-Katsuse [10]has extended it and demonstrated its effectiveness forspeech recognition in realistic noisy environments.
Acknowledgments: I thank Shotaro Akaho and Hideki
Asoh for their valuable discussions.

2.3. PreFEst-back-end: Sequential F0 tracking by
multiple-agent architecture
A simple way to identify the most predominant F0 is tofind the frequency that maximizes the F0’s PDF. This result is not always stable, 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 continuityof the F0 by using a multiple-agent architecture [5, 6, 8]
in which agents track different temporal trajectories of the
F0. The final F0 output is determined on the basis of themost dominant and stable F0 trajectory.

3. Experimental results
The PreFEst has been implemented in a real-time systemthat takes a musical audio signal as input and outputs thedetected melody and bass lines in several forms, such as
audio signals and computer graphics. The current implementation uses two adaptive tone models and main parameter values are described in [7, 8]. The system wastested on excerpts from 10 musical pieces in the popular,
jazz, and orchestral genres. The 20-s-long input monaural
audio signals — each containing a single-tone melody andthe sounds of several instruments — were sampled fromcompact discs. We evaluated the detection rates by comparing the estimated F0s with the correct F0s that were
hand-labeled using our F0 editor program [6, 8].
In our experiment the system correctly detected, formost parts of each audio sample, the melody lines provided by a voice or a single-tone mid-range instrument

5. References[1] Alain de Cheveigné, “Separation of concurrent harmonic sounds:
Fundamental frequency estimation and a time-domain cancellation model of auditory processing,” J. Acoust. Soc. Am.,
93(6):3271–3290, 1993.
[2] Alain de Cheveigné and Hideki Kawahara, “Multiple period estimation and pitch perception model,” Speech Communication,
27(3–4):175–185, 1999.
[3] Tero Tolonen and Matti Karjalainen, “A computationally efficient multipitch analysis model,” IEEE Trans. on Speech andAudio Processing, 8(6):708–716, 2000.
[4] Anssi P. Klapuri, “Multipitch estimation and sound separation by
the spectral smoothness principle,” In Proc. ICASSP 2001, 2001.
[5] Masataka 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 Computational Auditory
Scene Analysis, pp. 31–40, 1999.
[6] Masataka Goto, “A robust predominant-F0 estimation method
for real-time detection of melody and bass lines in CD recordings,” In Proc. ICASSP 2000, pp. II–757–760, 2000.
[7] Masataka Goto, “A predominant-F0 estimation method for CD
recordings: MAP estimation using EM algorithm for adaptive
tone models,” In Proc. ICASSP 2001, pp. V–3365–3368, 2001.
[8] Masataka Goto, “A real-time music scene description system:
Predominant-F0 estimation for detecting melody and bass lines
in real-world audio signals,” Speech Communication, 2004, (accepted).
[9] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. Roy. Stat.
Soc. B, 39(1):1–38, 1977.
[10] Ikuyo Masuda-Katsuse, “A new method for speech recognitionin the presence of non-stationary, unpredictable and high-level
noise,” In Proc. Eurospeech 2001, pp. 1119–1122, 2001.

II - 1088