11th International Society for Music Information Retrieval Conference (ISMIR 2010)

INFINITE LATENT HARMONIC ALLOCATION: A NONPARAMETRICBAYESIAN APPROACH TO MULTIPITCH ANALYSIS
Kazuyoshi Yoshii
Masataka Goto
National Institute of Advanced Industrial Science and Technology (AIST), Japan{k.yoshii, m.goto}@aist.go.jp

ABSTRACT

Our method (iLHA)

Conventional methods
Assume the number of spectral bases
Assume the number of harmonic partials Input:
Carefully tune prior distributions
PolyphonicThis paper presents a statistical method called Infinite Latent Harmonic Allocation (iLHA) for detecting multiple
fundamental frequencies in polyphonic audio signals. Conventional methods face a crucial problem known as model
selection because they assume that the observed spectra aresuperpositions of a certain fixed number of bases (soundsources and/or finer parts). iLHA avoids this problem byassuming that the observed spectra are superpositions ofa stochastically-distributed unbounded (theoretically infinite) number of bases. Such uncertainty can be treated in a
principled way by leveraging the state-of-the-art paradigmof machine-learning called Bayesian nonparametrics. To
represent a set of time-sliced spectral strips, we formulated
nested infinite Gaussian mixture models (GMMs) basedon hierarchical and generalized Dirichlet processes. Each
strip is allowed to contain an unbounded number of soundsources (GMMs), each of which is allowed to contain an
unbounded number of harmonic partials (Gaussians). To
train the nested infinite GMMs efficiently, we used a modern inference technique called collapsed variational Bayes(CVB). Our experiments using audio recordings of real piano and guitar performances showed that fully automatediLHA based on noninformative priors performed as well as
optimally tuned conventional methods.

audio signals
Model #1

Model #2

Select an optimal model…

RequireNo assumptions
No prior tuning
UnifiedModel

Output:Fundamental frequencies (F0s)

Figure 1. Methodological advantage of our method.
A crucial problem in these methods, known as model
selection, is that they perform best only if an appropriate
model complexity (the number of bases) is specified in advance. One might think that the optimal number of bases
must be equal to the number of sound sources, but it is notclear how many bases are most suited to represent a single
source if the spectral shape varies through time. Although
uncertainty is inherent in model selection, conventional
methods assume that a certain complexity exists uniquely
as an oracle. As shown in Figure 1, they require possiblemodels to be examined separately and exhaustively and the
optimal model selected in retrospect. Such a deterministic framework is not easy-to-use in practice although optimally tuned methods can achieve good performance.
To avoid model selection, we propose a novel statisticalmethod called Infinite Latent Harmonic Allocation (iLHA)
based on a modern paradigm of machine learning called
Bayesian nonparametrics. Note that the term “nonparametric” means that we do not have to fix model complexity
uniquely. We assume that an unbounded but finite numberof bases stochastically appears in a limited amount of available data although an infinite number of bases theoreticallyexists in the universe. Uncertainty in model selection canbe treated reasonably in a probabilistic framework.
iLHA can be derived by taking the infinite limit of conventional finite models [3, 4]. Conventionally, each spectral basis is often parameterized by means of a Gaussian
mixture model (GMM) in which a fixed number of Gaussians corresponds to the spectral peaks of harmonic partials, and a time-sliced polyphonic spectral strip is modeledby mixing a fixed number of GMMs. Here, we consider
both the number of bases and the number of partials to approach infinity, where most are regarded as unnecessary
and automatically removed through statistical inference.
A fundamental and practically-important advantage ofiLHA is that precise prior knowledge is not required. Conventional methods [3–5] heavily rely on prior distributions
regarding the relative strengths of harmonic partials, which

1. INTRODUCTION
Multipitch analysis of polyphonic audio signals [1–11] is
one of the most important issues because it is the basisof many applications such as music transcription, chord
recognition, and musical instrument recognition. We focus
on principled methods based on machine learning, whichhave recently yielded promising results. Some researchers,
for example, have proposed generative probabilistic models that explain how multiple spectral/signal bases (compositional units) are mixed to form polyphonic music [3–6].
The model parameters can be trained by means of statistical inference. Others have used nonnegative matrix factorization (NMF) to decompose polyphonic spectra into individual spectral bases [7–11]. NMF can be interpreted fromthe viewpoint of statistical inference [10–12].
Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that copies
bear this notice and the full citation on the first page.
c 2010 International Society for Music Information Retrieval.
309

11th International Society for Music Information Retrieval Conference (ISMIR 2010)
Density

have too much impact on performance, and forced us totune priors and their weighting factors by hand according
to the properties of target sound sources. iLHA, in contrast,
can be fully automated by layering noninformative hyperpriors on influential priors in a hierarchical Bayesian manner. This is consistent with the fact that humans can adaptively distinguish individual notes of various instruments.
One of major contributions of this study is to embody thefundamental Bayesian principle “Let the data speak for itself” in the context of multipitch analysis.
The rest of this paper is organized as follows: Section 2
describes statistical interpretation of polyphonic spectra.
Section 3 discusses related work. Sections 4 and 5 explain
finite models (LHA) and infinite models (iLHA). Section 6
reports our experiments. Section 7 concludes this paper.

Λ−k1

τ k1

τk2

µk

τ km

µ k + o2 L μk + om L

Frequency(Cent)

Density

Figure 2. Probabilistic model of a single basis.

Frequency(Cent)

Mix two GMMs

Density

Assumption: Amplitude additivity
in the frequency domain is satisfied
when multiple sounds are mixed

2. STATISTICAL INTERPRETATION
We interpret polyphonic spectra as histograms of observed
frequencies that independently occur. This interpretation
basically follows conventional studies [3–5].

This density distribution corresponds toobserved frequency histogram (spectral strip)

2.1 Assumptions

Frequency(Cent)

Figure 3. Probabilistic model of mixed multiple bases.

Suppose given polyphonic audio signals are generated fromK bases, each of which consists of M harmonic partials
located on a linear frequency scale at integral multiples ofthe fundamental frequency (F0). Note that each basis canbe associated with multiple sounds of different temporal
positions if these sounds are derived from the same pitchof the same instrument. We transform the audio signals
into wavelet spectra. Let D be the number of frames. If aspectral strip at frame d (1 ≤ d ≤ D) has amplitude a atfrequency f , we assume that frequency f was observed atimes in frame d. Assuming that amplitudes are additive,
we can consider each observed frequency to be generated
from one of M partials in one of K bases.
These notations are for the finite case. In Bayesian nonparametrics, we take the limit as K and M go to infinity.

3.1 Statistical Inference
Statistical methods [3–6] assume probabilistic models using a limited number of parameters to represent the generative process of observed spectra (audio signals). F0 estimation directly corresponds to finding model parametersthat provide the best explanations of the given data.
Goto [3] first proposed probabilistic models of harmonicsounds by regarding frequency spectra as probabilistic densities (histograms of observed frequencies).
As shown in Figure 2, the spectral distribution of basisk (1 ≤ k ≤ K) is modeled by a harmonic GMM:
Mk (x) =M


m=1

2.2 Observed and Latent Variables ,
τkm N xμk + om , Λ−1k(1)where x is a one-dimensional vector that indicates an observed frequency [cents]. 1 The Gaussian parameters, mean
μk and precision Λk , indicate F0 [cents] of basis k and a
degree of energy concentration around the F0 in the frequency domain. τkm is a relative strength of the m-th harmonic partial (1 ≤ m ≤ M ) in basis k. We set om to[1200 log2 m]. This means M Gaussians are located to have
harmonic relationships on the logarithmic frequency scale.
As shown in Figure 3, the spectral strip of frame d ismodeled by mixing K harmonic GMMs as follows:

Let the total observed variables over all D frames be represented by X = {X 1 , · · · , X D }, where X d is a set ofobserved frequencies X d = {xd1 , · · · , xdNd } in frame d.
Nd is the number of frequency observations. That is, Nd isequal to the sum of spectral amplitudes over all frequency
bins in frame d. xdn (1 ≤ n ≤ Nd ) is an one-dimensional
vector that represents an observed frequency.
Let the total latent variables corresponding to X be similarly represented by Z = {Z 1 , · · · , Z D }, where Z d ={z d1 , · · · , z dNd }. z dn is a KM -dimensional vector in which
only one entry, zdnkm , takes a value of 1 and the others takevalues of 0 when frequency xdn is generated from partialm (1 ≤ m ≤ M ) of basis k (1 ≤ k ≤ K).

Md (x) =Kπdk Mk (x)

(2)

k=1

where πdk is a relative strength of basis k in frame d. Therefore, the polyphonic spectral strip is represented by nested
finite Gaussian mixture models.
Several inference methods that have been proposed forparameter estimation are listed in Table 1. Goto [3] pro-

3. COMPARISON WITH RELATED WORK
The properties of iLHA are intermediate between those of
two successful approaches–statistical inference and NMF–
which are discussed here for comparison and to clarify the
positioning of our approach.

1 Linear frequency f in hertz can be converted to logarithmic freh
3
quency fc in cents as fc = 1200 log2 (fh /(440 12 −5 )).

310

11th International Society for Music Information Retrieval Conference (ISMIR 2010)
PreFEst [3]HC [4]HTC [5]
NMF [7]iLHA#(bases)
Fixed
Inferred
Fixed
Fixed
Infinite#(partials)
Fixed
Fixed
Fixed
Not used
InfiniteTemporal modeling
None
None
Continuity treatedExchangeable
Exchangeable

πd

α

ν

β

υ

Table 1. Comparison of multipitch analysis methods.

p(π) =p(τ ) =
p(μ, Λ) =

µk

Λk K

m0 ,b0W0 ,c0Dir(τ k |βυ) ∝

K


K
D 


πdk k

(6)

M
K 


βυm −1
τkm

(7)

αν −1

d=1 k=1

k=1 m=1

  

 
N μk m0 , (b0 Λk )−1 W Λk W 0 , c0 (8)where p(π) and p(τ ) are products of Dirichlet distributions
and p(μ, Λ) is a product of Gaussian-Wishart distributions.
αν and βυ are hyperparameters and α and β are called
concentration parameters when ν and υ sum to unity. m0 ,b0 , W 0 , and c0 are also hyperparameters; W 0 is a scale
matrix and c0 is a degree of freedom.
4.2 Variational Bayesian Inference
The objective of Bayesian inference is to compute a true
posterior distribution of all variables: p(Z, π, τ , μ, Λ|X).
Because analytical calculation of the posterior distribution
is intractable, we instead approximate it by using iterative
inference techniques such as variational Bayes (VB) andMarkov chain Monte Carlo (MCMC). Although MCMC
is considered to be more accurate in general, we use VB
because it converges much faster.
In the VB framework, we introduce a variational posterior distribution q(Z, π, τ , μ, Λ) and make it close to thetrue posterior p(Z, π, τ , μ, Λ|X) iteratively. Here, we assume that the variational distribution can be factorized asq(Z, π, τ , μ, Λ) = q(Z)q(π, τ , μ, Λ)(9)To optimize q(Z, π, τ , μ, Λ), we use a variational version of the Expectation-Maximization (EM) algorithm [13].
We iterate VB-E and VB-M steps until a variational lower
bound of evidence p(X) converges as follows:q ∗ (Z) ∝ exp Eπ ,τ ,μ,Λ [log p(X, Z, π, τ , μ, Λ)] (10)4.1 Model Formulationq ∗ (π, τ , μ, Λ) ∝ exp (EZ [log p(X, Z, π, τ , μ, Λ)]) (11)Figure 4 illustrates a graphical representation of the LHA
model. The full joint distribution is given by

4.3 Updating Formula
We derive the formulas for updating variational posterior
distributions according to Eqns. (10) and (11).

p(X, Z, π, τ , μ, Λ)= p(X|Z, μ, Λ)p(Z|π, τ )p(π)p(τ )p(μ, Λ)(3)4.3.1 VB-E Step

where the first two terms on the right-hand side are likelihood functions and the other three terms are prior distributions. The likelihood functions are defined as

An optimal variational posterior distribution of latent variables Z can be computed as follows:
log q ∗ (Z) = Eπ,τ ,μ,Λ [log p(X, Z, π, τ , μ, Λ)] + const.

zdnkm
N xdn μk + om , Λ−1(4)k(πdk τkm )

K


k=1

This section explains LHA, the finite version of iLHA, as apreliminary step to deriving iLHA. We formulate the conventional nested finite GMMs in a Bayesian manner.

p(Z|π, τ ) =Dir(πd |αν) ∝k=1

4. LATENT HARMONIC ALLOCATION



D


d=1

NMF-based methods [7–12] factorize observed frequency
spectra into the product of spectral bases and time-varying
envelopes under the nonnegativity constraint. K bases are
estimated by sweeping all frames of the given spectra. Although several methods [10, 12] take temporal continuity
into account, standard methods are based on temporal exchangeability. In other words, exchange of arbitrary frames
does not affect the factorized results. Although such temporal modeling is not sufficient, it is known to work wellin practice. Therefore, iLHA adopted the exchangeability.

dnkm

M

D

Then, we introduce conjugate priors as follows:

3.2 Nonnegative Matrix Factorization



τkm



Figure 4. A graphical representation of LHA.

posed a method called PreFEst that estimates only relative
strengths τ and π while μ and Λ are fixed by allocating
many GMMs to cover the entire frequency range as F0 candidates. Kameoka et al. [4] then proposed harmonic clustering (HC), which estimates all the parameters and selects
the optimal number of bases by using the Akaike information criterion (AIC). Although these methods yielded thepromising results, they analyze the spectral strips of different frames independently. Thus, Kameoka et al. [5] proposed harmonic-temporal-structured clustering (HTC) that
captures temporal continuity of spectral bases. Note that
all these methods are based on maximum-likelihood andmaximum-a-posteriori training of the parameters by introducing prior distributions of relative strengths τ , whichhave a strong impact on the accuracy of F0 estimation.
Our method called iLHA is based on hierarchical nonparametric Bayesian modeling that requires no prior tuning
and avoids specifying K and M in advance. More specifically, the limit of the conventional nested finite GMMs is
considered as K and M diverge to infinity.

p(X|Z, μ, Λ) =

xdn

zdn

zdnkm

= Eμ,Λ [log p(X|Z, μ, Λ)] + Eπ,τ [log p(Z|π, τ )] + const.
=
zdnkm log ρdnkm + const.
(12)

(5)

dnkm

dnkm

311

11th International Society for Music Information Retrieval Conference (ISMIR 2010)

ν1
ν~1

:

Recursively breakthe stick of length 1

L

ν2

ν3

ν~2

1−ν~1
:
1−ν~2
ν~ : 1−ν~

5.1 Model Formulation
First we let K approach infinity, where the infinite numberof harmonic GMMs is assumed to exist in the universe.
More specifically, the dimensionality of the Dirichlet distributions in Eqn. (6) is considered to be infinite. At eachframe d, π d is an infinite vector of normalized probabilities (mixing weights) drawn from the infinite-dimensional
Dirichlet prior. Such stochastic process is called a Dirichlet
process (DP). Every time frequency xdn is generated, oneof the infinite number of harmonic GMMs is drawn according to πd . Note that most entries of πd take extremely tiny
values because all entries sum to unity. If we can observe
the infinite number of frequencies (Nd → ∞), the infinite
number of harmonic GMMs can be drawn. However, Nd
is finite in practice. Therefore, only the finite number of
harmonic GMMs, K+  ∞, is drawn at frame d. Here,a problem is that harmonic GMMs that are actually drawn
at frame d are completely disjointed from those drawn atanother frame d . This is not a reasonable situation.
To solve this problem, we use the hierarchical Dirichlet
Process (HDP) [14]. More specifically, we assume that
infinite-dimensional hyperparameter ν in Eqn. (6), whichis shared among all D frames, is a draw from a top-level
DP. A generative interpretation is that after an unbounded
number of harmonic GMMs is initially drawn from the toplevel DP, an unbounded subset is further drawn according
to the local DP at each frame. This effectively ties frame dto another frame d . As shown in Figure 5, ν is known tofollow the stick-breaking construction [14] as follows:3

3

Figure 5. Stick-breaking construction of Dirichlet process.
where ρdnkm is defined aslog ρdnkm = Eπd [log πdk ] + Eτ k [log τkm ]+Eμk ,Λk log N xdn μk + om , Λ−1(13)kq ∗ (Z) is obtained as multinomial distributions given by zdnkm
q ∗ (Z) =
γdnkm
(14)
dnkm

where γdnkm is given by γdnkm =  ρdnkm
and is calledkm ρdnkma responsibility that indicates how likely it is that observed
frequency xdn is generated from harmonic partial m of basis k. Here, let ndkm be an observation count that indicates
how many frequencies were generated from harmonic partial m of basis k in frame d. ndkm and its expected value
can be calculated as follows:
ndkm =

E[ndkm ] =

zdnkm

n



γdnkm

(15)

n

For convenience in executing the VB-M step, we compute several sufficient statistics as follows:
Sk [1] ≡



γdnkm



γdnkm xdnm xTdnm

Sk [x] ≡dnm

Sk [xxT ] ≡γdnkm xdnm (16)

dnm

(17)

νk = ν̃kdnm

Consequently, an optimal variational posterior distributionof parameters π, τ , μ, Λ is shown to be given byq ∗ (π d )

d=1

K


k=1

q ∗ (τ k )Kq ∗ (μk , Λk ) (18)k=1

τkm = τ̃km

Since we use conjugate priors, each posterior has the same
form of the corresponding prior as follows:q ∗ (π d ) = Dir(π d |αd )q ∗ (τ k ) = Dir(τ k |β k )
(20)
 
  

q (μk , Λk ) = N μk mk , (bk Λk )−1 W Λk W k , ck (21)where the variational parameters are given by
bk = b0 + Sk [1]ck = c0 + Sk [1]b0 m0 + Sk [x]b0 m0 + Sk [x]=mk =b0 + Sk [1]
bk

−1
T
T
T
W −1
k = W 0 + b0 m0 m0 + Sk [xx ] − bk mk mkm−1(1 − τ̃km ) τ̃km ∼ Beta(βλ1 , βλ2 )(27)where β is a positive scalar and λ1 + λ2 = 1.
Because α, β, γ and λ are influential hyperparameters,
we put Gamma and Beta hyperpriors on them as follows:(19)

βkm = βυm + E[n·km ]

(26)

m =1

∗

αdk = ανk + E[ndk· ]ν̃k ∼ Beta(1, γ)

where γ is a concentration parameter of the top-level DP.
Therefore ν can be converted into ν̃ .
Now we let M approach infinity, where each harmonic
GMM consists of the infinite number of harmonic partials.
To put effective priors on τ , we use generalized DPs called
Beta two-parameter processes as follows:

4.3.2 VB-M StepD

(1 − ν̃k )k =1where xdnm is defined as xdnm = xdn − om .

q ∗ (π, τ , μ, Λ) =

k−1


(22)

p(α) = Gam(α|aα , bα ) p(γ) = Gam(γ|aγ , bγ )(28)p(β) = Gam(β|aβ , bβ ) p(λ) = Beta(λ|u1 , u2 )(29)

where a{α,β,γ} and b{α,β,γ} are shape and rate parameters.
Figure 6 shows a graphical representation of the iLHA
model. The full joint distribution is given by

(23)
(24)

p(X, Z, π, τ̃ , μ, Λ, α, β, γ, λ, ν̃) = p(X|Z, μ, Λ)p(μ, Λ)(25)p(Z|π, τ̃ )p(π|α, ν̃)p(τ̃ |β, λ)p(α)p(β)p(γ)p(λ)p(ν̃|γ) (30)Here dot (·) denotes the sum over that index.

where p(Z|π, τ̃ ) is obtained by plugging Eqn. (27) into
Eqn. (5) and p(π|α, ν) is the same as Eqn. (6). p(ν̃|γ) andp(τ̃ |β, λ) are defined according to Eqns. (26) and (27) as

5. INFINITE LATENT HARMONIC ALLOCATIONThis section derives hierarchical nonparametric Bayesianmodels, i.e., nested infinite GMMs for polyphonic spectra.

p(ν̃|γ) =


k

312

Beta(τ̃k |1, γ) p(τ̃ |β, λ) =km

Beta(τ̃km |βλ)(31)11th International Society for Music Information Retrieval Conference (ISMIR 2010)

xdn

zdn

πd

α

aα , bαaγ , bγ

γ

ν~k

τ~km

K →∞
M →∞

aβ , bβ

β

λ



µk

Λk

u m0 ,b0

W0 ,c0

M

sd

ηd

α

aα , bα

xdn

zdn

D
K

ν~k

γ

aγ , bγ

K →∞M → ∞ aβ , bβ



ξkm
tkm M

β

u m0 ,b0λDK

W0 ,c0Figure 7. A collapsed model with auxiliary variables.

Figure 6. A graphical representation of iLHA.
5.2 Collapsed Variational Bayesian Inference

5.3.1 CVB-E Step

To train the HDP model we use a sophisticated version of
VB called collapsed variational Bayes (CVB) [15]. CVB
enables more accurate posterior approximation in the spaceof latent variables where parameters are integrated out.
Figure 7 shows a collapsed iLHA model. By integratingout π, τ̃ , μ, Λ, we obtain the marginal distribution given by

A variational probability of zdnkm = 1 is given by
log q ∗ (zdnkm = 1) = Ez¬dn log G[ανk ] + n¬dn
dk·
 

m−1
 G[βλ2 ] + n¬dnG[βλ1 ] + n¬dn
·k>m
·km
+ Ez¬dn log¬dnE[β] + n¬dn·k≥m m =1 E[β] + n·k≥m¬dn ¬dn(37)+ Ez¬dn log S(xdnm |m¬dn
zk , Lzk , czk ) + const.

p(X, Z, α, β, γ, λ, ν̃)= p(X|Z)p(Z|α, β, λ, ν̃)p(α)p(β)p(γ)p(λ)p(ν̃|γ)

where subscript ¬dn denotes a set of indices without d andn, G[x] denotes the geometric average exp(E[log x]), andS is the Student-t distribution. L¬dn
is given by L¬dn
=
zk
zk

(32)

where the first two terms are calculated as follows:1n···p(X|Z) = (2π)− 2p(Z|α, β, λ, ν̃) =db0 2 B(W 0 , c0 )bzk
B(W zk , czk )kΓ(α)Γ(ανk + ndk· )Γ(α + nd·· )k¬dn¬dn ¬dn, where m¬dn
are obzk , bzk , W zk , czk
tained by substituting zdnkm for γdnkm required by Eqns.
(23), (24), and (25) and calculating sum without z dn . Eachterm of Eqn. (37) can be calculated efficiently [15, 16].

Γ(ανk ) Γ(β)Γ(βλ1 + n·km )Γ(βλ2 + n·k>m )
km

b¬dn
¬dn
zk
c¬dn
zk W zk
1+b¬dnzk(33)

Γ(βλ1 )Γ(βλ2 )Γ(β + n·k≥m )5.3.2 CVB-M Step

(34)

First, α, β and γ are Gamma distributed as follows:

where bzk , W zk , czk are obtained by substituting zdnkm for
γdnkm in calculating Eqns. (23) and (25).
Because CVB cannot be applied directly to Eqn. (32),
we introduce auxiliary variables by using a technique called
data augmentation [15]. Let ηd and ξkm be Beta-distributed
variables and sdk and tkm be positive integers that satisfy
1 ≤ sdk ≤ ndk· , 1 ≤ tkm1 ≤ n·km , and 1 ≤ tkm2 ≤ n·k>m .
Eqn. (34) can be augmented asp(Z, η, ξ, s, t|α, β, λ, ν̃) =q(α) ∝ αaα +E[s·· ]−1 e−α(bα − d E[log ηd ])
q(β) ∝ β aβ +E[t··· ]−1 e−β (bβ − km E[log ξkm ])q(γ) ∝ γkq∗u +E[t··1 ]−1 u2 +E[t··2 ]−1
λ2
1+E[s·k ]−1(ν̃k ) ∝ ν̃k(1 − ν̃k )E[γ]+E[s·>k ]−1where [ ] denotes a Stirling number of the first kind. The
augmented marginal distribution is given byE[α]−1q ∗ (ηd ) ∝ ηdq∗∗

(1 − ηd )nd·· −1

E[β]−1(ξkm |Z) ∝ ξkm (1 − ξkm )n·k≥m −1

ndk·
s



G[ανk ]
nq ∗ (tkm1 = t|Z) ∝ ·km
G[βλ1 ]ttnq ∗ (tkm2 = t|Z) ∝ ·k>mG[βλ2 ]ttq (sdk = s|Z) ∝p(X, Z, η, ξ, s, t, α, β, γ, λ, ν̃)= p(Z, η, ξ, s, t|α, β, λ, ν̃)p(α)p(β)p(γ)p(λ)p(ν̃|γ) (35)In the CVB framework, we assume that the variational
posterior distribution can be factorized as follows:

(40)

(41)
(42)

s

(43)
(44)
(45)
(46)
(47)

To calculate E[sdk ] (average sdk over Z ), we exactly treat
the case ndk· = 0 and apply second-order approximation
when ndk· > 0 (see details in [15]). E[log ξkm ], E[tkm1 ],and E[tkm2 ] can be calculated in the same way.
To estimate F0s, we need explicitly compute the variational posteriors of the integrated-out parameters μ, Λ. To
do this, we execute the standard VB-M step once by using
the responsibilities q(Z) obtained in the CVB-E step.

q(Z, η, ξ, s, t, α, β, γ, λ, ν̃)q(z dn )(39)

Finally, the variational posteriors of η, ξ, s, t are given bykm

k E[log(1−ν̃k )])q ∗ (λ) ∝ λ1 1 ξ β−1 (1 − ξkm )n·k≥m −1  n
n·k>m
·km
km
(βλ1 )tkm1(βλ2 )tkm2
tkm1
tkm2
Γ(n·k≥m )= q(α, β, γ, λ)q(ν̃)q(η, ξ, s, t|Z)e

Then, λ and τ̃ are Beta distributed as follows: η α−1 (1 − ηd )nd·· −1   ndk· d(ανk )sdk
sdk
Γ(nd·· )daγ +K−1 −γ (bγ −

(38)

(36)

dn

We also use an approximation technique called variational posterior truncation. More specifically, we assumeq(zdnkm ) = 0 when k > K and m > M . In practice, it isenough that K and M are set to sufficiently large integers.

6. EVALUATION

5.3 Updating Formula
We descrive the formulas for updating variational posterior
distributions.

This section reports our comparative experiments evaluating the performance of iLHA.

313

11th International Society for Music Information Retrieval Conference (ISMIR 2010)
Piece number
RWC-MDBJ-2001 No.1
J-2001 No.2
J-2001 No.6
J-2001 No.7
J-2001 No.8
J-2001 No.9
C-2001 No.30
C-2001 No.35
Total

Optimally tuned
PreFEst [3] HTC [5]
75.8
79.0
78.5
78.0
70.4
78.3
83.0
86.0
85.7
84.4
85.9
89.5
76.0
83.6
72.8
76.0
79.4
82.0

Fully automated
LHA
iLHA
70.7
82.2
69.1
77.9
49.8
71.2
70.2
85.5
55.9
84.6
68.9
84.7
81.4
81.6
58.9
79.6
65.8
81.7

Table 2. Frame-level F-measures of F0 detection.
6.1 Experimental Conditions
We evaluated LHA and iLHA on the same test set used in[5], which consisted of nine pieces of piano and guitar solo
performances excerpted from the RWC music database [17].
The first 23 [s] of each piece were used for evaluation. The
spectral analysis was conducted by the wavelet transformusing Gabor wavelets with a time resolution of 16 [ms].
The values and temporal positions of actual F0s were prepared by hand as ground truth. We evaluated performancein terms of frame-level F-measures. The priors and hyperpriors of LHA and iLHA were set to noninformative uniform distributions. K and M were set to sufficiently large
numbers, 60 and 15. iLHA is not sensitive to these values.
No other tuning was required. To output F0s at each frame,
we extracted bases whose expected weights π were over a
threshold, which was optimized as in [5].
For comparison, we referred to the experimental resultsof PreFEst and HTC reported in [5]. Although the groundtruth data was slightly different from ours, it would be sufficient for roughly evaluating performance comparatively.
The number of bases, priors, and weighting factors were
carefully tuned by using the ground-truth data to optimize
the results. Although this is not realistic, the upper boundsof potential performance were investigated in [5].
6.2 Experimental Results
The results listed in Table 2 show that the performance ofiLHA approached and sometimes surpassed that of HTC.
This is consistent with the empirical findings of many studies on Bayesian nonparametrics that nonparametric models
were competitive against optimally-tuned parametric models. HTC outperformed PreFEst because HTC can appropriately deal with temporal continuity of spectral bases.
This implies that incorporating temporal modeling would
improve the performance of iLHA.
The results of LHA were worse than those of iLHA because LHA is not based on hierarchical Bayesian modelingand requires precise priors. In fact, we confirmed that theresults of PreFEst and HTC based on MAP estimation weredrastically degraded when we used noninformative priors.
In contrast, iLHA stably showed the good performance.
7. CONCLUSION
This paper presented a novel statistical method for detecting multiple F0s in polyphonic audio signals. The method
allows polyphonic spectra to contain an unbounded number of spectral bases, each of which can consist of an un-

314bounded number of harmonic partials. These numbers canbe statistically inferred at the same time that F0s are estimated. Even in experimental evaluation using noninformative priors, our automated method performed well or better
than conventional methods manually optimized by trial anderror. To our knowledge, this is the first attempt to apply
Bayesian nonparametrics to multipitch analysis.
Bayesian nonparametrics is an ultimate methodologicalframework avoiding the model selection problem faced invarious areas of MIR. For example, how many sections
should one use for structuring a musical piece? How manygroups should one use for clustering listeners according totheir tastes or musical pieces according to their contents?
We are freed from these problems by assuming that in theory there is an infinite number of classes behind observed
data. Unnecessary classes are automatically removed fromconsideration through statistical inference. We plan to use
this powerful framework in a wide range of applications.
Acknowledgement: This study was partially supported by
JST CREST and JSPS KAKENHI 20800084.

8. REFERENCES[1] A. Klapuri. Multipitch analysis of polyphonic music andspeech signals using an auditory model. IEEE Trans. on
ASLP, Vol. 16, No. 2, pp. 255–266, 2008.
[2] M. Marolt. A connectionist approach to transcription of polyphonic piano music. IEEE Trans. on Multimedia, Vol. 6,
No. 3, pp. 439–449, 2004.
[3] M. Goto. A real-time music scene description system:
Predominant-F0 estimation for detecting melody and basslines in real-world audio signals. Speech Communication,
Vol. 43, No. 4, pp. 311–329, 2004.
[4] H. Kameoka et al. Separation of harmonic structures basedon tied Gaussian mixture model and information criterion for
concurrent sounds. ICASSP, Vol. 4, pp. 297–300, 2004.
[5] H. Kameoka et al. A multipitch analyzer based on harmonic temporal structured clustering. IEEE Trans. on ASLP,
Vol. 15, No. 3, pp. 982–994, 2007.
[6] A. Cemgil et al. A generative model for music transcription.
IEEE Trans. on ASLP, Vol. 14, No. 2, pp. 679–694, 2006.
[7] P. Smaragdis and J. Brown. Nonnegative matrix factorization
for polyphonic music transcription. WASPAA, 2003.
[8] A. Cont. Realtime multiple pitch observation using sparse
non-negative constraints. ISMIR, pp. 206–211, 2006.
[9] E. Vincent et al. Adaptive harmonic spectral decomposition
for multiple pitch estimation. IEEE Trans. on ASLP, Vol. 18,
No. 3, pp. 528–537, 2010.
[10] N. Bertin et al. Enforcing harmonicity and smoothness inBayesian non-negative matrix factorization applied to polyphonic music transcription. IEEE Trans. on ASLP, Vol. 18,
No. 3, pp. 538–549, 2010.
[11] P. Peeling et al. Generative spectrogram factorization models for polyphonic piano transcription. IEEE Trans. on ASLP,
Vol. 18, No. 3, pp. 519–527, 2010.
[12] T. Virtanen et al. Bayesian extensions to nonnegative matrix
factorisation for audio signal modelling. ICASSP, 2008.
[13] H. Attias. A variational Bayesian framework for graphicalmodels. NIPS, pp. 209–215, 2000.
[14] Y. W. Teh et al. Hierarchical Dirichlet processes. J. of Am.
Stat. Assoc., Vol. 101, No. 476, pp. 1566–1581, 2006.
[15] Y. W. Teh et al. Collapsed variational inference for HDP.
NIPS, Vol. 20, 2008.
[16] J. Sung et al. Latent-space variational Bayes. IEEE Trans. on
PAMI, Vol. 30, No. 12, pp. 2236–2242, 2008.
[17] M. Goto et al. RWC music database: Popular, classical, andjazz music database. ISMIR, pp. 287–288, 2002.