13th International Society for Music Information Retrieval Conference (ISMIR 2012)

INFINITE COMPOSITE AUTOREGRESSIVE MODELSFOR MUSIC SIGNAL ANALYSIS
Kazuyoshi Yoshii
Masataka Goto
National Institute of Advanced Industrial Science and Technology (AIST), Japan{k.yoshii, m.goto}@aist.go.jp
ABSTRACT

Comb-shaped excitations Combinatorial All-pole transfer functions
products⊗

This paper presents novel probabilistic models that can beused to estimate multiple fundamental frequencies (F0s)from polyphonic audio signals. These models are nonparametric Bayesian extensions of nonnegative matrix factorization (NMF) based on the source-filter paradigm, and inthem an amplitude or power spectrogram is decomposedas the product of two kinds of spectral atoms (sources andfilters) and time-varying gains of source-filter pairs. In thisstudy we model musical instruments as autoregressive systems that combine two types of sources—periodic signals(comb-shaped densities) and white noise (flat density)—
with all-pole filters representing resonance characteristics.
One of the main problems with such composite autoregressive models (CARMs) is that the numbers of sources
and filters should be given in advance. To solve this problem, we propose nonparametric Bayesian models based on
gamma processes and efficient variational and multiplicative learning algorithms. These infinite CARMs (iCARMs)
can discover appropriate numbers of sources and filters ina data-driven manner. We report the experimental resultsof multipitch analysis on the MAPS piano database.

IJ bases⊕…

…

White noise (flat density)⊗

IJ gains

J filters
I sources(latent variables) Observed spectrum (latent variables)

Figure 1. Overview of composite autoregressive models:The combinatorial products of I sources and J filters yield
IJ spectral bases, which are activated according to the corresponding time-varying gains at each frame. We take theinfinite limit as both I and J diverge to infinity.
bases are needed to fully represent the timbral variations of instrument spectra (e.g., envelopes) even if
these spectra share the same fundamental frequency(F0). Such an unconstrained increase of model complexity is likely to result in optimization algorithms
easily getting stuck in bad local optima.
2. A post-processing step for estimating the F0s fromindividual bases is required because the F0s are notparameterized for representing the spectral bases. If
the shapes of spectral bases are unconstrained, theresulting bases often deviate from natural harmonic
spectra. This makes F0 estimation difficult and weneed to judge the existence of an F0.
3. The number of bases (model complexity) should be
carefully specified in advance because it has a strong
impact on the decomposition results. Note that thislimitation is closely related to the first. A naive solution is to exhaustively test all possible complexities
and find an optimal value, but such model selection
is often computationally impractical.
As noted above, unconstrained NMF is too flexible for aset of musically-meaningful bases to be induced automatically. Although these limitations have partially been dealt
with in previous studies [1–13], no study has overcome allof them simultaneously in a principled manner.
In this paper we propose infinite composite autoregressive models (iCARMs) (Fig. 1) developed for fusing thefollowing techniques into a unified Bayesian framework:
1. Source-filter factorization (inspired by [1])We further factorize the spectral bases as the combinatorial products of sources and all-pole filters. This

1. INTRODUCTION
Multiple fundamental frequency estimation (a.k.a. multipitch analysis) is the basis of various kinds of music content
analysis. Recently, nonnegative matrix factorization (NMF)
has gained a lot of popularity [1–13]. The standard NMF
approximates an amplitude or power spectrogram (nonnegative matrix) as the product of two nonnegative matrices,
one of which is a compact set of spectral bases and theother of which is a set of the corresponding time-varying
gains [15, 16]. Such low-rank approximation is well justified by the fact that each musical piece consists of only limited kinds of sounds that repeatedly appear. In addition, apractical advantage of NMF is that the bases and gains canbe alternately optimized by using efficient iterative algorithms called multiplicative update (MU) rules. The standard NMF, however, has three fundamental limitations:
1. The spectral bases are time-invariant, and only theirgains vary over time. A large number of independent
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 2012 International Society for Music Information Retrieval.
7913th International Society for Music Information Retrieval Conference (ISMIR 2012)

idea originates in the autoregressive (AR) modelingof speech signals: various vowels can be generated
by changing the shape of the vocal tract (filter) while
keeping the same F0 (source). This factorization enables us to represent a wide variety of instrumental
sounds in terms of two separate aspects (timbre andF0) with reasonable complexity.
2. Harmonicity modeling (inspired by [7] and [9])We represent each source as a comb-shaped functionthat uses an F0 parameter for representing equallyspanned harmonic partials of the same weight. Sincesuch sources are multiplied by acoustically-inspired
AR filters, the relative weights of partials of the bases
are constrained to take realistic values as natural harmonic sounds. In addition, we can directly optimize
the values of the F0s jointly with decomposition.
As proposed in [7, 14], we additionally introduce a
special source representing white noise (flat density).
This enables us to deal with percussive and transient
sounds having widely distributed spectra. Their timbres (envelopes) are characterized by AR filters.
3. Bayesian nonparametrics (inspired by [12])We build nonparametric Bayesian models that can
automatically adjust the numbers of sources and filters needed to factorize a given spectrogram. Rather
than these numbers being specified, the infinite limitof the conventional source-filter NMF [1] is taken as
the numbers of sources and filters diverge to infinity.
We perform sparse learning by introducing infinitedimensional priors in such a way that only limited
numbers of sources and filters are actually activated.

MethodKameoka [1]Badeau [2]Durrieu [3]Virtanen [4]Carabias-Orti [5]Heittola [6]Yasuraoka [7]Hennequin [8]Proposed

Divergence
IS
IS
IS
KL
KL
KL
KL
Beta (0.5)
KL or IS

Sources #
I
H
I
(H)
I
I
H
I
I∗
H+N
I∗
I
H+N ∞

Filters #
AR
J
MA
J
J
J
J
J
AR
J
ARMA J ∗AR∞(H: harmonic sources, N: noise source, -: others, ∗: varying over time)

Table 1. Several variants of source-filter NMF
As shown in Table 1, for example, Kullback-Leibler (KL)[15] and Itakura-Saito (IS) [16] divergences have been used
intensively. Some studies used beta divergence [17], which
includes KL and IS divergences as special cases.
In the context of audio modeling, although IS-NMF isjustified in theory (see [16] and Section 3.2.1), KL-NMF
often yields better results in a maximum likelihood estimation setting. One main reason is that the nonconvexity of IS
divergence makes it difficult for gradient-descent-type optimization algorithms to find global optima. Note that nocomparative tests have been conducted under a Bayesianestimation setting. In this paper we formulate two kinds ofiCARMs, i.e., KL-iCARM and IS-iCARM.
2.1.2 Source-Filter Factorization
One extension is obtained with the source-filter paradigm,
as listed in Table1. Kameoka and Kashino [1], for example,
originally proposed the idea of the composite autoregressive model (CARM) using fixed numbers of unconstrained
sources and autoregressive (AR) filters (all-pole transfer
functions). Although similar models were devised by someresearchers [3–5], filters were not acoustically constrained.
Badeau et al. [2] used moving-average (MA) filters (allzero transfer functions) with harmonic sources.
Some NMF variants allow sources or filters to vary overtime to richly capture temporal variations of spectral basesat the cost of increasing complexity. Heittola et al. [6] andYasuraoka and Okuno [7] used time-varying sources while
a fixed number of filters was shared over time. Hennequin
et al. [8], on the other hand, used time-varying ARMA filters that could be estimated by efficient MU rules.

To optimize the iCARMs, we propose a new class of iterative algorithms that integrates a variational Bayesian (VB)
technique with standard MU rules [8, 9].
The rest of this paper is organized as follows: Section 2
discusses the positioning of this study. Section 3 presents
the iCARMs. Section 4 describes the evaluation. Section 5
concludes the paper with a mention of future work.
2. RELATED WORKThis section introduces two machine-learning (ML) techniques, i.e., NMF and Bayesian nonparametrics.

2.1.3 Harmonicity Modeling
Another extension is based on harmonicity constraints onspectral bases or sources. For example, Vincent et al. [10]and Bertin et al. [11] assumed each basis as a weighted
sum of narrowband template spectra consisting of a few
adjacent harmonic partials. In the source-filter paradigm,
Badeau et al. [2] represented each source as a binary vector
whose elements are determined by independent Bernoulli
trials, where particular elements corresponding to harmonic
partials are more likely to take the value of 1. Yasuraoka
and Okuno [7] and Hennequin et al. [9] represented eachsource as a parametric function based on a (weighted) sum
of atomic functions (e.g., Gaussian functions) corresponding to harmonic partials. Carabias-Orti et al. [5] proposed
to further factorize a set of partials’ weights as a weighted
sum of several patterns. Efficient MU rules for estimating
the parameters of the function were proposed in [5, 9].

2.1 Nonnegative Matrix Factorization
NMF is a powerful tool for sparse decomposition of nonnegative matrix data [15]. It was first used for representing
face images as linear combinations of a compact set of basis images corresponding to “local parts” such as eyes andnoses. Such parts-based sparse representation is spontaneously induced by the nonnegativity constraint that allows
only summation of basis images. Therefore, NMF fits naturally into audio spectrogram decomposition because theenergy of harmonic sounds is concentrated at the discrete
frequencies of harmonic partials.
2.1.1 Optimization Criteria
To perform NMF, we need some criterion for evaluating
the “goodness-of-fit” of reconstructed data (linear combinations of spectral bases) to observed data (a spectrogram).

80

13th International Society for Music Information Retrieval Conference (ISMIR 2012)

M
N
I
J
Xmn
Ymn
θi
φj
Wim
Ajm
Hijn

A key feature of [7] is to consider an additional source
having a flat density. This idea was inspired by the speech
production mechanism. Excitation signals produced by vocal cords are roughly categorized into periodic signals (harmonic comb-shaped spectra) and white noise (flat spectra).
These signals are then articulated by the vocal tract whose
resonance characteristics can be represented by AR filters.
This assumption is widely accepted as reasonable to someextent for music signal modeling. In this study we model
this generative process in a Bayesian framework.

Number of frequency bins
Number of frames
Number of sources (diverges to infinity)
Number of filters (diverges to infinity)
Amplitude (power) at m-th bin and n-th frame
Reconstructed value at m-th bin and n-th frame
Global gain of i-th source
Global gain of j-th filter
Amplitude (power) of i-th source at m-th bin
Gain of j-th filter at m-th bin
Gain of i-th source & j-th filter pair at n-th frame

Table 2. Definition of mathematical symbols

2.2 Bayesian Nonparametrics
number of filters, which are assumed to go to infinity. Let
the lower-case letters m, n, i, and j indicate the indices.
In this paper we aim to factorize a nonnegative representation of X (amplitude or power spectrogram) into
three kinds of “factors” W , A, and H as follows:
I,J→∞

2
|Xmn | or |Xmn | ≈θi φj Wim Ajm Hijn
(1)

Another emerging ML technique is Bayesian nonparametrics [18], which is a generalization of the classical Bayesian
technique. In the typical Bayesian framework, we put priordistributions on unknown random variables of interests andthen, given observed data, estimate a posterior distribution
over those variables. However, this framework cannot be
used for determining model complexities (the numbers ofsources and filters in this study) because these complexities are simply treated as hyperparameters. We thus have to
use an expensive grid search for combinatorial model selection. Bayesian nonparametrics enables us to treat model
complexities as random variables and estimate their optimal values jointly with posterior computation.
Bayesian modeling is being used in music signal analysis, and Bayesian extensions of NMF [19] have been usedwith great success for audio decomposition (source separation). An especially important breakthrough was recently
made by Hoffman et al. [12]. They proposed a nonparametric Bayesian extension called the gamma-process NMF(GaP-NMF) that in theory allows an observed spectrogram
to contain an infinite number of bases. A limited effective
number of bases can be obtained by using an efficient variational inference algorithm. This extension is the basis ofa more elaborate model that can consider infinite kinds of
temporal variations of each basis [13].

i,jwhere Wim , Ajm , and Hijn respectively indicate the amplitude (power) of the i-th source at the m-th bin, the gainof the j-th filter at the m-th bin, and the gain of the i-thsource and j-th filter pair at the n-th frame. In addition,
two kinds of variables, θi and φj , are introduced to respectively indicate the global gain of the i-th source and theglobal gain of j-th filter over all N frames. Even when I
and J diverge to infinity, finite numbers of the elements ofθ and φ are expected to be substantially greater than zero
while all other elements are negligibly small. This makes it
possible for the “effective” numbers of sources and filters,
I + and J + , to be estimated in a data-driven manner.
Our goal is, given the spectrogram X, to compute a posterior distribution p(θ, φ, H|X; W , A) over random variables and estimate parameters that represent W and A. We
will discuss concrete forms of priors p(θ), p(φ), p(H),likelihood p(X|θ, φ, H; W , A), and parametric functionsof W and A according to KL or IS divergence.
3.2 Mathematical Formulation

3. PROPOSED MODELS

We explain the different formulations of iCARMs based on
KL and IS divergences.

This section presents new nonparametric Bayesian models
called infinite composite autoregressive models (iCARMs).
The essential concept of these models is inspired by a composite autoregressive model (CARM) [1] that decomposesa power spectrogram into fixed numbers of sources and AR
filters by using IS divergence as an optimization criterion.
We formulate another CARM that decomposes an amplitude spectrogram by using KL divergence. To enforce harmonicity we explicitly represent each source—except for asingle source that has a flat spectral density (white noise)—as a parametric comb-shaped function as proposed in [7].
Finally, both KL-CARM and IS-CARM are extended to intheory contain infinite numbers of sources and filters byusing gamma processes as suggested in [12].

3.2.1 Observation Likelihoods for X
We use KL oras an optimization criterion.
IS divergenceij
ij
Let Ymn be ij Ymn, where Ymn= θi φj Wim Ajm Hijn .
We aim to optimize Ymn such that the KL or IS divergence
between Xmn and Ymn is minimized, as shown in Eq.(1).
This is known to be equivalent to maximum likelihood estimation of a Poisson or exponential distribution having Ymn
as its parameter, given an observation Xmn [16]. We hereij
introduce a complex-valued latent variable Xmn
that indicates the contribution of thei-thsourceandj-thfilterpair

ij
in Xmn such that Xmn = ij Xmn
is satisfied.
The KL-iCARM is based
amplitude-additivity as on an
ij
sumption; i.e., |Xmn | = ij |Xmn|. This is obviously inij
ij
correct but is useful in practice. If |Xmn
| ∼ Poisson(Ymn),
the reproductive propertyofthePoissondistribution
leads

ij
to |Xmn | ∼ Poisson( ij Ymn), which means

3.1 Overall Framework
We first define mathematical symbols as shown in Table 2.
Let X be an M × N complex-valued spectrogram, where
M is the number of frequency bins and N is the numberof frames. Let I be the number of sources and J be the|Xmn | ∼ Poisson (Ymn )

81

(2)

13th International Society for Music Information Retrieval Conference (ISMIR 2012)

where si ≡ {sit }2M
t=1 is a signal of the i-th source andaj ≡ {aj0 , · · · , ajp }T is a coefficient vector of the j-th filteri(aj0 = 1). Let wi ≡ {wti }2M
t=1 be the autocorrelation of s
2Mand {Wim }m=1 be a complex (amplitude) spectrum density obtained by discrete Fourier transform (DFT) of w i .
ij 2MLet {Xm}m=1 be a complex spectrum density obtained by
DFT of xij . If the source signal si is a stationary Gaussianijnoise, each Xm
is independently distributed as a complex
ij
Gaussian Nc (0, Σijm ), where Σm = Wim Ajm and

The IS-iCARM is based on a complex-domain additivity
ij
ij
assumption (see Section 3.2.5). If Xmn∼ Nc (0, Ymn),
the reproductive property
ofthecomplex
Gaussianleads

ij
to Xmn ∼ Nc (0, ij Ymn). This assumption, however,may be violated when the sources are not stationary Gaussian noise (see Section 3.2.4). We nonetheless assume|Xmn |2 ∼ Exponential (Ymn )(3)

3.2.2 Gamma Process Priors on θ and φ
We put gamma process (GaP) priors on infinite-dimensional
vectors θ and φ. More specifically, we introduce independent gamma priors on elements of θ and φ as follows:
α 
γ 
θi ∼ Gamma, α , φj ∼ Gamma
,γ
(4)
I
J
As the truncation level I diverges to infinity, the vector θapproximates an infinite sequence drawn from a GaP withshape parameter α. It is proven that the effective numberof elements, I + , such that θi >  for some number  > 0
is almost surely finite. If we set I to be sufficiently larger
than α, we can expect that only a few of the I elements of θ
will be substantially greater than zero. This condensation
property enables sparse learning in an infinite space. The
same reasoning can be applied to the GaP on φ.

1
1
Ajm = 
2 = T P
aU
j −2π m pi 
m aj
j
2M
 p=0 ap e

U m is a (P + 1) × (P + 1) Toeplitz matrix with [U m ]pq =mij 2
cos(2π 2M(p − q)). This means that |Xm| is distributed
as an exponential distribution having Wim Ajm as its scaleparameter. In other words, maximum likelihood estimationof aj for xij is equivalent to minimizing the IS divergenceij 2 M2
between {|Xm| }m=1 and {Wim Ajm }M
m=1 .
In the iCARMs based on KL and IS divergences, theabove discussion leads to the following formulations:
AKL
jm =

3.2.3 Gamma Chain Priors on H
To impose smooth transitions on H, we put a gamma chainprior [20] on a temporal sequence of gains of each sourcefilter pair. More specifically, H ij is modeled as follows:

or AISjm =1
aTj U m aj

(9)

3.3 Variational and Multiplicative Optimization
The posterior over random variables p(θ, φ, H|X; W , A)
and W and A (parameters μ, σ, and a) are determinedsuch that the log-evidence log p(X; W , A) is maximized.
Since this cannot be analytically computed, we use an approximate method called variational Bayes (VB), which restricts the posterior to a factorized form given by

Gijn ∼ Gamma (β, βHijn−1 )(5)where β is a hyperparameter that controls the strength ofthe priors (degree of smoothness) and Gijn is an auxiliary
variable that imposes a positive correlation between temporally adjacent gains Hijn−1 and Hijn ( Eprior [Gijn ] =−1
Hijn−1and Eprior [Hijn ] = G−1
ijn ). Marginalizing Gijn
out, we obtain a positively correlated Markovian transition(Hijn−1 Hijn )β−1
kernel as p(Hijn |Hijn−1 ) = Γ(2β)2Γ(β) (Hijn−1 +Hijn )2β Hijn .

q(θ, φ, H) =q(θi )iq(φj )jq(Hijn )

(10)

ijn

and iteratively updates this form by monotonically increasing a lower bound3 of the log-evidence, L, given bylog p(X; W , A) ≥ E[log p(X|θ, φ, H; W , A)]+ E[log p(θ)] + E[log p(φ)] + E[log p(H)]3.2.4 Comb-shaped Functions for W− E[log q(θ)] − E[log q(φ)] − E[log q(H)] ≡ L (11)We represent each harmonic source W i as a comb-shaped
function that is the sum of H Gaussian functions, where H
is the number of harmonic partials. Specifically,H(m − hμi )2
Wim =
(6)
exp −
2σ 2

The iterative update rules areq(θ) ∝ exp(Eq(φ,H) [log p(X, θ, φ, H; W , A)])q(φ) ∝ exp(Eq(θ,H) [log p(X, θ, φ, H; W , A)])q(H) ∝ exp(Eq(θ,φ) [log p(X, θ, φ, H; W , A)]) (12)h=1

To optimize W and A (μ, σ, and a), we use multiplicative update (MU) rules inspired by [8, 9]. A general rule isobtained from the partial derivative of a “cost” function,e.g., −L. For example, if we can write the derivative as thedifference of two positive terms, i.e., −∂L∂μi = Gμi −Fμi , anwhere μi indicates F01 and σ indicates an energy diffusion
around the frequencies of partials. Note that only the lastsource is reserved as white noise, i.e., WIm = 1.
3.2.5 All-pole Transfer Functions for A

F

2M
We assume each basis signal xij ≡ {xijt }t=1 in a frame tobe represented as a P -order AR process as follows:
P

i
xij
ajp xij(7)t =−t−p + stupdate rule for μi is given by μi ← μi × Gμμi . Note that μi
i
becomes constant if the derivative is zero, and is updated in
the opposite direction of the derivative. We omit detailed
derivations and only describe update rules below.

p=1
1 When the value of F0 is given by μ̃

1
aTj U m aj

A reason for taking the “root” in the KL-iCARM is that weassume an “amplitude” spectrogram as observed data.

Hij1 ∼ Gamma (β, β/d)
Hijn ∼ Gamma (β, βGijn )(8)

2 In linear predictive coding (LPC), the source signal si is generally
limited to white noise (Wim = 1). This is a conventional assumption.
3 More specifically, a further lower bound of L should be computed.

i [Hz], μi = μ̃i /(r/2M ) [bins],

where r is a sampling rate and 2M is a window size of frequency analysis.

82

13th International Society for Music Information Retrieval Conference (ISMIR 2012)

3.3.1 Variational Updates for KL-iCARM

4. EVALUATION

The variational posterior of each random variable is set tobe the same family as its prior distribution as follows:q(θi ) = Gamma(aθi , bθi ), q(φj ) = Gamma(aφj , bφj )Hq(Hijn ) = Gamma(aHijn , bijn )We report comparative experiments that were conducted toevaluate the performance of the iCARMs based on KL and
IS divergences as multipitch analyzers.

(13)

4.1 Experimental Conditions

The variational parameters are given byα aθi = + mnj |Xmn |λmnijI bθi = α + mnj E[φj Wim Ajm Hijn ]γ 
aφj = + mni |Xmn |λmnij
J 
bφj = γ + mni E[θi Wim Ajm Hijn ]

aH
(14)
ijn = 2β +m |Xmn |λmnijHbijn = βE[Gijn +Gijn+1 ] + m E[θi φj Wim Ajm ]where λmnij is an auxiliary variable given by
λmnij ∝ exp(E[log(θi φj Wim Ajm Hijn )])We used thirty pieces of “ENSTDkCl” subset included in
the MAPS piano database [21]. We truncated each piece to
30 s as in [5,11] and converted the original CD-quality signals into monaural signals sampled at 16 [kHz]. The spectrograms were obtained with short-time Fourier transform(STFT) with a window size of 2048 samples and a shifting
interval of 10 [ms], i.e., M = 1024 and N = 3000. The
amplitude or power spectrogram of each piece was scaled
such that M1N mn |Xmn | = 1 or maxmn |Xmn |2 = 1.
The hyperparameters were specified as I = 88+1, J = 10,
α = 1, β = γ = 0.1, H = 20, P = 4, and d = Eemp [|Xmn |]or Eemp [|Xmn |2 ]. Although J = 10 was too small to accurately approximate the GaP, it was sufficiently large in our
experiments because the audio signals contain only pianosounds. We initialized {μi }88
i=1 as the frequencies corresponding to the 88 keys of the standard piano. The otherparameters were initialized randomly.
Multiple F0s were detected at each frame in a thresholding process. If the gain of the i-th source, j θi φj Hijn ,was larger than the threshold, we judged that the n-th frame
includes an F0 indicated by μi . The threshold was globally
determined such that the frame-level precision and recall
rates were balanced to yield the best average F-measure.

(15)

3.3.2 Variational Updates for IS-iCARM
As proposed in [12], the variational posterior of each variable is given by a generalized inverse-Gaussian (GIG) distribution (see the Appendix) as follows:q(θi ) = GIG(aθi , bθi , cθi ), q(φj ) = GIG(aφj , bφj , cφj )HHq(Hijn ) = GIG(aHijn , bijn , cijn )(16)

The variational parameters are given byαE[φj Wim Ajm Hijn ]aθi = , bθi = α + mnj
ξmn
I

θ
2 2
ci = mnj |Xmn | ηmnij E φj Wim A1jm Hijn
γ
E[θ W A H]
aφj = , bφj = γ + mni i imξmnjm ijn
J2cφj = mni |Xmn |2 ηmnij
E θi Wim A1jm Hijn

H
1
2 2
aH
ijn = 2β, cijn =m |Xmn | ηmnij E θi φj Wim Ajm
 E[θi φj Wim Ajm ]bH(17)
ijn = βE[Gijn +Gijn+1 ] +m
ξmn

4.2 Experimental Results
We first tested our models on toy data obtained by extracting the first 4.9 s (490 frames) of the piece “alb se2,” whichcontains five different F0s and a polyphony level that increases one by one up to five (D4, +C#4, +C4, +A3, +F#3).
As shown in Fig. 2, the KL-iCARM could successfully discover the correct number of sources (five harmonic sources+ one white-noise source) in a data-driven manner. In addition, we could separate X into harmonicand noise compo
inents by computing E[Ymn] = j E[θi φj Wim Ajm Hijn ],which represents the component of the i-th source at them-th bin and n-th frame.
As shown in Fig. 3, on the other hand, the IS-iCARM
overestimated the numbers of sources and filters and mademany octave errors (half-F0 errors). One reason is that IS
divergence permits a reconstructed power to exceed an observed power with a smaller penalty. It is therefore difficult
to reduce false alarms of harmonic partials.
We then used the 30 pieces for evaluation. The KL- and
IS-iCARMs achieved the frame-level F-measures of 48.4%
and 35.1% respectively. Although these preliminary results
were not really impressive compared with the state-of-theart results [5,11], we consider our framework to be promising because of its elegant nature of sparse learning over an
infinite space. A main limitation of the KL-iCARM is that
we still need to resort a thresholding process for temporalgains although limited numbers of sources and filters canbe obtained by using GaPs. One solution would be to introduce binary latent variables that indicate note existences.

where ηmnij and ξmn are auxiliary variables given byηmnij ∝ E[ θi φj Wim1Ajm Hijn ]−1 s.t. ij ηmnij = 1
ξmn = ij E[θi φj Wim Ajm Hijn ](18)
3.3.3 Multiplicative Updates for KL- and IS-iCARMs
The MU rules for μ, σ, and a are given by μi ← G−1
μ i F μ i μi ,−122
σ 2 ← G−1
F
σ
,
and
a
←
G
F
a
,
where
j
aj j
aj
σ2 σ
2

i)
F
G
Fμi = mnjh h(mVmnij+ hμi Vmnij) exp − (m−hμ
2σ 2
2

i)
G
F
Gμi = mnjh h(mVmnij+ hμi Vmnij) exp − (m−hμ
2σ 2
2

i)
F
Fσ2 = mnijh Vmnij(m − hμi )2 exp − (m−hμ
2σ 2
2

i)
G
Gσ2 = mnijh Vmnij(m − hμi )2 exp − (m−hμ
2σ 2

3
F KL
aj =
mni θi φj Wim Hijn Ajm U m

2
GKL
aj =
mni |Xmn |λmnij Ajm U mE[θi φj Wim Hijn ] 2
F IS
Ajm U m
aj =
mni
ξmn

IS
2 2
Gaj = mni |Xmn | ηmnij E θi φj W1im Hijn U m
(19)
F
G
F
and Vmnij
are given by Vmnij= E[θi φj Ajm Hijn ]
Vmnij
−1
G
and Vmnij = |Xmn |λmnij Wim in the KL-iCARM. On
E[θ φ AH]

i j jm ijn
FG=
and Vmnij=
the other hand, Vmnij
ξmn
1
2 2
|Xmn | ηmnij E θi φj W 2 Ajm Hijn in the IS-iCARM.
im

83

13th International Society for Music Information Retrieval Conference (ISMIR 2012)

Observation X

Reconstruction Y

Observation X

Reconstruction Y

Harmonic components in Y

Noise components in Y

Harmonic components in Y

Noise components in Y

Global gains of sources E[θ]Effective AR filters in A

Global gains of sources E[θ]Effective AR filters in A

C#4
D4

F#3

Overestimated sources
with heavy-tailed gains

Unnecessary sources with

A3 C4 sufficienly small gains

Global gains of filters E[φ]Global gains of filters E[φ]Unnecessary filters withvanishingly small gains

Overestimated filterswith heavy-tailed gains

Figure 2. Decomposition results obtained by KL-iCARM

Figure 3. Decomposition results obtained by IS-iCARM

5. CONCLUSION[8] R. Hennequin, R. Badeau, and B. David. NMF with timefrequency activations to model nonstationary audio events.
IEEE Trans. on ASLP, 19(4):744–753, 2011.
[9] R. Hennequin, R. Badeau, and B. David. Time-dependent
parametric and harmonic templates in non-negative matrix
factorization. DAFx, pp. 1–8, 2010.
[10] E. Vincent, N. Bertin, and R. Badeau. Adaptive harmonic
spectral decomposition for multiple pitch estimation. IEEE
Trans. on ASLP, 18(3):528–537, 2010.
[11] N. Bertin, R. Badeau, and E. Vincent. Enforcing harmonicity
and smoothness in Bayesian non-negative matrix factorization applied to polyphonic music transcription. IEEE Trans.
on ASLP, 18(3):538–549, 2010.
[12] M. Hoffman, D. Blei, and P. Cook. Bayesian nonparametricmatrix factorization for recorded music. ICML, 2010.
[13] M. Nakano et al. Bayesian nonparametric spectrogram modeling based on infinite factorial infinite hidden Markov
model. WASPAA, pp. 325–328, 2011.
[14] B. Niedermayer. Improving accuracy of polyphonic musicto-score alignment. ISMIR, pp. 585–590, 2009.
[15] D. Lee and H. Seung. Algorithms for non-negative matrix
factorization. NIPS, pp. 556–562, 2000.
[16] C. Févotte, N. Bertin, and J.-L. Durrieu. Nonnegative matrix
factorization with the Itakura-Saito divergence: With application to music analysis. NECO, 21(3):793–830, 2009.
[17] C. Févotte and J. Idier. Algorithms for nonnegative matrix
factorization with the beta-divergence. NECO, 23(9):2421–
2456, 2011.
[18] P. Orbanz and Y. W. Teh. Bayesian nonparametric models.
Encyclopedia of Machine Learning. Springer, 2010.
[19] A. T. Cemgil. Bayesian inference for nonnegative matrix factorisation models. Computational Intelligence and Neuroscience, 2009:Article ID 785152, 2009.
[20] A. T. Cemgil and O. Dikmen. Conjugate gamma Markov random
fields for modelling nonstationary sources. ICA, 2007.
[21] V. Emiya, R. Badeau, and B. David. Multipitch estimation of
piano sounds using a new probabilistic spectral smoothness
principle. IEEE Trans. on ASLP, 18(6):1643–1654, 2010.

We presented nonparametric Bayesian models called infinite composite autoregressive models (iCARMs) that decompose an observed spectrogram into three kinds of factors, i.e., sources, filters, and time-varying gains of sourcefilter pairs. The experimental results showed that appropriate numbers of sources and filters can be discovered in adata-driven manner by using gamma processes for sparselearning. To improve the accuracy of multipitch analysis,
we are considering the use of log-frequency spectrograms
obtained by constant-Q or wavelet transform. We also plan
to use these models for “timbre-based” source separation
by distinguishing different resonance characteristics of instrument and vocal sounds by AR filters.
Acknowledgment: This study was supported by JSPS
KAKENHI 23700184 and JST OngaCREST project.
6. REFERENCES[1] H. Kameoka and K. Kashino. Composite autoregressive system for sparse source-filter representation of speech. ICASSP,
pp. 2477–2480, 2009.
[2] R. Badeau, V. Emiya, and B. David. Expectation-maximization
algorithm for multi-pitch estimation and separation of overlapping harmonic spectra. ICASSP, pp. 3073–3076, 2009.
[3] J.-L. Durrieu, G. Richard, B. David, and C. Févotte. Source/Filter
model for unsupervised main melody extraction from polyphonicaudio signals. IEEE Trans. on ASLP, 18(3):564–575, 2010.
[4] T. Virtanen and A. Klapuri. Analysis of polyphonic audio
using source-filter model and non-negative matrix factorization. NIPS Ws. on Adv. in Mod. for Acoust. Proc., 2009.
[5] J. J. Carabias-Orti, T. Virtanen, P. Vera-Candeas, N. Ruiz-Reyes,
and F. J. Cañadas-Quesada. Musical instrument sound multiexcitation model for non-negative spectrogram factorization.
IEEE J. of Sel. Top. in Sig. Proc., 5(6):1144–1158, 2011.
[6] T. Heittola, A. Klapuri, and T. Virtanen. Musical instrument
recognition in polyphonic audio using source-filter model forsound separation. ISMIR, pp. 327–332, 2009.
[7] N. Yasuraoka and H. G. Okuno. Musical audio signal modeling for joint estimation of harmonic, inharmonic, and timbral
structure and its application to source sepatation. SIG Technical Reports, volume 2012-MUS-94, pp. 1–8, 2012.

A. PROBABILITY DISTRIBUTIONSx
1ba a−1 −bx
eExponential(x|λ) = e− λGamma(x|a, b) =
x
Γ(a)
λ
a

Poisson(x|λ) =

84

λx −λ
(b/c) 2 xa−1 −(bx+ c )x√e
GIG(x|a, b, c) =ex!
2Ka (2 bc)