IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

717

A Nonparametric Bayesian Multipitch Analyzer
Based on Infinite Latent Harmonic Allocation
Kazuyoshi Yoshii, Member, IEEE, and Masataka Goto

Abstract—The statistical multipitch analyzer described in thispaper estimates multiple fundamental frequencies (F0s) in polyphonic music audio signals produced by pitched instruments.
It is based on hierarchical nonparametric Bayesian models thatcan deal with uncertainty of unknown random variables such asmodel complexities (e.g., the number of F0s and the number of
harmonic partials), model parameters (e.g., the values of F0s andthe relative weights of harmonic partials), and hyperparameters(i.e., prior knowledge on complexities and parameters). Using
these models, we propose a statistical method called infinite latentharmonic allocation (iLHA). To avoid model-complexity control,
we allow the observed spectra to contain an unbounded numberof sound sources (F0s), each of which is allowed to contain an
unbounded number of harmonic partials. More specifically, tomodel a set of time-sliced spectra, we formulated nested infinite
Gaussian mixture models based on hierarchical and generalized
Dirichlet processes. To avoid manual tuning of influential hyperparameters, we put noninformative hyperprior distributions onthem in a hierarchical manner. For efficient Bayesian inference,
we used a modern technique called collapsed variational Bayes.
In comparative experiments using audio recordings of piano andguitar solo performances, iLHA yielded promising results andwe found that there would be room for improvement based onmodeling of temporal continuity and spectral smoothness.
Index Terms—Bayesian nonparametrics, Dirichlet process, infinite latent harmonic allocation (iLHA), multipitch analysis.

I. INTRODUCTION
NCERTAINTY is inherent in music analysis. A musicalpiece about which we have little prior knowledge canoften be interpreted in various ways. One might, for example,
have various degrees of belief in different possible interpretations of tempo and semantic structures, and when we try to transcribe the music we hear in an audio recording, we often finddifficult to identify the notes with absolute confidence. Even ifin the end we need to determine which interpretation or transcription is the most reasonable, during the analysis it is important to keep all possibilities open with various degrees of belief.
We should therefore take an approach that can evaluate, propagate, and integrate the uncertainties of interdependent musicalelements or musical notes.

U

Manuscript received November 19, 2010; revised March 29, 2011; accepted
July 26, 2011. Date of publication August 15, 2011; date of current version January 11, 2012. This work was supported in part by CREST, JST, and KAKENHI
20800084. The associate editor coordinating the review of this manuscript andapproving it for publication was Dr. Daniel Ellis.
The authors are with the National Institute of Advanced Industrial Scienceand Technology (AIST), Tsukuba 305-8568, Japan (e-mail: k.yoshii@aist.go.jp;
m.goto@aist.go.jp).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TASL.2011.2164530

A natural way to manage uncertainty is to take a Bayesian
approach and use Bayesian probabilities to indicate degrees ofbelief. For example, suppose we have a distorted die. If the probabilities of getting the numbers 1, 2, 3, , 6 (called parameters)
are known, we can evaluate the likelihood for a set of numbers(called observed data) obtained by casting the die many times.
Note that the true values of the parameters do not vary stochastically. When the parameters are unknown, a probabilistic distribution is used as a means of representing how strongly possible
values are believed to be the true values. Such degrees of belief vary according to the amount of observed data. Before we
get observed data, prior distributions tend to be widely spread.
The more data we get, the sharper the peaks of posterior distributions become. That is, the degree of belief on a certain possibility increases. The objective of Bayesian inference is to calculate posterior distributions of unknown variables by formulating
probabilistic models defined by likelihood functions and priordistributions.
A critical problem in the conventional Bayesian approachis that we have to specify the complexity of the probabilistic
models in advance (complexity means the number of mixturesin Gaussian mixture models (GMMs) and the number of statesin hidden Markov models (HMMs)). If model complexities are
unknown, both the uncertainty of model complexities and that ofmodel parameters should be dealt with appropriately. The conventional approach, however, forces us to train many models
of different complexities independently and then select one according to some criteria. Such fine-comb model selection, ormodel-complexity control, is often impractical, especially in the
optimization of combinatorial-complexity models.
A nonparametric Bayesian approach avoiding the model selection problem has recently attracted a lot of attention [1]. Herethe term “nonparametric” means that the size of a parameter
space (complexity) is not fixed and in theory an infinite numberof parameters (infinite complexity) are considered. If an infinite amount of observed data were available, an infinite numberof parameters would be needed to represent variety of the data.
Actually, however, only a limited number of parameters areneeded because the amount of observed data is limited. The effective complexities of nonparametric models can be automatically adjusted according to observed data. Such nonparametricmodels are essentially different from conventional parametricmodels. In a single nonparametric model, an infinite number ofparametric models with different complexities are stochasticallyoverlapped.
In this paper, we propose a nonparametric Bayesian method
for multipitch analysis, which is the basis of music transcription
and music information retrieval (MIR). The method is called infinite latent harmonic allocation (iLHA), and our goal is to esti-

1558-7916/$26.00 © 2011 IEEE

718

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

Fig. 1. Advantage of our method: We are not required to specify the numberof spectral bases and the number of harmonic partials in advance. In addition,
we do not have to adjust hyperparameters carefully.

mate multiple fundamental frequencies (F0s) from polyphonicaudio signals. Instead of determining the values of F0s (parameters) and the number of them (complexity) uniquely, our method
estimates a joint posterior distribution of all unknown variables
when amplitude spectra of musical audio signals are given as observed data. We formulate nested infinite GMMs for observed
spectra by using nonparametric priors called Dirichlet processes(DPs). These models can be obtained by taking the limit of thenested finite GMMs proposed by Goto [2] and Kameoka et al.
[3] as the number of mixtures goes to infinity. More specifically,each spectral strip is allowed to contain an unbounded number ofsound sources (harmonic structures), each of which is allowed
to contain an unbounded number of harmonic partials. An important problem is that the parameters of the DPs (called hyperparameters) should be given appropriately because they affectthe effective number of mixtures.
To avoid hyperparameter tuning, our models are formulated
in a hierarchical Bayesian manner by putting prior distributions(called hyperprior distributions) on influential hyperparameters. Conventionally, we need to specify the hyperparametersof Dirichlet prior distributions on the relative weights ofharmonic partials [2], [3]. Although these hyperparametersstrongly impact the accuracy of F0 estimation, it is difficult
to optimize them by hand. We instead put noninformative
hyperprior distributions on the hyperparameters of DP priorsof the infinite number of F0s and harmonic partials. This isreasonable because we have little knowledge of the hyperparameters. As shown in Fig. 1, we can completely automate
iLHA by leveraging natural Bayesian treatment of parameters,
complexities, and influential hyperparameters.
The reminder is organized as follows. Section II introduces
related work. Section III compares parametric models of conventional methods and nonparametric models of our method.
Section IV describes a finite version of our method (LHA) andSection V explains our method (iLHA). Section VI reports our
experiments. Section VII concludes this paper.
II. RELATED WORK
Many researchers have applied probabilistic models to multipitch analysis. Goto [2] proposed a probabilistic model for asingle-frame amplitude spectrum (spectral strip) that contains
multiple harmonic structures (see Section III) and used it to estimate the F0s of melody and bass lines from polyphonic audio
signals. Kameoka et al. [3] estimated multiple F0s by using asimilar model for grouping frequency components into multiple

sound sources. Kameoka et al. [4] extended the model by capturing the temporal continuity of harmonic structures. Raphael[5] formulated a HMM based on a large number of chord hypotheses. Cemgil et al. [6] used a dynamic Bayesian network(DBN) to represent the sound generation process, i.e., to associate a music-score level with an audio-signal level. Raczyński
et al. [7] also used a DBN to model temporal dependencies between musical notes. Emiya et al. [8] proposed a probabilistic
model that jointly represents spectral envelopes and harmonic
partials.
Recently, nonnegative matrix factorization (NMF) [9] hasbeen considered to be promising. It regards time–frequency
spectra as a nonnegative matrix and decomposes it into the
product of two nonnegative matrices, one corresponding to a setof spectral bases and the other corresponding to a set of temporal
activations. Smaragdis et al. [10] pioneered the use of NMF formusic transcription. Virtanen et al. [11] and Peeling et al. [12]proposed Bayesian extensions of NMF. Raczyński et al. [13] and
FitzGerald et al. [14] proposed harmonicity constraints for spectral bases, and Bertin et al. [15] further introduced smoothness
constraints for temporal activations. Vincent et al. [16] proposeda method of training spectral bases from audio signals of isolated
tones and adapting them to target polyphonic audio signals.
Cont [17] developed NMF with sparsity constraints for real-time
pitch tracking. Several variants of NMF—such as the complex
NMF proposed by Kameoka et al. [18], the Itakura–Saito (IS)divergence NMF proposed by Févotte et al. [19], and the gamma
process NMF proposed by Hoffman et al. [20]—have been
applied to spectrogram decomposition, but F0s have not been
estimated from the spectral bases thus obtained.
Many other approaches have been also proposed (see [21]for a review). For example, Marolt [22] and Klapuri [23]proposed auditory-model-based methods that use a peripheral
hearing model. Computationally efficient approaches basedon harmonic sums [24] and correlograms [25] have also been
investigated. Pertusa and Iñesta [26] proposed a spectral-peakclustering method. Bello et al. [27] tackled grouping of frequency components by using a heuristic set of rules.
There have been attempt to estimate F0 contours of melody
lines (vocal parts) from polyphonic audio signals. Dressler [28]used instantaneous frequency estimation, sinusoidal extraction,
psychoacoustics, and auditory stream segregation. Ryynänen
and Klapuri [29] formulated a HMM based on acoustic and musicological modeling, and Durrieu et al. [30] proposed a statistical method of extracting the main melody by using source/
filter models. Poliner et al. [31] have reported a comparativeevaluation of several approaches.
Most methods mentioned above can achieve good results if
the number of sound sources and/or manual parameters are appropriately specified. However, it is difficult to always bring out
the full potential of these methods in practice.
III. PROBABILISTIC MODELSOur method is based on nonparametric Bayesian extensionof conventional finite mixture models proposed by Goto [2] andKameoka et al. [3]. Here we explain the conventional models forobserved spectra and then derive our infinite mixture models byextending the conventional models.

YOSHII AND GOTO: NONPARAMETRIC BAYESIAN MULTIPITCH ANALYZER BASED ON iLHA

719

Fig. 2. Gaussian mixture model for the k th basis (single basis). Each Gaussiancorresponds to a harmonic partial, and the mixing weights represent the relative
strengths of M harmonic partials.

A. Notations
Suppose that given polyphonic audio signals contain bases,
harmonic partials located at inteeach of which consists ofgral multiples of the F0 on a linear frequency scale. Each basis
can be associated with multiple sounds of different temporal positions if these sounds are derived from the same pitch of thesame instrument. We transform the audio signals into wavelet
be the number of frames. Note thatandspectra. Let
are finite integers that in conventional methods are specified in
go to infinity.
advance. Our method considers that andFig. 3. Nested Gaussian mixture model for mixed multiple bases. It is obtainedby mixing multiple Gaussian mixture models in a weighted manner under theassumption of amplitude additivity.
TABLE IMULTIPITCH ANALYSIS METHODS

B. Conventional Finite Models and MAP Estimation
Probabilistic models can evaluate how likely observed data is
to be generated by using a limited number of parameters. Therefore, estimation of multiple F0s corresponds directly to finding
model parameters that give the highest probability to the generation of the observed data (called model training).
Goto [2] first proposed probabilistic models of harmonic
structures by regarding an amplitude spectrum (a spectral stripof a single frame) as a probability density function. As shown
canin Fig. 2, the amplitude distribution of basis
be modeled by a harmonic GMM as follows:(1)
where is a one-dimensional vector indicating a logarithmic
frequency [cents].1 The Gaussian parameters (mean and vari) represent the F0 of basis and the degree of enanceis the relative strength of theergy spread around the F0.
th harmonic partial
in basis . We setto. This means that Gaussians are located to havethe harmonic relationship on the logarithmic frequency scale.
can be precomputed beOne might think that the value ofsinusoidal signals (see Apcause the basis sound consists ofpendix I in [4]). This is true if these sinusoidal signals are stationary, but frequency-modulated sounds (e.g., vibrato) result
because of the uncertainty principle ofin a larger value oftime–frequency resolution.
As shown in Fig. 3, the spectral strip of frame is modeledby mixing harmonic GMMs as follows:(2)
1Linear frequency f in hertz can be converted to logarithmic frequency fin cents as follows: f
f =.

= 1200log ( (440(2

))

where
is a relative strength of basis in frame . Consequently, the polyphonic spectral strip can be represented bymeans of a nested finite GMM.
Several methods that have been proposed for parameter estimation are listed in Table I. Goto [2] proposed a predominant-F0 estimation method (PreFEst) that estimates only relative strengths and by allocating many GMMs ( andare fixed) to cover the entire frequency range as F0 candidates.
Kameoka et al. [3] proposed harmonic clustering (HC), which
estimates all the parameters and selects the optimal number ofbases by using the Bayesian information criterion. Althoughthese methods yielded the promising results, they analyze thespectral strips of different frames independently. Kameoka etal. [4] therefore proposed harmonic-temporal-structured clustering (HTC), which captures the temporal continuity of spectral
bases. Because all the above methods use a maximum a posteriori (MAP) estimation strategy to train the finite models, a priordistribution of relative strengths has a large effect on the accuracy of F0 estimation.
C. Our Infinite Models and Bayesian InferenceandWe would like to discuss the limit of (1) and (2) asdiverge to infinity. There is a reason that taking the infinite
limit is reasonable even though there are a finite number of
discrete pitches (e.g., the standard piano has 88 keys). The
F0s and spectral shapes of many instruments (strings, woodwinds, brasses, etc.) vary infinitely according to playing styles(vibrato, marcato, legato, staccato, etc.), and it is difficult tocapture these variations when using a parametric model of fixedcomplexity.
Although there are theoretically infinite number of mixing
and, in theweights

720

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

finite amount of observed data in practice there are a finite
number of bases and a finite number of harmonic partials. Mostof mixing weights must therefore be almost equal to zero. In
other words, only a limited number of bases and a limited
number of harmonic partials are allowed to become active.
To realize such “sparse” GMMs, we put nonparametric priordistributions on mixing weights as sparsity constraints. We
developed a method of Bayesian inference called iLHA to train
the nested infinite GMMs (see Section V).
1) Definition of Observed Data: In the context of Bayesianinference we need to explicitly define the observed data from thestatistical viewpoint. More specifically, we regard each spectral
strip as a histogram of observed frequencies as in [32]. If a spectral strip at frame
has amplitude at frequency ,times in frame
we assume that frequency was observed, where is a scaling factor of wavelet spectra. In other words,
we suppose there are countable frequency “particles” (soundquanta), each corresponding to an independent and identically
distributed (i.i.d.) observation. Note that there is a nontrivial
issue in determining the value of (see Section III-C3). Assuming that amplitudes are additive, we can consider each obpartials in one ofservation to be generated from one ofbases.
Let the total observations over all
frames be represented, whereis a set of observed freby
in frame .
is the number
quenciesof frequency observations (i.e., the sum of spectral amplitudes
is aover all frequency bins in frame ) and
one-dimensional vector that represents an observed frequency.
be the total number of observations overWe letall frames.
be
Let the total latent variables corresponding to, wheresimilarly represented by.
is a-dimensional vector, takes a value of 1 and thein which only one entry,
others take values of 0 when frequency
is generated fromof basis.
partial2) Positioning of Our Method: Our method can be viewed
as an extension of a well-known topic modeling method called
latent Dirichlet allocation (LDA) [33]. LDA was developed asa Bayesian extension of probabilistic latent semantic analysis(pLSA) [34] in the field of natural language processing. In LDA,
each document is represented as a weighted mixture of multiple
topics that are shared over all documents contained in observed
data. Our method similarly represents frames as weighted mixtures of bases. An important difference between our method and
LDA, however, is that iLHA represents each basis as a continuous distribution (a GMM) on the frequency space while LDA
represents each topic as a discrete distribution over words (a setof unigram probabilities).
Another relevant extension of pLSA is probabilistic latentcomponent analysis (PLCA) [35]. PLCA has been applied tosource separation by assuming the time–frequency spectrogram
to be a two-dimensional histogram of sound quanta. A major difference between our method and PLCA is that iLHA is based ona continuous distribution on the frequency space at each frame
while PLCA is based on a two-dimensional discrete distribution
on the space of frame-frequency pairs.

Our method is also similar to the standard NMF [10] basedon temporal exchangeability of spectral strips (see Table I). Our
method simultaneously trains GMMs of all frames contained in
the observed spectra. In other words, if we permute a temporalsequence of spectral strips, the same results would be obtained.
Although such temporal modeling is not appropriate for music,
it is known to work well in practice.
As discussed above, we fuse the topic modeling framework
into the NMF-style decomposition. This is reasonable because
any (local) maximum-likelihood solution of pLSA is proven to
be a solution of NMF that uses Kullback–Leibler (KL) divergence as a cost function [36]. In addition, we propose a nonparametric Bayesian extension.
3) Limitations of Our Method: The amplitude quantizationand i.i.d. assumption are not justified in a physical sense. The
amplitudes at the integral multiples of a F0 are correlated toeach other when they were generated from a single harmonicsound. Besides this, there is arbitrariness in determining the total(the scaling factor multiplied tonumber of observations
raw wavelet spectra). The larger is, the more observations wehave, resulting in a more compact posterior distribution becauseof reduced uncertainty. This criticism can be applied not only to
topic models like [32], [35] but also to probabilistic models of
NMF with KL divergence. This NMF assumes the value of amplitude to follow a Poisson distribution that is defined over nonnegative integers and has no scale parameter. Note that another
NMF with IS divergence [19] does not have such a problem
because it assumes the value of power (squared amplitude) tofollow an exponential distribution that is defined over nonnegative real numbers and has a scale parameter. Therefore, NMF
with IS divergence is scale invariant.
This is more problematic in the context of “nonparametric”
Bayesian inference because the larger number of observations
allows iLHA to activate more relatively small mixture components (i.e., bases and harmonic partials). We therefore need toperform a thresholding process according to the value of after
training the weights of bases. In our experiments, the accuracyof multipitch analysis little varied if we changed the value of(see Section VI).
Another limitation is that our method represents harmonicsounds in an oversimplified manner. We assume that harmonicsounds consist only of several sinusoidal signals corresponding
to harmonic partials. Actually, however, measurable noisy components are widely distributed along the frequency axis even if
the target musical pieces are played only by pitched instruments.
iLHA is thus forced to use too many harmonic GMMs to represent such noisy components. This is another reason that weneed the thresholding process in the end.
IV. LATENT HARMONIC ALLOCATIONThis section explains LHA, the finite version of iLHA, as apreliminary step to deriving iLHA. LHA deals with nested finite
GMMs described in Section III in a Bayesian manner. First,
we mathematically represent the LHA model by putting priordistributions on unknown variables. Then, we explain a training
method of estimating posterior distributions.

YOSHII AND GOTO: NONPARAMETRIC BAYESIAN MULTIPITCH ANALYZER BASED ON iLHA

Fig. 4. Graphical representation of nested finite Gaussian mixture models for
LHA. First, finite sets of mixing weights,  and  , are stochastically generated
according to Dirichlet prior distributions. At the same time, KM Gaussian distributions are stochastically generated according to a Gaussian–Wishart priordistribution. Then one of M harmonic partials in one of K bases is stochastically selected as a latent variable z according to multinomial distributions
defined by  and  . Finally, frequency x is stochastically generated according
to a Gaussian distribution specified by z .

A. Model Formulation
Fig. 4 shows a graphical representation of the LHA model.
The full joint distribution is given by(3)
where the first two terms are likelihood functions and the otherthree terms are prior distributions. The likelihood functions aredefined as
(4)

721

(MCMC) [38], which sequentially generates samples (the concrete values of unknown variables) from the true posterior distribution in a stochastic way by constructing a Markov chainthat has the target distribution as its equilibrium distribution. It
is generally difficult, however, to tell whether or not a Markov
chain has reached a stationary distribution from which we canget samples within an acceptable error.
In the VB framework, we introduce a variational posterior
and make it close to the true posdistribution
iteratively. Here we assume that theteriorvariational distribution can be factorized as(9), we use a variational version of
To optimize
the expectation–maximization (EM) algorithm [37]. We iterate
VB-E and VB-M steps alternately until a variational lower
converges as follows:bound of evidence
(10)
(11)
C. Variational Posterior Distributions
We derive the formulas for updating variational posterior distributions according to (10) and (11).
1) VB-E Step: An optimal variational posterior distributionof latent variables can be computed as follows:(5)
Then, we introduce conjugate priors as follows:
(12)
(6)
where

is defined as

(7)
(13)
(8)
and
are products of Dirichlet distributions andis a product of Gaussian–Wishart distributions.
andare normalization factors, andandare hyperparameters. We let and sum to unity, respectively. andare often called concentration parameters.
, ,
, and
are also hyperparameters:
is a Gaussian mean, is a scaling
is a scale matrix, and is afactor of the precision matrix,
degree of freedom.
where

B. Variational Bayesian Inference
The goal of Bayesian inference is to compute a true posterior distribution of all unknown variables:.
Because analytical calculation of the posterior distribution is
intractable, we use an approximation technique, called variational Bayes (VB) [37], that limits the posterior distribution toan analytical form and optimizes it iteratively in a deterministic
way. Another possible technique is Markov chain Monte Carlo

where. Consequently,
multinomial distributions given by

is obtained as
(14)

where
is called a responsibility
isthat indicates how likely it is that observed frequency
be thegenerated from harmonic partial of basis . Let
number of frequencies that were generated from harmonic partial of basis in frame . This number and its expected value
can be calculated as follows:(15)
For convenience in executing the VB-M step, we calculate
several sufficient statistics as follows:
(16)

722

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

(17)
(18)
2) VB-M Step: Similarly, an optimal variational posterior
is given by
distribution of parameters

where is the digamma function, which is defined as the logarithmic derivative of the gamma function.
D. Variational Lower Bound
To judge convergence, we examine the increase of the variational lower bound. Its maximization is inextricably linked withminimization of the KL divergence between the true and variational posteriors. Let be the lower bound given by(19)
This distribution can be factorized into the product of posterior
distributions of respective parameters as follows:
(34)
(20)
Since our model is based on the conjugate prior distributions,each posterior distribution has the same form of the corresponding prior distribution as follows:
(21)
(22)
(23)
where the variational parameters are given by
(24)
(25)
(26)
(27)
(28)
(29)
where we introduced a dot notation for improved readability.
We let dot “ ” denote the sum over that index. For conveniencein the subsequent sections, we also introduce notations using
comparison operators ( and ). For example, we write(30)
The three terms of (13) can therefore be calculated as follows:(31)

The calculation of these terms is described in Appendix I.
V. INFINITE LATENT HARMONIC ALLOCATIONOur goal is to formulate and train nested infinite GMMs
without model selection and hyperparameter tuning. To do this,we consider the limit of the nested finite GMMs described inandapproach infinity. In addition,
Section IV as both
we put noninformative hyperprior distributions on influential
hyperparameters in a hierarchical Bayesian manner and thencalculate the posterior distributions of those hyperparameters.
As a result of Bayesian inference, likely values are given
large posterior probabilistic densities and unlikely values are
given small densities. Such informative posterior distributions naturally emerge from noninformative prior distributions
as polyphonic spectra are observed. That is, uncertainty isdecreased by getting additional information. In the end we estimate F0s by taking MAP values of the posterior distributions.
A. Mathematical Preparation
We explain the Dirichlet process (DP) and the hierarchical
Dirichlet process (HDP), which can be used as nonparametric
Bayesian priors in our infinite models. In this section, mathematical symbols are defined according to the custom. Therefore,
the definition is valid only in this section.
1) Dirichlet Process: The DP and its extensions play important roles in the theory of Bayesian nonparametrics [39]. Formally introduced by Ferguson [40] in 1973, in the past 10 years
it has often been used as a building block of infinite mixture
models.
A formal definition of the DP is that its marginal distributions must be Dirichlet distributed [40]. Let be a positive real
number andbe a distribution over a sample space . We saya random distribution over is DP distributed with concenif
tration parameter and base measure

(32)
(35)
for any finite measurable partitionof . The
DP is thus a distribution over distributions. This is written
(33)

(36)

YOSHII AND GOTO: NONPARAMETRIC BAYESIAN MULTIPITCH ANALYZER BASED ON iLHA

723

Fig. 5. Stick-breaking construction of the Dirichlet process. Starting with astick of length 1, we break it at 
Beta(1; ) and assign  to be the lengthof the stick we just broke off. We obtain the infinite number of mixing weights, ;  ; . . . ;  , by breaking the remaining portion recursively.



f

g

Fig. 7. Overview of hierarchical Dirichlet process. G becomes an infinite-dimensional discrete distribution when H is a continuous distribution. The smalleris, the fewer atoms in G occupy most of its total probability mass. Thismeans that G becomes more sparse.

Fig. 6. Discretization property of the Dirichlet process. G becomes an infinite-dimensional discrete distribution when G is a continuous distribution. The
smaller is, the fewer atoms in G occupy most of its total probability mass. Thismeans that G becomes more sparse.

Then, a concrete sample

is drawn fromas follows:(37)

An alternative constructive definition of the DP is known as
the stick-breaking construction (SBC) [41]. As illustrated inFig. 5, a random distribution can be written explicitly as acountably infinite sum of point masses (“atoms”):
(38)
(39)
is the Dirac delta function that diverges to positive
where, is otherwise equal to 0, and integrates to 1
infinity atof is given bywith respect to . The point mass
(40)
(41)
The distribution on the infinite number of mixing weights
is often written, where theletters stand for Griffiths, Engen, and McCloskey.
must be a disAn important property of the DP is thatcrete distribution. As shown in Fig. 6, is an infinite-dimenis a continuous distribusional discrete distribution whention. The DP can therefore be used as a prior distribution to formulate an infinite mixture model. In case of an infinite GMM(iGMM), for example, is a space of Gaussians (i.e., a space ofis usually set to a Gaussian–Wishartmeans and variances).
distribution, which is a conjugate prior distribution over Gaussians. drawn from the DP is also a distribution over Gaussians. Every time an observation is generated, Gaussian

is drawn fromof Gaussians, whereis selected from the infinite number
according to their probabilities. This is a straightforward extension of a conventional finite GMM.
Several extensions increasing a degree of freedom of the standard DP have been proposed. For example, a beta two-parameter
process [42] is obtained when
(42)

where positive real numbers and are adjustable parametersof the beta distribution.
2) Hierarchical Dirichlet Process: We discuss how to simultaneously train tied infinite mixture models when observed dataconsists of multiple groups, e.g., spectral strips (frames). Herea set of component distributions should be shared across mixture models trained for different groups. Such parameter tying
enables us to directly compare compositions of different groups
in terms of mixing weights of component distributions. This issimilar to vector quantization (VQ) [43]. Let be the numberof groups. In this setting, it is natural to use a DP for modelingobserved data of each group as follows:
(43)
where
is a random distribution on for group .
is a continuous distribution, atoms
A problem is that iffor generating ob(component distributions) drawn fromservations are almost surely disjointed from those drawn from. This is because DPs can independently determine the positions of the countably infinite number of discrete(cardinality ) from the uncountably
atoms
infinite continuous space (cardinality ).
To solve this problem, we use a HDP [44] as a nonparametricprior distribution. As shown in Fig. 7, we consider the base meaitself to be distributed according to a top-level DP as
surefollows:
(44)

724

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

where is a concentration parameter and a base measure over
. In this model,
always becomes a discrete distribution. The
SBC of the top-level DP is given by
(45)
(46)
and
are the point
where
masses and positions of atoms and we have.
Similarly, the SBC of a lower-level DP is given by
(47)
(48)
and eachis selected from . Note thatwhere
can be equal toifbecauseis a discrete distribution. Another direct representation based on determinedby the top-level DP is as follows:(49)
and the hyperparameter controls the
where
and
. Therefore, only point masses
difference between(mixture weights) differ between groups while positions(component distributions) are shared across groups.
A remaining problem is how to adjust the influential hyperparameters and . This problem is often solved by putting vague
gamma hyperprior distributions on these hyperparameters andinferring the posterior distributions.
B. Model Formulation
We explain how to formulate nested infinite GMMs basedon a HDP and generalized DPs by extending the nested finite
GMMs described in Section IV.
First we discuss. An important requirement is thatbasis models [harmonic GMMs represented by (1)] should beframes because each basis
shared as a global set across all
sound has a duration and may appear in different frames while
only its weight varies. The HDP can satisfy this requirement
and we can explain the HDP from the generative point of view.
After an unbounded number of bases are initially generated according to a top-level DP, an unbounded number of bases are selected in each frame according to a frame-specific DP. In practice, a limited number of bases are used to represent a spectral strip because the number of observed frequency particles
is limited. Mathematically speaking, in (6) we consider
infinite-dimensional Dirichlet distributions, which are equivalent to the frame-specific DPs, and assume hyperparameter tobe distributed according to the top-level DP as follows:

where

is a concentration parameter of the top-level DP.

Fig. 8. Graphical representation of nested infinite Gaussian mixture models for
iLHA. First the infinite sets of mixing weights  and  are stochastically generated according to a HDP and beta two-parameter processes (generalized DPs).
At the same time, the infinite number of Gaussian distributions are stochastically generated according to a Gaussian–Wishart prior distribution. Then one of
the harmonic partials contained in one of the bases is stochastically selected asa latent variable z according to multinomial distributions defined by  and . Finally, frequency x is stochastically generated according to a Gaussian
distribution specified by z .

Now we discuss. Because each basis is allowed toconsist of a unique infinite set of harmonic partials (basis models
are independent of each other), instead of (7) we can use beta
two-parameter processes as follows:
(52)
(53)
where is a positive real number and we letand
sum tounity. Note that we used the size-biased permutation propertyof the SBC to encourage lower harmonic partials to have larger
weights because roughly speaking, the weights of harmonic partials of an instrument sound decrease exponentially.
Because hyperparameters , , , and are influential, we
put hyperprior distributions on them as follows:
(54)
(55)
(56)
(57)
and
are shape and rate parameters of
where
the gamma distributions.
andare parameters of the beta
distribution. These distributions are set to be vague (,, andin our experiments
described in Section VI).
Fig. 8 shows a graphical representation of the iLHA model.
The full joint distribution is given by

(58)

(50)

where
is given by plugging (52) into (5) andis given by (6).
andare defined according to(50) and (52) as follows:

(51)

(59)
(60)

YOSHII AND GOTO: NONPARAMETRIC BAYESIAN MULTIPITCH ANALYZER BASED ON iLHA

725

We then introduce auxiliary variables by using a techandbenique called data augmentation [45]. Letandbe positive intebeta-distributed variables and,, andgers that satisfy. We can augment (63) as follows:

Fig. 9. Graphical representation of collapsed nested infinite mixture models
for iLHA. After the original parameters  , ~ ,  , and 3 are integrated out, theauxiliary variables  ,  , s , and t are introduced to set up conjugacy between
hyperprior distributions and a marginalized likelihood function.

(64)

C. Collapsed Variational Bayesian Inference
There are two problems in training the HDP mixture model.
The first problem is that VB needs to assume the independence
between latent variables and parameters to factorize a posterior
distribution as in (9). This assumption is sometimes too strong
and leads to incorrect posterior approximation. The second
problem is that applying VB to hierarchical Bayesian models
that have no conjugacy between priors and hyperpriors is
generally difficult.
To solve these problems, we use a sophisticated version of
VB called collapsed variational Bayes (CVB) [45]. It instead
assumes independence between individual latent variables ina “collapsed” space in which parameters are integrated out(marginalized out). This is reasonable because the dependence
between individual latent variables in the collapsed space isgenerally much weaker than the dependence between a set ofparameters and a set of latent variables in the non-collapsed
space. In addition, we introduce auxiliary variables to apply
CVB to hierarchical Bayesian models.
Fig. 9 shows a graphical representation of a collapsed iLHA
model. Integrating out , , , and , we obtain the marginal
distribution given by(61)
The first term of (61) can be easily calculated by leveraging
conjugacy between
andas follows:
(62)
where
and
are normalization factors,,of prior and posterior Gaussian–Wishart distributions.
are obtained by substitutingforin caland
culating (26), (27), and (29). Similarly, the second term of (61)can be calculated by leveraging conjugacy between
andas follows:

where [] denotes a Stirling number of the first kind. We canconfirm that (64) reduces to (63) by marginalizing out auxiliary
variables , , , and . The augmented marginal distribution is
given by(65)
To apply CVB to approximate the true posterior distribution, we assume that the variational
posterior distribution can be factorized as follows:(66)

where we assumed independence between hyperparameters,
auxiliary variables, and elements of . We also use an approximation technique called variational posterior truncation. Morewhenandspecifically, we assume that
. In practice, we set
andto sufficiently large
integers. This does not mean that effective model complexities
are fixed in advance. The larger the truncation levels we use,
the more the accurate approximations we obtain.
, we use a variational
To optimize
EM algorithm that iterates the following steps:
(67)
(68)
(69)
(70)
where

denotes a set of indices withoutand .

D. Variational Posterior Distributions
(63)
where

is the gamma function.

We derive the formulas for updating variational posterior distributions according to (67)–(70).
1) CVB-E Step: An optimal variational distribution of canbe obtained as the product of multinomial distributions. The

726

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

posterior probability thatwas generated from themonic partial of basis is given by

th har-(71)
whereis the geometric averageand is the Student-t distribution defined by the three parame,, and.
is given by
ters

(82)
To optimize the variational posterior distributions, we need to
calculate the expectations of these variables. If a random variwith shape parameter and rate
able follows
andparameter , its expectations are given by. If follows
with parameters and , its expectations are given byand. Note that the distributions given
by Equations (79)–(82) are conditioned by . The expectations
must therefore be averaged over . For example, we now havethe following conditional expectation:
(83)

(72)
where
,
,
, and
are obtained according tois substituted forand the sums(26)–(29) in which
are calculated without. Each term of (71) can be approximated efficiently by using first-order and second-order Taylor
expansions [45]–[47].
Equation (71) calculates the geometric averages of three
predictive distributions under posterior distributions. These
predictive distributions are derived from an infinite-dimensional Dirichlet distribution (a DP for an infinite mixture of
iGMMs), stick-breaking construction (a DP for an iGMM), anda Gaussian distribution. Interestingly, this corresponds to (13)based on the geometric averages of three likelihood functions
under posterior distributions. This implies that CVB is morerobust to the local-optima problem than standard VB is.
2) CVB-M Step: We can optimize the variational posterior
distributions of the hyperparameters analytically by optimizing
those of the auxiliary variables. First, , , and are gamma
distributed as follows:over , but
We use Taylor expansion to average
diverges to negative infinity muchthe digamma function
faster than the logarithmic function does in the vicinity of theorigin. To solve this problem, we use a method that treats theexactly and applies second-order approxicase[45]. We can similarly average themation when
following conditional expectations:
(84)
(85)
(86)
To estimate F0s in the end, we explicitly compute the variational posterior distributions of the integrated-out parametersand . To do this, we need to execute the standard VB-M stepobtained in the CVB-E step.
once using
E. Variational Lower Bound
As in LHA, we monitor the increase of the variational lower, which is given bybound of evidence

(73)
(74)
(75)
and

and

are beta distributed as follows:
(87)

(76)
(77)
Then

and

are beta distributed as follows:

VI. EVALUATION
(78)
(79)

and

The calculation of these terms is described in Appendix II.

and are multinomial distributed as follows:
(80)
(81)

This section reports the results of two comparative evaluation
experiments. We compared LHA and iLHA with PreFEst andHTC because these four methods are based on the same idea
for modeling harmonic structures. Using a different data set,
we then compared iLHA with NMF-based methods and othermethods. In the latter experiment, we investigated how significantly the value of the scaling factor (i.e., how many frequency particles are assumed to be observed in total) affects theaccuracy of multipitch analysis.

YOSHII AND GOTO: NONPARAMETRIC BAYESIAN MULTIPITCH ANALYZER BASED ON iLHA

TABLE IIFRAME-LEVEL F-MEASURES OF F0 DETECTION

727

TABLE III
FRAME-LEVEL ACCURACY OF F0 DETECTION

A. Comparison with Conventional Parametric Methods
1) Experimental Conditions: We evaluated LHA and iLHA
on a test that was used in [4] and consisted of eight pieces of
piano and guitar solo performances excerpted from the RWC
music database [48]. The first 23 s of each piece were used forevaluation. Spectral analysis with a 16-ms time resolution was
conducted using a wavelet transform with Gabor wavelets. The
correct values and temporal positions of actual F0s were prethe
pared by hand as ground truth. Denoting by , , andrespective numbers of ground-truth, estimated, and correct F0s
on frame , we calculated the following frame-level recall andprecision rates and F-measure for each piece:(88)
and we averaged each of these measures over all pieces.
The prior and hyperprior distributions of LHA and iLHA
were set to noninformative distributions. In LHA,
andwere set to 60 and 15. In iLHA,
andwere also set to 60
and 15. iLHA is not sensitive to these values, and no other tuning
was needed for either method. To output F0s at each frame, we
extracted bases whose expected weights were over a thresholdthat was optimized as in [4].
For comparison, we referred to the PreFEst and HTC experimental results reported in [4]. Although the ground-truth data inthat study was slightly different from ours, it was close enoughfor roughly evaluating performance comparatively. The numberof bases, priors, and weighting factors of the PreFEst and HTC
were carefully tuned to optimize the results. Although this is notrealistic, the upper bounds of potential performance were investigated in the literature.
2) Experimental Results: The results listed in Table IIshow that the performance of iLHA closely approached andsometimes surpassed that of HTC. This is consistent with theempirical findings of many studies on Bayesian nonparametrics
that nonparametric models were competitive with optimallytuned parametric models. HTC outperformed PreFEst because
HTC can appropriately deal with temporal continuity of spectral bases. This implies that incorporating temporal modelingwould improve the performance of iLHA.
The results of LHA were worse than those of iLHA because
LHA is not formulated in a hierarchical Bayesian manner andrequires precise priors. In fact, we confirmed that the results of
PreFEst and HTC based on MAP estimation were drasticallydegraded when using noninformative priors. Automated iLHA,
in contrast, stably showed the good performance.

We found that model flexibility can be greatly enhanced bymaking time-consuming fine tuning unnecessary. Conventional
studies assumed that appropriate prior knowledge is required
to constrain flexibility (called regularization). By using a truly
flexible hierarchical model based on Bayesian nonparametrics,
however, we can let the data speak for itself. This naturally results in optimal performance.
B. Comparison With NMF-Based Methods and Other Methods
1) Experimental Conditions: We then evaluated iLHA on atest set that was used in [16] and consisted of 50 pieces of pianosolo performances excerpted from the MAPS piano database [8].
The first 30 s of each piece were used for evaluation. Spectral analysis with a 10-ms time resolution was conducted using a Gaborwavelet transform. The value of
was increased to 88, (thenumber of notes in a standard piano) because the piano pieces
were much sophisticated than those used in first experiment. The
were equal to those used intime resolution and the value of[16], and performance was evaluated in terms of F-measures.
For comparison, we referred to the experimental results ofseven methods reported in [16]. We compared iLHA with four
NMF-based methods: one using no constraints, one using harmonicity constraints (a subset of [13]), one using harmonicity
and source-filter constraints [14], and one using harmonicity andspectral smoothness constraints [16]. Note that only the last one
was manually tuned to yield the best results (the effect of hyperparameter tuning was investigated in [16]). We also compared
it with a method based on harmonic sums [24], a method basedon correlograms [25], and a method based on spectral peak clustering [26].
2) Experimental Results: The results listed in Table III showthat iLHA was the second best among the seven methods. Although the best variant of NMF gained the better F-measure(67.0%) than iLHA did (61.2%), we can say that well-automatediLHA is still competitive because it is reported that non-optimal
settings deteriorated the performance of NMF moderately [16].
The F-measure of iLHA (61.2%) was close to that of NMF
using only harmonicity constraints (60.5%). As discussed in
Section III-C2, pLSA and PLCA are proven to have a close connection to NMF. Therefore, the similarity between iLHA basedon harmonic GMMs and NMF based on harmonicity constraints
was experimentally and theoretically supported. In addition, thedifference between NMF using only harmonicity constraints
and NMF adding spectral smoothness constraints implies thatthe performance of iLHA would be improved by incorporating
spectral smoothness modeling.

728

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012

It is interesting that in almost all methods the was higher
than the . This means that there were many F0s that were hardto detect because of the complex overlapping of multiple F0s.
To solve this problem, more accurate spectral modeling would
be required by removing the assumption of amplitude additivity
that forms a basis of iLHA and NMF.
3) Impact of Scaling Factor: We investigated the impact of
the scaling factor described in Section III-C. We tested three. The similarity of respective
different values:
F-measures—61.2%, 60.6%, and 60.1%—indicates that the results are not sensitive to the value of the scaling factor . The
automatic optimization of would be an interesting research
topic that tackles the limitation of many methods based on theassumption of amplitude quantization.
VII. CONCLUSION
This paper presented a novel statistical method for detecting
multiple F0s in polyphonic music audio signals. In this method,which is called iLHA and is the first to apply Bayesian nonparametrics to multipitch analysis, we formulated nested infinite
GMMs that represent polyphonic spectral strips in a hierarchicalnonparametric Bayesian manner. More specifically, each spectral strip is allowed to contain an unbounded number of spectral bases, each of which can contain an unbounded number of
harmonic partials. The method was fully automated by putting
noninformative hyperprior distributions on influential hyperparameters except for the final thresholding process. The joint posterior distribution of all unknown variables can be inferred efficiently according to the VB framework. In our experiments
comparing iLHA with the state-of-the-art methods manually optimized by trial and error, we found that iLHA is competitive
enough and there is room for improvement based on modelingof temporal continuity and spectral smoothness. One interesting
future direction is to use MCMC methods such as Gibbs sampling and more efficient variants for training the iLHA model.
Bayesian nonparametrics is a powerful framework avoiding
the model selection problem faced in various areas of music
information retrieval (MIR). For example, how many sections
are required for structuring a musical piece? How many groups
are required for clustering listeners according to their tastes ormusical pieces according to their contents? We can avoid these
problems by assuming that in theory there is an infinite numberof objects (sections or groups) behind available observed data.
Unnecessary objects are automatically removed from consideration through statistical inference. Hoffman et al. recently successfully applied this framework to the calculation of musicalsimilarity [49] and the detection of repeated patterns [32], and
we also plan to use this powerful framework in a wide range ofapplications.
APPENDIX I
The nine terms of the variational lower bound of LHA in (34)
can be calculated as follows:

where the fifth and last terms can be obtained as follows:

APPENDIX II
The 13 terms of the variational lower bound of iLHA in (87)can be calculated as follows:YOSHII AND GOTO: NONPARAMETRIC BAYESIAN MULTIPITCH ANALYZER BASED ON iLHA

where
and
mean the first-order and second-order approximations based on Taylor expansion (see [45]–[47]).
ACKNOWLEDGMENT
The authors would like to thank Dr. H. Kameoka (The University of Tokyo/NTT, Japan) for providing the ground-truthtranscriptions of eight pieces included in the RWC music
databases [48]. They would also like to thank Dr. V. Emiya(INRIA, France) for allowing them to use the valuable MAPS
piano database [8].
REFERENCES[1] P. Orbanz and Y. W. Teh, “Bayesian nonparametric models,” in Encyclopedia of Machine Learning. New York: Springer, 2010.
[2] M. Goto, “A real-time music scene description system: Predominant-F0 estimation for detecting melody and bass lines in real-world
audio signals,” Speech Commun., vol. 43, no. 4, pp. 311–329, 2004.
[3] H. Kameoka, T. Nishimoto, and S. Sagayama, “Separation of harmonic
structures based on tied Gaussian mixture model and information criterion for concurrent sounds,” in Proc. IEEE Int. Conf. Acoust., Speech,
Signal Process. (ICASSP), 2004, vol. 4, pp. 297–300.

729[4] H. Kameoka, T. Nishimoto, and S. Sagayama, “A multipitch analyzer
based on harmonic temporal structured clustering,” IEEE Trans. Audio,
Speech, Lang. Process., vol. 15, no. 3, pp. 982–994, Mar. 2007.
[5] C. Raphael, “Automatic transcription of piano music,” in Proc. 3rd Int.
Conf. Music Inf. Retrieval (ISMIR), 2002, pp. 161–166.
[6] A. T. Cemgil, H. J. Kappen, and D. Barber, “A generative model formusic transcription,” IEEE Trans. Audio, Speech, Lang. Process., vol.
14, no. 2, pp. 679–694, Mar. 2006.
[7] S. A. Raczyński, E. Vincent, F. Bimbot, and S. Sagayama, “Multiple
pitch transcription using DBN-based musicological models,” in Proc.
11th Int. Conf. Music Inf. Retrieval (ISMIR), 2010, pp. 363–368.
[8] V. Emiya, R. Badeau, and B. David, “Multipitch estimation of pianosounds using a new probabilistic spectral smoothness principle,” IEEE
Trans. Audio, Speech, Lang. Process., vol. 18, no. 6, pp. 1643–1654,
Aug. 2010.
[9] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in Proc. Adv. Neural Inf. Process. Syst. (NIPS), 2000, pp.
556–562.
[10] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for
polyphonic music transcription,” in Proc. IEEE Workshop Applicat.
Signal Process. Audio Acoust. (WASPAA), 2003, pp. 177–180.
[11] T. O. Virtanen, A. T. Cemgil, and S. J. Godsill, “Bayesian extensions tononnegative matrix factorisation for audio signal modelling,” in Proc.
IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2008, pp.
45–48.
[12] P. H. Peeling, A. T. Cemgil, and S. J. Godsill, “Generative spectrogram
factorization models for polyphonic piano transcription,” IEEE Trans.
Audio, Speech, Lang. Process., vol. 18, no. 3, pp. 519–527, Mar. 2010.
[13] S. A. Raczyński, N. Ono, and S. Sagayama, “Multipitch analysis withharmonic non-negative matrix approximation,” in Proc. 6th Int. Conf.
Music Inf. Retrieval (ISMIR), 2007, pp. 381–386.
[14] D. FitzGerald, M. Cranitch, and E. Coyle, “Extended nonnegative
tensor factorisation models for musical sound source separation,”
Comput. Intell. Neurosci., vol. 2008, 2008.
[15] N. Bertin, R. Badeau, and E. Vincent, “Enforcing harmonicity andsmoothness in Bayesian non-negative matrix factorization applied topolyphonic music transcription,” IEEE Trans. Audio, Speech, Lang.
Process., vol. 18, no. 3, pp. 538–549, Mar. 2010.
[16] E. Vincent, N. Bertin, and R. Badeau, “Adaptive harmonic spectral decomposition for multiple pitch estimation,” IEEE Trans. Audio, Speech,
Lang. Process., vol. 18, no. 3, pp. 528–537, Mar. 2010.
[17] A. Cont, “Realtime multiple pitch observation using sparse non-negative constraints,” in Proc. 7th Int. Conf. Music Inf. Retrieval (ISMIR),
2006, pp. 206–211.
[18] H. Kameoka, T. Nishimoto, and S. Sagayama, “Complex NMF: A newsparse representation for acoustic signals,” in Proc. IEEE Int. Conf.
Acoust., Speech, Signal Process. (ICASSP), 2009, pp. 45–48.
[19] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music
analysis,” Neural Comput., vol. 21, no. 3, pp. 793–830, 2009.
[20] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametric matrix
factorization for recorded music,” in Proc. 27th Int. Conf. Mach. Learn.
(ICML), 2010, pp. 439–446.
[21] Signal Processing Methods for Music Transcription, A. Klapuri and M.
Davy, Eds. New York: Springer, 2010.
[22] M. Marolt, “A connectionist approach to transcription of polyphonicpiano music,” IEEE Trans. Multimedia, vol. 6, no. 3, pp. 439–449,
2004.
[23] A. Klapuri, “Multipitch analysis of polyphonic music and speech
signals using an auditory model,” IEEE Trans. Audio, Speech, Lang.
Process., vol. 16, no. 2, pp. 255–266, 2008.
[24] A. Klapuri, “Multiple fundamental frequency estimation by summing
harmonic amplitudes,” in Proc. 7th Int. Conf. Music Inf. Retrieval(ISMIR), 2006, pp. 216–221.
[25] T. Tolonen and M. Karjalainen, “A computationally efficient multipitch
analysis model,” IEEE Trans. Speech Audio Process., vol. 8, no. 6, pp.
708–716, Nov. 2000.
[26] A. Pertusa and J. M. Iñesta, “Multiple fundamental frequency estimation using Gaussian smoothness,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Process. (ICASSP), 2008, pp. 105–108.
[27] J. P. Bello, L. Daudet, and M. B. Sandler, “Automatic piano transcription using frequency and time-domain information,” IEEE Trans. Audio,
Speech, Lang. Process., vol. 14, no. 6, pp. 2242–2251, Nov. 2006.
[28] K. Dressler, “Extraction of the melody pitch contour from polyphonic audio,” in Proc. 2nd Music Inf. Retrieval Eval. eXchange(MIREX), 2005 [Online]. Available: http://www.musicir.org/evaluation/mirex-results/articles/melody/dressler.pdf

730

IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. 20, NO. 3, MARCH 2012[29] M. P. Ryynänen and A. P. Klapuri, “Transcription of the singing
melody in polyphonic music,” in 7th Int. Conf. Music Inf. Retrieval(ISMIR), 2006, pp. 206–211.
[30] J.-L. Durrieu, G. Richard, B. David, and C. Févotte, “Source/filter
model for unsupervised main melody extraction from polyphonicaudio signals,” IEEE Trans. Audio, Speech, Lang. Process., vol. 18,
no. 3, pp. 564–575, Mar. 2010.
[31] G. E. Poliner, D. P. Ellis, A. F. Ehmann, E. Gómez, S. Streich, and B.
Ong, “Melody transcription from music audio: Approaches and evaluation,” IEEE Trans. Audio, Speech, Lang. Process., vol. 15, no. 4, pp.
1247–1256, May 2007.
[32] M. Hoffman, D. Blei, and P. Cook, “Finding latent sources in recorded
music with a shift-invariant HDP,” in Proc. 12th Int. Conf. Digital
Audio Effects (DAFX), 2009.
[33] D. Blei, A. Ng, and M. Jordan, “Latent Dirichlet allocation,” Mach.
Learn. Res., vol. 3, pp. 993–1022, 2003.
[34] T. Hofmann and J. Puzicha, “Probabilistic latent semantic indexing,”
in Proc. 22nd Int. Conf. Res. Develop. Inf. Retrieval (SIGIR), 1999, pp.
50–57.
[35] M. Shashanka, B. Raj, and P. Smaragdis, “Probabilistic latent variable
models as non-negative factorizations,” Comput. Intell. Neuosci., vol.
2008, 2008.
[36] E. Gaussier and C. Goutte, “Relation between pLSA and NMF andimplications,” in Proc. 28th Int. Conf. Res. Develop. Inf. Retrieval(SIGIR), 2005, pp. 601–602.
[37] H. Attias, “A variational Bayesian framework for graphical models,”
in Adv. Neural Inf. Process. Syst. (NIPS), 2000, pp. 209–215.
[38] R. M. Neal, “Probabilistic inference using Markov chain Monte Carlo
methods,” Dept. of Comput. Sci., Univ. of Toronto, Toronto, ON,
Canada, Tech. Rep. CRG-TR-93-1, 1993.
[39] Y. W. Teh, “Dirichlet processes,” in Encyclopedia of Machine
Learning. New York: Springer, 2010.
[40] T. Ferguson, “Bayesian analysis of some nonparametric problems,”
Ann. Statist., vol. 1, no. 2, pp. 209–230, 1973.
[41] J. Sethuraman, “A constructive definition of Dirichlet priors,” Statist.
Sinica, vol. 4, pp. 639–650, 1994.
[42] H. Ishwaran and L. F. James, “Gibbs sampling methods for
stick-breaking priors,” J. Amer. Statist. Assoc., vol. 96, no. 453,
pp. 161–173, 2001.
[43] R. Gray, “Vector quantization,” IEEE ASSP Mag., vol. 1, no. 2, pp.
4–29, Apr. 1984.
[44] Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei, “Hierarchical
Dirichlet processes,” J. Amer. Statist. Assoc., vol. 101, no. 476, pp.
1566–1581, 2006.
[45] Y. W. Teh, K. Kurihara, and M. Welling, “Collapsed variational inference for HDP,” in Proc. Adv. Neural Inf. Process. Syst. (NIPS), 2007.
[46] J. Sung, Z. Ghahramani, and S.-Y. Bang, “Latent-space variational
Bayes,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 30, no. 12, pp.
2236–2242, Dec. 2008.

[47] J. Sung, Z. Ghahramani, and S.-Y. Bang, “Second-order latent-spacevariational Bayes for approximate Bayesian inference,” IEEE Signal
Process. Lett., vol. 15, pp. 918–921, 2008.
[48] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC music database: Popular, classical, and jazz music database,” in Proc. 3th Int.
Conf. Music Inf. Retrieval (ISMIR), 2002, pp. 287–288.
[49] M. Hoffman, D. Blei, and P. Cook, “Content-based musical similarity
computation using the hierarchical Dirichlet process,” in Proc. 9th Int.
Conf. Music Inf. Retrieval (ISMIR), 2008, pp. 349–354.

Kazuyoshi Yoshii (M’08) received the Ph.D. degree
in informatics from Kyoto University, Kyoto, Japan,
in 2008.
He is currently a Research Scientist at the National
Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Japan. His research interests include probabilistic music analysis, blind source
separation, and Bayesian nonparametrics.
Dr. Yoshii has received several awards including
the IPSJ Yamashita SIG Research Award and theBest-in-Class Award of MIREX 2005. He is amember of the Information Processing Society of Japan (IPSJ) and Institute ofElectronics, Information, and Communication Engineers (IEICE).

Masataka Goto received the D.Eng. degree fromWaseda University, Tokyo, Japan, in 1998.
He is currently the leader of the Media Interaction
Group, National Institute of Advanced Industrial
Science and Technology (AIST), Tsukuba, Japan.
He serves concurrently as a Visiting Professor at
the Institute of Statistical Mathematics, an Associate
Professor (Cooperative Graduate School Program)in the Graduate School of Systems and Information
Engineering, University of Tsukuba, and a ProjectManager of the MITOH Program (the Exploratory
IT Human Resources Project) Youth division by the Information Technology
Promotion Agency (IPA).
Dr. Goto received 25 awards over the past 19 years, including the Commendation for Science and Technology by the Minister of MEXT “Young Scientists’ Prize,” the DoCoMo Mobile Science Awards “Excellence Award in Fundamental Science,” the IPSJ Nagao Special Researcher Award, and the IPSJ Best
Paper Award.