INFINITE SUPERIMPOSED DISCRETE ALL-POLE MODELINGFOR MULTIPITCH ANALYSIS OF WAVELET SPECTROGRAMS
Kazuyoshi Yoshii1
Katsutoshi Itoyama1
Masataka Goto2
1
Graduate School of Informatics, Kyoto University, Japan
2
National Institute of Advanced Industrial Science and Technology (AIST), Japan{yoshii,itoyama}@kuis.kyoto-u.ac.jp m.goto@aist.go.jp
ABSTRACT

Statistical modeling of music audio signals based on machine learning techniques is a hot topic in the field of musicsignal analysis. In particular, nonnegative matrix factorization (NMF) has often been used for multiple fundamental
frequency (F0) estimation (multipitch analysis) and source
separation [1–4, 7, 14, 15, 17, 21–26]. The standard NMF
approximates a nonnegative spectrogram (matrix) as the
product of two nonnegative matrices: a set of basis spectra and a set of the corresponding activations. An efficient
multiplicative-updating (MU) algorithm was proposed forminimizing a cost function that measures the approximation error [18]. This was later found to be maximum likelihood estimation of a particular probabilistic model [5].
Statistical source-filter models, which were inspired by
the simplified model of the speech production mechanism,c Kazuyoshi Yoshii, Katsutoshi Itoyama, Masataka Goto.
Licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0). Attribution: Kazuyoshi Yoshii, Katsutoshi Itoyama,
Masataka Goto. “Infinite Superimposed Discrete All-pole Modeling forMultipitch Analysis of Wavelet Spectrograms”, 16th International Society for Music Information Retrieval Conference, 2015.

86

bases

⊗
gains
Log frequency [cents]

Log frequency [cents]⊕…

1. INTRODUCTION⊗…

This paper presents a statistical multipich analyzer basedon a source-filter model that decomposes a target music
audio signal in terms of three major kinds of sound quantities: pitch (fundamental frequency: F0), timbre (spectralenvelope), and intensity (amplitude). If the spectral envelope of an isolated sound is represented by an all-pole filter,
linear predictive coding (LPC) can be used for filter estimation in the linear-frequency domain. The main problemof LPC is that although only the amplitudes of harmonic
partials are reliable samples drawn from the spectral envelope, the whole spectrum is used for filter estimation. To
solve this problem, we propose an infinite superimposed
discrete all-pole (iSDAP) model that, given a music signal,
can estimate an appropriate number of superimposed harmonic structures whose harmonic partials are drawn froma limited number of spectral envelopes. Our nonparametric Bayesian source-filter model is formulated in the logfrequency domain that better suits the frequency characteristics of human audition. Experimental results showed thatthe proposed model outperformed the counterpart model
formulated in the linear frequency domain.

Only the weights of harmonic partials are scaled according to a filter
DAP
Harmonic structures
All-pole transfer functions

sources (pitches)filters (timbres)
Observed spectrum of each frame

Figure 1. Overview of infinite superimposed discrete allpole (iSDAP) modeling: We take the infinite limit as boththe numbers of sources and filters, I and J, go to infinity.
have often been proposed for representing musical instrument sounds [7,14,25]. The pitches and timbres of musicalinstrument sounds are well characterized by fine structures(sources) and spectral envelopes (filters) in the frequency
domain. Since the human auditory system is sensitive tospectral peaks and formants, the spectral envelope of eachframe is usually modeled by an all-pole frequency transfer
function (frequency response of an autoregressive (AR) filter) [14]. A classical method of all-pole spectral envelope
estimation called linear predictive coding (LPC) [16] corresponds to maximum likelihood estimation of a particularprobabilistic model under a strong assumption that source
signals have the flat spectrum (white noise).
The composite autoregressive (CAR) modeling [17] isa promising statistical approach that overcomes the limitation of classical source-filter modeling in the frameworkof NMF. A given audio spectrogram is decomposed into
specified numbers of fine structures (sources) and spectralenvelopes (filters). A key feature of this approach is thatsource spectra themselves can be estimated (not limited
to white noise) at the same time as all-pole spectral envelope estimation. The probabilistic interpretation of sourcefilter NMF makes it possible to formulate a nonparametric Bayesian extension called infinite CAR (iCAR) modeling that can automatically choose the appropriate numbersof sources and filters according to a given audio spectrogram [26]. Another useful extension is to restrict source
spectra to harmonic structures by using parametric functions [26]. The F0s of source spectra can be estimated in a
principled maximum-likelihood framework.
Conventional methods of source-filter NMF including
CAR [7,14,17,25,26] have two major problems as follows:Proceedings of the 16th ISMIR Conference, Málaga, Spain, October 26-30, 2015
1. All the frequency bins are taken into account forspectral envelope estimation although only the amplitudes of harmonic partials can be regarded as reliable samples from spectral envelopes.
2. Linear-frequency spectrograms given by short-time
Fourier transform (STFT) are used for all-pole modeling although log-frequency spectrograms given by
wavelet or constant-Q transform better suit the frequency characteristics of human audition.
To solve these problems, we propose a new statistical
approach to source-filter NMF called infinite superimposed
discrete-all pole (iSDAP) modeling. Our approach is basedon a well-known technique called discrete all-pole (DAP)
modeling [8] that takes into account only the peaks of harmonic partials for spectral envelope estimation. To deal
with polyphonic music audio signals, however, we need toseparate individual harmonic structures and estimate theirF0s (positions of discrete harmonic partials). A major contribution of this study is to extend the DAP modeling for
dealing with an arbitrary number of superimposed harmonic
structures in a similar way to the iCAR modeling. This enables us to decompose a log-frequency spectrogram into
appropriate numbers of pitches (F0s), timbres (spectral envelopes), and their volumes by leveraging the frequencyscale-free characteristics of the DAP modeling.

20
0

[dB]

2.1 Linear Predictive Coding (All-pole Modeling)
The linear predictive coding (LPC) [16] is a signal modeling method that can be used for estimating the spectral envelope of an observed spectrum. The underlying assumption is that the corresponding audio signal x = {xm }1m= 1(a local signal {xm }M
m=1 is infinitely repeated) follows a
P -order autoregressive (AR) process as follows:
xm =

P
X

ap xm p + sm i.e.,

p=1

P
X

ap xm p = sm , (1)p=0where a = [a0 , · · · , aP ]T is a vector of AR coefficients(a0 = 1) and {sm }M
m=1 is a set of prediction errors. Eq. (1)can be interpreted in terms of source-filter modeling, i.e.,
when x is a speech signal, s is an excitation signal generated by the vocal cords (source) and a represents the resonance characteristics of the vocal tract (filter).
Eq. (1) can be regarded as a linear system (governed bya) that takes s as input and then gives x as output. Since
Eq. (1) is a convolution of a with x, we can say
A(z)X(z) = S(z) i.e., X(z) = S(z)F (z),(2)

where X(z) and S(z) are the z-transforms of x and s, respectively, which are given by
X(z) =

1X

m= 1

xm z m and S(z) =

1
X

m= 1

sm z m , (3)DAP is not biased towards harmonic partials
Discrete all-pole modeling (DAP)

-20
-40
-60
-80
-100
-120

Linear predictive coding (LPC)

Observed spectrum(harmonic structure)

-140

1

2

3

4

5

6

7

8 [kHz]Figure 2. Spectral envelopes estimated by LPC and DAP.
def

1
is an all-pole transfer function given by
and F (z) = A(Z)F (z) =

1
1
= PP
.
p
A(z)
p=0 ap z(4)m
Letting 2⇡ M
= !m and substituting z = ei!m into Eq. (2),we get the Fourier-transform representation as follows:

X(ei!m ) = S(ei!m )F (ei!m ),(5)

where {X(ei!m )}M
m=1 is the complex spectrum of the observed signal x, {S(ei!m )}M
m=1 is that of the source signal
isthefrequency characteristics ofs, and {F (ei!m )}M
m=1
the all-pole transfer function.
The goal of LPC is to estimate a set of AR coefficientsa under a strong unrealistic assumption that the source signal s is Gaussian white noise. This means that S(ei!m ) iscomplex Gaussian-distributed as follows:

2. RELATED WORKThis section reviews probabilistic models of source-filter
decomposition, NMF, and source-filter NMF as a basis offormulating the iSDAP model. Most conventional models
are formulated in the linear frequency (STFT) domain.

87

S(ei!m ) ⇠ Nc (0,

2

(6)

),

where
is the power of the white spectrum of the source
signal s. Using Eq. (5) and Eq. (6), we get2

X(ei!m ) ⇠ Nc (0,2|F (ei!m )|2 ).

(7)

Letting Xm = |X(e )| and Fm = |F (ebriefly rewrite Eq. (7) as follows:
i!m

2

i!m

Xm ⇠ Exponential( 2 Fm ),)| , we2(8)

where {Xm }M
m=1 is the power spectrum of the observedMsignal x and {Fm }M
m=1 is the spectral envelope of {Xm }m=1 ,as shown in Figure 2. Eq. (8) defines the probabilistic model2of LPC. {Fm }M
can be estimated in a
m=1 (i.e., a) andmaximum-likelihood manner [16].
The main problem of LPC is that if we analyze a pitchedsound derived from a periodic source signal (e.g., vibrationof strings), the estimated envelope {Fm }M
m=1 loosely fits
and
its
peaks
(formants)
the observed spectrum {Xm }M
m=1
tend to be biased to the positions of harmonic partials. Thisis because all M frequency bins are used for all-pole modeling although in reality only the amplitudes of harmonic
partials can be considered to be reliable samples from thespectral envelope.
2.2 Discrete All-pole Modeling
The discrete all-pole (DAP) modeling [8] is a well-known
spectral envelope estimation method that was proposed forsolving the problem of LPC. Since DAP is an extension of
LPC, the probabilistic model of DAP has the same form
as Eq. (8). A key feature of DAP is that Eq. (8) is definedover only a partial set of frequency bins, ⌦, as follows:
Xm ⇠ Exponential( 2 Fm ) m 2 ⌦,

(9)

88

Proceedings of the 16th ISMIR Conference, Málaga, Spain, October 26-30, 2015

The composite autoregressive (CAR) modeling [17] is avariant of source-filter NMF that is used for decomposing alinear-frequency mixture spectrogram into I fine structures(sources) and J spectral envelopes (filters), as shown in
Figure 3. Let X be an M ⇥ N power spectrogram, where
M is the number of frequency bins and N is the number offrames. The nonnegative matrix X is factorized into threekinds of “factors” S, F , and H as follows:

⊗
bases

⊗
gains
Linear frequency [Hz]Linear frequency [Hz]⊕…

2.3 Composite Autoregressive Modeling

The whole harmonic structure is scaled according to a filter
LPC
Harmonic structures
All-pole transfer functions…

where if ⌦ = {1, · · · , M }, DAP reduces to LPC. To estimate the spectral envelope of a harmonic spectrum, we cantake into account only the discrete frequencies of harmonic
partials. The estimated envelope passes close to the peaksof harmonic partials (Figure 2). To maximize the likelihood given by Eq. (9), an efficient algorithm was proposed
for alternately optimizing a and 2 [8]. It was later foundas a multiplicative updating algorithm [1, 14].
The main limitation of DAP is that the F0 and its overtones of an observed spectrum {Xm }M
m=1 should be given
in advance for defining a set of discrete frequencies to beconsidered, ⌦. To analyze polyphonic music audio signals
consisting of superimposed harmonic structures, we needto separate harmonic structures and estimate their F0s.

sources (pitches)filters (timbres)
Observed spectrum of each frame

Figure 3. Overview of composite autoregressive (CAR)
modeling defined in the linear frequency domain.
2.3.2 Several Extensions
Another probabilistic model of CAR was proposed by using the Kullback-Leibler (KL) divergence instead of the IS
divergence as a cost function for a practical reason [26].
Instead of Eq. (14), we assume
Xmn ⇠ Poisson (Ymn ) ,(15)

M
where E[Xmn ] = Ymn holds. {Xmn }M
m=1 and {Ymn }m=1

The probabilistic model of CAR can be formulated by precisely modeling source signals in a statistical manner. To
avoid the unrealistic assumption of LPC that each source
signal is Gaussian white noise (Eq. (6)), we assume

are the amplitude spectra of frame n because KL-NMF
models are usually formulated in the amplitude domain byassuming the amplitude additivity [10, 18].
A nonparametric Bayesian extension called infinite CAR
enables us to automatically estimate appropriate numbersof sources and filters according to the observation X [26].
This technique is based on gamma process NMF [15].
Another extension of CAR is to force the amplitude
spectrum of each source {Sim }M
m=1 to have a harmonic
structure [26]. If the source signal is a train of periodic impulses (an idealized model of the vocal chords), {Sim }M
m=1
has a harmonic structure consisting of equally-spaced harmonic partials with the same weight. The optimal valueof the F0 can be estimated such that the likelihood given
by Eq. (15) is maximized. This technique of F0 estimation
has a potential to solve the limitation of DAP.

Si (ei!m ) ⇠ Nc (0, Sim ) ,3. PROPOSED MODEL

Xmn ⇡

I X
J
X

def

Sim Fjm Hnij = Ymn ,(10)

i=1 j=1

where {Sim }M
m=1 is the linear-frequency power spectrumof source i, {Fjm }M
m=1 is that of filter j, and Hnij is thegain of a pair of source i and filter j at frame n. All thesevariables should be estimated from X.
2.3.1 Original Formulation(11)

where {Si (ei. Using Eq. (5) and Eq. (11), we geti!m)}M
m=1 is the complex spectrum of source

Xijmn (ei!m ) ⇠ Nc (0, Sim Fjm Hnij ) ,(12)

Xmn (ei!m ) ⇠ Nc (0, Ymn ) ,(13)

Xmn ⇠ Exponential (Ymn ) ,(14)

where {Xijmn (egenerated from source i and filter j at frame n. Using the
reproducing property of the Gaussian and Eq. (10), we geti!m)}M
m=1 is a latent complex spectrum

where {Xmn (etrum at frame n. Eq. (13) is equivalent toi!m)}M
m=1 is the observed complex spec-

where E[Xmn ] = Ymn is satisfied and {Xmn }M
m=1 and{Ymn }M
m=1 are the power spectra of frame n.

This means that the Itakura-Saito (IS) divergence is theoretically justified as a cost function that evaluates the error
between Xm and Ym in Eq. (10) [17]. In general, however,
optimization algorithms tend to get stuck in bad local minima because the IS divergence is not convex w.r.t. Ymn .

This section presents a nonparametric Bayesian approachcalled infinite superimposed discrete all-pole (iSDAP) modeling for source-filter decomposition of wavelet spectrograms. Our model can estimate multiple F0s at each frameand discover several kinds of instrument timbres (all-pole
spectral envelopes) from polyphonic music audio signals.
To achieve this, we integrate the technique of discrete allpole (DAP) modeling [8] into the framework of composite
autoregressive (CAR) modeling [17, 26] in a probabilisticmanner. The iSDAP model can be regarded as a Bayesian
extension of log-frequency source-filter NMF based on asingle filter [19], and has all of the following features:
1. Superimposed DAP modeling: Our model can estimate the spectral envelope of each of harmonic structure contained in mixed sounds. The original DAP
model can deal with only isolated sounds [8].
2. Precise F0 modeling: Each frame is allowed to contain a unique set of F0s (sources) for capturing fineProceedings of the 16th ISMIR Conference, Málaga, Spain, October 26-30, 2015

3. Log-frequency modeling: Our source-filter model
can deal with wavelet spectrograms that suit the characteristics of human audition by leveraging an advantage of DAP modeling that only discrete frequencies are required for spectral envelope estimation.
4. Bayesian nonparametrics: Our model can estimate
effective numbers of sources and filters according toa given spectrogram by allowing unbounded (infinite
in theory) numbers of sources and filters to be used.
3.1 Model Formulation
We explain a probabilistic model of iSDAP. Let X be an
M ⇥ N log-frequency amplitude spectrogram with M frequency bins and N frames. The nonnegative matrix X isfactorized in a similar way to Eq. (10) as follows:
1
0
I!1
X J!1
X
Xmn ⇠ Poisson@✓ni j Wnijm Hnij A, (16)
i=1

j=1

where ✓ni is the local weight of source i at frame n, j
is the global weight of filter j, and Hnij is the gain of apair of source i and filter j at frame n. {Wnijm }M
m=1 isthe amplitude spectrum derived from the source-filter pairat frame n. Note that ✓ni and Wnijm are allowed to varyover time to represent the F0 fluctuation unlike Eq. (10).
We aim to perform sparse learning of weight vectors ✓ n =[✓n1 , · · · , ✓nI ]T and = [ 1 , · · · , J ]T when the numberof sources I and the number of filters J go to infinity.

The all-pole spectral envelope determines the weights of harmonic partials
Amplitude

fluctuations of F0s (e.g., vibrato). The original CAR
models [17, 26] assume that a common set of source
spectra (semitone-level F0s) is shared over all frames.

89

・・・

・・・

Log-frequency [cents]

Figure 4. Composition of source i and filter j at frame nin the log-frequency domain.
(KL-NMF) is considered to fit the amplitude domain rather
than the power domain [19].
3.1.2 Prior Distributions
We put gamma process (GaP) priors on infinite-dimensional
vectors ✓n and as in [15, 26]. Specifically, we put independent gamma priors on elements of ✓n and as follows:
⇣↵
⇣↵
⌘
⌘
✓
, ↵✓ , j ⇠ Gamma, ↵ , (20)✓ni ⇠ Gamma
I
J
where ↵✓ and ↵ are hyperparameters called concentration
parameters. As J diverges to infinity, the vector approximates an infinite vector drawn from a GaP with ↵ . It is
proven that the effective number of filters, J + , such thatj > ✏ for some number ✏ > 0 is almost surely finite [15].
If J is sufficiently larger than ↵ (J is often called a truncation level in weak-limit approximation to infinite modeling), the GaP can be well approximated. The same reasoning can be applied to the GaP on ✓n . On the other hand,
we put a standard Gamma prior on Hnij as follows:
Hnij ⇠ Gamma(aH , bH ),(21)

where aH and bH are shape and rate hyperparameters.

3.1.1 Parametric Functions

3.2 Variational Inference

As shown in Figure 4, we force the amplitude spectrum{Wnijm }M
m=1 to have a harmonic structure as follows:

The posterior over random variables p(✓, , H|X;µ,a)and parameters µ and a are determined such that a lower
bound L of the log-evidence log p(X;µ,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 byYYYq(✓, , H) =q(✓ni )q( j )
q(Hnij ). (22)

Wnijm =R
X

Smnir Fnijr ,

(17)

r=1

where R is the number of harmonic partials and {Smnir }M
m=1
is the monomodal spectrum of the r-th harmonic partial ofsource i at frame n given by
✓
◆
1
2
(f
(µ
+
1200
log
r))
, (18)
Smnir = exp
m
ni
2
2 2

where µni is the F0 [cents] of source i at frame n, fm is thelog-frequency [cents] corresponding to the m-th bin, and2
indicates energy diffusion around harmonic partials.
We then represent the weights of discrete harmonic partials, {Fnijr }R
r=1 , by using an all-pole transfer function in
the log frequency domain as follows:1
1= aTj U (!nir )aj 2, (19)Fnijr = PP!nir pi
p=0 ajp ewhere aj ⌘ [aj0 , · · · , ajP ]T , !nir is a normalized frequency [rad] corresponding to the r-th harmonic partial ofsource i at frame n, and U (!) is a (P +1)⇥(P +1) matrix
with [U (!)]pq = cos(!(p q)). Note that Fnijr indicatesthe value of amplitude (not power). The Poisson likelihood

ni

j

nij

Iteratively updating this posterior, we can monotonicallyincrease a lower bound of the log-evidence given bylog p(X; µ, a)E[log p(X|✓, , H; µ, a)]+ E[log p(✓)] + E[log p( )] + E[log p(H)]E[log q(✓)]

E[log q( )]E[log q(H)] ⌘ L0 , (23)where the first term can be further lower bounded by Jensen’s
inequality on the concave logarithmic function as follows:E[log p(X|✓, , H; µ, a)]h PiP= mn Xmn E log ijr ✓ni j Smnir Fnijr Hnij⇤⇥P
mnijr E ✓ni j Smnir Fnijr Hnij + const.
P

mnijr

h

mnijr Xmn E log✓nij Smnir Fnijr Hnij
mnijr

i

⇤
⇥
mnijr E ✓ni j Smnir Fnijr Hnij + const. (24)P

90

Proceedings of the 16th ISMIR Conference, Málaga, Spain, October 26-30, 2015
where
P
ijr

mnijr is a normalized auxiliary variable such thatmnijr = 1. The equality holds (i.e., the lower boundof L0 is maximized) if and only if

mnijr / exp(E[log(✓ni j Smnir Fnijr Hnij )]).

(25)

Using Eq. (24), the objective function of our model to be
maximized, L, is obtained as the lower bound of L0 . For
convenience, we define Xmnijr and Ymnijr as
mnijr Xmn ,(26)

Ymnijr = E[✓ni j Smnir Fnijr Hnij ].

(27)

Xmnijr =As in [1, 14], we then represent the partial derivative of@L= (Rj R0j )aj , where Rj
L with respect to aj as @ aj0
and Rj are positive definite matrices given by
X
2
Xmnijr Fnijr
U (!nir ),
(33)
Rj =
mnir

R0j =

X

2
Ymnijr Fnijr
U (!nir ).

(34)

mnir

The updating rule of aj is given byaj

Rj 1 R0j aj .

(35)

Finally, we forcibly adjust the scale of the filter Fnijr suchthat ↵j0 = 1 for normalizing the filter. Although this step
violates the convergence of the optimization algorithm, it
was empirically found to work well.

3.3 Variational Bayesian Updating of ✓, , and H
The VB updating rules are given byq(✓) / exp(Eq( ,H) [log p(X, ✓, , H; µ, a)]),q( ) / exp(Eq(✓,H) [log p(X, ✓, , H; µ, a)]), (28)

3.5 Binary Piano-roll Estimation

q(H) / exp(Eq(✓, ) [log p(X, ✓, , H; µ, a)]).

To perform multipitch analysis, i.e., make a binary pianoroll representation, we need to judge the existence of each
The variational posterior of each random variable is set tosemitone-level pitch at each frame. Using a trained model,
be the same family as its prior distribution as follows:
we calculate an activation matrix V = {Vkn }88,N
k=1,n=1 overq(✓ni ) = Gamma(a✓ni , b✓ni ), q( j ) = Gamma(aj , bj ),pitch k and frame n (continuous-valued piano-roll repreH
q(Hnij ) = Gamma(aH(29)
sentation e.g., the middle figure of Figure 5) by accumunij , bnij ).
latingP the expected amplitude of the first partial of source
The variational parameters are given byi, j E[✓ni j Fnij1 Hnij ], into Vkn indicated by µni . FiX
X↵✓a✓ni =
Xmnijr , b✓ni = ↵✓ +E[ j Hnij ]Wnijm , nally, the activation matrix V P+
is normalized such that all
I
mjr
mjr
the elements sum to unity, i.e., kn Vkn = 1.
X
X
↵
+
aj =
Xmnijr , bj = ↵ +E[✓ni Hnij ]Wnijm , There are several approaches to binary piano-roll estiJ mnirmation. The common approach is to make a binary decimnir
X
X
sion based on a threshold ⌘. Another approach is to defineaH
Xmnijr , bH
E[✓ni j ]Wnijm . a hidden Markov model (HMM) and use the Viterbi-searchnij = aH +nij = bH +mrmralgorithm for estimating a sequence of hidden binary states
To estimate the effective number of filters J + , we per{Zkn }N
n=1 from a sequence of pitch-existence likelihoods
p Nform sparse learning. If E[ j ] becomes sufficiently small{Vkn}n=1 for each pitch k, where p controls the dynamic
for some filter j, we degenerate it and J
J 1. A simirange. In our implementation, p = 0.2 and the transition
lar treatment is applied to E[✓ni ]. The proposed variational
matrix is [0.8, 0.2; 0.01, 0.99] in the Matlab notation.
algorithm is gradually accelerated per iteration.
4. EVALUATION

3.4 Multiplicative Updating of µ and a
To estimate parameters µ and a, we use the multiplicative
update (MU) algorithm as in [1,14]. In general, to optimize
x, we represent the partial derivative of a “cost” function
C with respect to x as the difference of two positive terms,
i.e., @C
R0 . An updating rule of x is then given by@x = R
R0x
R x. Note that x becomes constant if the derivative iszero, and is updated in the opposite direction of the derivative. In this study the cost function is the negative lower
bound of the log-evidence, L.
First, we represent the partial derivative of L with re00spect to µmi as @µ@L
= Rni Rni
, where Rni and Rni
areni
positive terms given byXRni =(µni + 1200 log2 r)Xmnijr + fm Ymnijr , (30)
mjr

0
Rni
=

X

fm Xmnijr + (µni + 1200 log2 r)Ymnijr ,(31)mjr

The updating rule of µni is given byµni

0
Rni1 Rni
µni .

(32)

We report comparative experiments that were conducted
for evaluating the performance of the iSDAP model in multipitch analysis of piano music. Since the proposed model
assumes that input mixture signals contain only harmonicsounds, we also tested the use of harmonic and percussive
source separation (HPSS) [12] as a preprocessor.
4.1 Experimental Conditions
We used 30 pieces (labeled as ”ENSTDkCl”) selected from
the MAPS database [9] that contain stereo signals sampled at 44.1 [kHz]. The audio signals were converted tomonaural signals and truncated to 30 [s] from the beginning as in [2, 4, 21, 22, 24]. The amplitude spectrogram ofeach piece over the frequency bins ranging from 0 [cents](16.325 [Hz]) to 12000 [cents] (16717 [Hz]) was obtainedby performing the wavelet transform with a Gabor wavelet,
a frequency interval of 10 [cents], and a shifting intervalof 10 [ms], i.e., M = 1200 and N = 3000. The otherquantities were I = 88, J = 3, R = 20, P = 13, and= 25. The priors were set to be less informative, i.e.,Proceedings of the 16th ISMIR Conference, Málaga, Spain, October 26-30, 2015
Filter learning
Unsupervised

12000 [cents]
9600
7200
4800
2400
0
12000
9600
7200
4800
2400
0
108 [MIDI note number]

HMM
3

Observation

Supervised(open test)
Supervised(close test)

Reconstruction

Pitch activation

21
108

140

[dB]

100

Final output

F-measure: 69.4%

108

3
3
3
3
3
3

3
3
3

R
55.3
62.2
62.4
67.4
62.4
69.9
59.4
67.4

P
57.9
60.2
64.3
64.2
67.0
64.5
69.1
67.8

F
56.6
61.2
63.4
65.8
64.4
67.3
63.9
67.6

Table 1. Experimental results of multipitch analysis for 30
piano pieces labeled as ENSTDkCl.
120

21

HPSS

91

All-pole filter learned from“ENSKDkCl” isolated piano sounds

80
60

All-pole filter learned from
“MUS-mz̲333̲3̲ENSTDkCl” data
(superimposed piano sounds)

All-pole filter learned from
“StpkBGCl” isolated piano sounds

40
20
0

21

5

10

15

20

Ground truth25
30 [s]

Figure 5. Analysis of MUS-mz 333 3 ENSTDkCl.
↵✓ = ↵ = aH = 1, and bH = Eemp [Xmn ] 1 . Since X
contained only piano sounds, the truncation level J = 3
worked well (two filters were degenerated in this experiment, i.e., J + = 1). The values of {µni }Ii=1 were initialized as the frequencies corresponding to the 88 keys ranging from 900 [cents] to 9600 [cents]. The value of each↵jp (1  p  P ) was drawn from a Gaussian with a zeromean and a small variance of 0.01. The variational posteriors were initialized as the corresponding priors.
The proposed model was tested under possible combinations of preprocessing (with or without HPSS) and postprocessing (thresholding or Viterbi decoding). HPSS was
performed in the log-frequency domain. The model with asingle filter (J = J + = 1) was also tested in a supervised
setting. A set of filter coefficients a1 was pretrained from
264 isolated sounds of the same or different piano (ENSTDkCl in a closed test or SptkBGCl in an open test) by using
LPC, and kept constant during multipitch analysis.
The estimation results were evaluated in terms of theframe-level recall/precision rates and F-measure as in [24]:
P
P
cn
cn
2RP
, (36)
R = Pn , P = P n , F =
r
e
R
+P
n
n
n
n
where rn , en , and cn are the numbers of ground truth, estimated and correct pitches on frame n, respectively. The
threshold ⌘ was determined as ⌘ = 10 1.3 without HPSSand ⌘ = 10 1.5 with HPSS.
4.2 Experimental Results
The experimental results shown in Figure 5 and Table 1 indicate the great potential of the iSDAP model. The model
supervised in the open condition (67.3%) significantly outperformed the iCAR model formulated in the linear frequency domain (48.4%) [26] and tied with the state-ofthe-art methods, e.g., harmonic NMF (67.7%) [24], NMF
with group sparsity (71.3%) [21], and NMF with Hellinger

-20

0

2400

4800

7200

9600 [cents] 12000

Figure 6. All-pole filters learned from isolated sounds or a
piano piece (mixed sounds) in the log-frequency domain.
sparse coding (66.5%) [22]). While many recent methods
need to pretrain a dictionary of basis spectra for reasonable
decomposition [2,4,21,24], our model works well (65.8%)
even in the completely unsupervised condition. As shownin Figure 6, a filter learned from a music signal dropped
faster than the pretrained filters because the model failed tocapture higher-order overtones even in the log-frequency
domain due to the strong inharmonicity of piano sounds.
Nonetheless, the learned filter acted as an effective constraint on the relative weights of harmonic partials.
There would be much room for improving the performance. KL-NMF [18] and IS-NMF [10] are special casesof -divergence NMF [11, 20] with = 1, 0, respectively.
It was reported that the use of an intermediate divergencewith = 0.5 significantly improves the performance by
about 5% [24]. Similar findings were reported in the context of source separation [13]. This calls for the use of the
Tweedie likelihood instead of the Poisson likelihood [6].
5. CONCLUSION
We presented a new nonparametric Bayesian approach tosource-filter NMF called infinite superimposed discrete allpole (iSDAP) modeling that can decompose a wavelet spectrogram into three kinds of factors, i.e., harmonic sources,
all-pole filters, and time-varying gains of source-filter pairs.
Our model clearly outperformed its counterpart called theiCAR model formulated in the linear frequency domain.
One important research direction is to build a unified modelof harmonic and percussive sounds. To bridge the gap between multipitch analysis and music transcription, we plan
to incorporate a prior distribution on the time-frequency
positions of musical notes into a Bayesian framework.
Acknowledgment: This study was partially supported by JST
OngaCREST Project, JSPS KAKENHI 24220006, 26700020,
and 26280089, and Kayamori Foundation.

92

Proceedings of the 16th ISMIR Conference, Málaga, Spain, October 26-30, 2015
6. REFERENCES[1] R. Badeau and A. Ozerov. Multiplicative updates formodeling mixtures of non-stationary signals in thetime-frequency domain. In European Signal Processing Conference (EUSIPCO), 2013.
[2] E. Benetos, R. Badeau, T. Weyde, and G. Richard.
Template adaptation for improving automatic musictranscription. In International Society for Music Information Retrieval Conf. (ISMIR), pages 175–180, 2014.
[3] N. J. Bryan, G. Mysore, and G. Wang. Source separation of polyphonic music with interactive userfeedback on a piano roll display. In International Society for Music Information Retrieval Conference (ISMIR), pages 119–124, 2013.
[4] J. J. Carabias-Orti, T. Virtanen, P. Vera-Candeas,
N. Ruiz-Reyes, and F. J. Cañadas-Quesada. Musical instrument sound multi-excitation model for nonnegative spectrogram factorization. IEEE Journal of
Selected Topics in Signal Proc., 5(6):1144–1158, 2011.
[5] A. T. Cemgil. Bayesian inference for nonnegative matrix factorisation models. Computational Intelligence
and Neuroscience, 2009:Article ID 785152, 2009.
[6] U. Şimşekli, A. Cemgil, and Y. K. Yılmaz. Learningthe -divergence in Tweedie compound Poisson matrix
factorization models. In International Conference on
Machine Learning (ICML), pages 1409–1417, 2013.
[7] J.-L. Durrieu, G. Richard, B. David, and C. Févotte.
Source/filter model for unsupervised main melody extraction from polyphonic audio signals. IEEE Transactions on Audio, Speech, and Language Processing,18(3):564–575, 2010.
[8] A. El-Jaroudi and J. Makhoul. Discrete all-polemodeling. IEEE Transactions on Signal Processing,39(2):411–423, 1991.
[9] V. Emiya, R. Badeau, and B. David. Multipitch estimation of piano sounds using a new probabilistic spectral smoothness principle. IEEE Transactions on Audio, Speech and Language Processing, 18(6):1643–
1654, 2010.

[10] C. Févotte, N. Bertin, and J.-L. Durrieu. Nonnegative
matrix factorization with the Itakura-Saito divergence:With application to music analysis. Neural Computation, 21(3):793–830, 2009.
[11] C. Févotte and J. Idier. Algorithms for nonnegative
matrix factorization with the beta-divergence. NeuralComputation, 23(9):2421–2456, 2011.
[12] D. FitzGerald. Harmonic/percussive separation using
median filtering. In International Conference on Digital Audio Effects (DAFx), 2010.
[13] D. FitzGerald, M. Cranitch, and E. Coyle. On the useof the beta divergence for musical source separation. In
Irish Signals and Systems Conf., pages 1–6, 2008.
[14] R. Hennequin, R. Badeau, and B. David. NMF with
time-frequency activations to model nonstationary audio events. IEEE Transactions on Audio, Speech, andLanguage Processing, 19(4):744–753, 2011.

[15] M. Hoffman, D. Blei, and P. Cook. Bayesian nonparametric matrix factorization for recorded music. In International Conference on Machine Learning (ICML),
pages 439–446, 2010.
[16] F. Itakura and S. Saito. Analysis synthesis telephonybased on the maximum likelihood method. In International Congress on Acoustics (ICA), pages C17–C20,
1968.
[17] H. Kameoka and K. Kashino. Composite autoregressive system for sparse source-filter representation of
speech. In IEEE International Symposium on Circuits
and Systems (ISCAS), pages 2477–2480, 2009.
[18] D. Lee and H. Seung. Algorithms for non-negative
matrix factorization. In Neural Information Processing
Systems (NIPS), pages 556–562, 2000.
[19] T. Nakamura, K. Shikata, N. Takamune, and
H. Kameoka. Harmonic-temporal factor decomposition incorporating music prior information for informed monaural source separation. In International
Society for Music Information Retrieval Conference(ISMIR), pages 623–628, 2014.
[20] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano,
N. Ono, and S. Sagayama. Convergence-guaranteed
multiplicative algorithms for non-negative matrix factorization with beta divergence. In International Workshop on Machine Learning for Signal Processing(MLSP), pages 283–288, 2010.
[21] K. O’Hanlon and M. D. Plumbley. Polyphonic pianotranscription using non-negative matrix factorisation
with group sparsity. In International Conference on
Acoustics, Speech, and Signal Processing (ICASSP),
pages 2214–2218, 2014.
[22] K. O’Hanlon, M. Sandler, and M. D. Plumbley. Matrix factorisation incorporating greedy Hellinger sparse
coding applied to polyphonic music transcription. In
International Conference on Acoustics, Speech, andSignal Processing (ICASSP), pages 3112–3116, 2015.
[23] P. Smaragdis, C. Févotte, G. Mysore, N. Mohammadiha, and M. Hoffman. Dynamic source separation
using nonnegative factorizations: A unified view. IEEE
Signal Processing Magazine, 31(3):66–75, 2014.
[24] E. Vincent, N. Bertin, and R. Badeau. Adaptive harmonic spectral decomposition for multiple pitch estimation. IEEE Transactions on Audio, Speech and Language Processing, 18(3):528–537, 2010.
[25] T. Virtanen and A. Klapuri. Analysis of polyphonic audio using source-filter model and non-negative matrix
factorization. In NIPS Workshop on Advances in Models for Acoustic Processing, 2009.
[26] K. Yoshii and M. Goto. Infinite composite autoregressive models for music signal analysis. In International
Society for Music Information Retrieval Conference(ISMIR), pages 79–84, 2012.