2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

October 16-19, 2011, New Paltz, NY

BAYESIAN NONPARAMETRIC SPECTROGRAM MODELINGBASED ON INFINITE FACTORIAL INFINITE HIDDEN MARKOV MODEL
Masahiro Nakano† , Jonathan Le Roux‡ , Hirokazu Kameoka‡ , Tomohiko Nakamura† , Nobutaka Ono† and Shigeki Sagayama††

Graduate School of Information Science and Technology, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan‡NTT Communication Science Laboratories, NTT Corporation,
3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198, Japan

ABSTRACTThis paper presents a Bayesian nonparametric latent source discovery method for music signal analysis. In audio signal analysis, an
important goal is to decompose music signals into individual notes,
with applications such as music transcription, source separation ornote-level manipulation. Recently, the use of latent variable decompositions, especially nonnegative matrix factorization (NMF), hasbeen a very active area of research. These methods are facing two,mutually dependent, problems: ﬁrst, instrument sounds often exhibit time-varying spectra, and grasping this time-varying nature is
an important factor to characterize the diversity of each instrument;
moreover, in many cases we do not know in advance the number ofsources and which instruments are played. Conventional decompositions generally fail to cope with these issues as they suffer from the
difﬁculties of automatically determining the number of sources andautomatically grouping spectra into single events. We address boththese problems by developing a Bayesian nonparametric fusion of
NMF and hidden Markov model (HMM). Our model decomposesmusic spectrograms in an automatically estimated number of components, each of which consisting in an HMM whose number ofstates is also automatically estimated from the data.
Index Terms— Nonnegative matrix factorization (NMF), Hierarchical Dirichlet process (HDP), Inﬁnite hidden Markov model(iHMM), Gamma process, Collapsed variational Bayes (CVB)
1. INTRODUCTION
In acoustic signal processing, the use of latent variable decompositions, especially NMF [1], has been a very active area of research.
In the context of audio signal analysis based on NMF, the observed
signal (magnitude or power spectrogram) is approximated by a linear combination of spectral bases which can be regarded as the representatives of the various spectral patterns in the signal.
This paper focuses on the two problems which standard NMF
faces in music signal analysis. Firstly, one needs to deal with
the non-stationarity of instrument sounds [2, 3]. Hopefully, each
spectral basis should be the representative of a single instrument
sound. In standard NMF, this implies that the elementary components (which one would like to correspond to one pitch of one instrument) of the analyzed sound are assumed to be nearly stationary.
However, real world sounds often exhibit non-stationary spectral
characteristics. For example, a piano note would be more accurately
characterized by a succession of several spectral patterns such as attack, decay, sustain and release. As another example, singing voices
and stringed instruments feature a particular musical effect, vibrato.

978-1-4577-0693-6/11/$26.00 c 2011 IEEE

Spectra of instrument sounds with such variations are typically difﬁcult to handle by standard NMF. Indeed, a single instrument soundwill tend to be modeled as the sum of several spectral bases, leading
to the difﬁcult problem of determining which sources each component belongs to. Standard NMF thus requires some post-processing
to group the bases into single notes. Secondly, it is difﬁcult to estimate the number of notes which are included in the observed polyphonic music signals [4, 5]. Most methods require that the numberof sources be speciﬁed in advance, or found with model selection
techniques [5]. This problem is related to the non-stationarity ofthe spectrum of instrumental sounds, because the variations in the
spectral patterns of a sound make it even more difﬁcult to estimate
the number of sources. These problems should thus be addressed
simultaneously.
We do so by developing a Bayesian nonparametric fusion of
NMF and hidden Markov model (HMM), which can be regarded as
an extension of the factorial hidden Markov model (FHMM). The
number of components and the number of states have a strong impact on the standard FHMM. To overcome this issue, we place herea prior over all their possible combinations via a Bayesian nonparametric framework. Our model decomposes music spectrograms in
an automatically estimated number of components, each of which
consisting in an HMM whose number of states is also automaticallyestimated from the data.
Conventional Bayesian nonparametric NMF, here reformulatedwith the generalized KL-divergence, is explained in section 2. Section 3 presents our model based on Bayesian nonparametric combined NMF and HMM. An efﬁcient inference algorithm based on
variational Bayes (and collapsed VB) is derived in Section 4. Section 5 discusses potential improvements to the model based on adifferent construction of the prior. Experimental results are shownin Section 6.
2. GAMMA PROCESSNONNEGATIVE MATRIX FACTORIZATIONAn NMF approach to audio signal analysis is typically based on
the assumption that a magnitude or power spectrogram Y =(Yω,t )Ω×T ∈ R≥0,Ω×T , where ω = 1, · · ·, Ω is a frequency bin
index, and t = 1, · · ·, T is a time frame index, can be modeled
as the product of two non-negative matrices, H = (Hω,d )Ω×D ∈R≥0,Ω×D and UP= (Ud,t )D×T ∈ R≥0,D×T . This can be written as Yω,t ≈d Hω,d Ud,t , where D is the number of bases
hd = [H1,d , · · ·, HΩ,d ]T . The term “component” is often used torefer to a spectral portion corresponding to hd and Ud,t . Although
the discrepancy measure between Y and HU can be deﬁned in3252011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

many ways, we will focus in this paper on the use of the generalized
Kullback-Leibler (KL) divergence in the magnitude domain, since
this choice has been shown to lead to good music source separation
performance [6].
The standard NMF assumes that the number of components D
is known in advance. In practice, we rarely know the adequate value
for this number and must resort to expensive techniques such asmodel selection to estimate it. Recently, Bayesian nonparametricapproaches have proven to be effective for the problem of automatically inferring model complexities [5, 7]. The Gamma process nonnegative matrix factorization (GaP-NMF) [5] introduces ahidden vector of nonnegative values „, where each element θd isthe overall gain of the corresponding source d. While the original GaP-NMF [5] is basically a nonparametric extension of the
Itakura-Saito divergence-based NMF, we start here by formulating a nonparametric extension of the generalized KL divergencebased NMF using the Gamma process prior. The observed spectrogram Y and the overall gains „ are expressed according to thefollowing generative process: Cω,t,d ∼ Poisson(θP d Hω,d Ud,t ),θd ∼ Gamma(η/D, ηλ), Yω,t ∼ δ(Yω,t − d Cω,t,d ) where
C 1 , . . . , C D can be interpreted as latent components. Note that D
works as the truncation level. As D increases towards inﬁnity, theoverall gains „ can be considered to be drawn from a Gamma process with shape parameter η and inverse-scale parameter ηλ [8, 5].
When the truncated level D is sufﬁciently large, we can expect thatan adequate number of components will be active and the rest will
be suppressed. We thus don’t need to determine the adequate number of components in advance. Instead, we only have to prepare a
sufﬁciently large number D of them.
3. BAYESIAN NONPARAMETRICCOMBINED NMF AND HMM
FOR MODELING MUSIC SPECTROGRAMS
GaP-NMF gives us an adequate number of components required by
the model to describe the data, but we would like each of these components to correspond to a single event, such as “C of the piano”
or “A of the violin”. This would likely be the case if we could
assume that the spectrum of a note of a musical instrument can be
represented through a single spectral basis whose amplitude is modulated in time, but its variations in time are actually much richer.
Recently, an HMM-based approach has been proposed toovercome this problem [2]. The spectrogram of the observed
signal is modeled under the assumption that it is composedof spectral patterns which are themselves composed of a limited number of Markov-chained states. To introduce this concept into GaP-NMF, we consider state-transition bases H =
(k)
(k)
(k)
{(Hω,1 )Ω×K , . . . , (Hω,D )Ω×K }, where (Hω,d )Ω denotes the kth possible state for the spectral basis of the d-th component andK denotes the number of states. If we let Z = (Zd,t )D×T ∈ N
denote which spectral basis state of the d-th component is acti(k)vated at time t, the generative model can be expressed as Hω,d ∼Gamma(aH , bH ), Ud,t | Wd,t ∼ Gamma(aU , aU Wd,t ), Wd,t |Ud,t−1 ∼ Gamma(aU , aU Ud,t−1 ), θd ∼ Gamma(η/D, ηλ) and“”(Zd,t )(k)Cω,t,d | (Hω,d )K
Ud,t ,
k=1 , Ud,t , Zd,t ∼ Poisson θd Hω,d“”X
Yω,t ∼ δ Yω,t −Cω,t,d ,(1)dwhere W = (Wd,t )D×T are auxiliary variables. The Gamma chainof U and W promotes temporal continuity, often encountered inreal-world sounds [13].

October 16-19, 2011, New Paltz, NY

The number of states is likely to be an important factor to characterize the diversity of each instrument. Hopefully, the adequate
number of states should be assigned in response to the tested instrumental sounds. To address this problem, we introduce a Bayesiannonparametric approach to the state-transition bases, similarly tothe inﬁnite HMM [9]. Our model is constructed based on theBayesian nonparametric HMM using hierarchical Dirichlet process (HDP) [9, 10, 11]. A two-level HDP can be used to develop
an HMM with an inﬁnite state space. It is a collection of DPs,
{G1 , G2 , . . . }, sharing a base distribution G0 : Gj ∼ DP(α, G0 ),for each j. G0 is also drawn from a DP with a base distribution F :G0 ∼ DP(γ, F ). α, γ are referred to as concentration parameters,
which can be interpreted as controlling the impact of the base distribution. The intuitive interpretation of the HDP in our model canbrieﬂy be described as follows: the base distribution F is a distribution over the nonnegative spectrum space and draws independently
the atoms of the top-level DP G0 . As a draw from DP, G0 is almost
surely discrete, leading to the discreteness of the HMM. The factthat all Gj share G0 enables them to share the same set of atoms
with only different weights. It ensures that transitions are allowed
only on the same set of atoms.
Various constructions of the HDP-HMM prior have been proposed [9, 12, 11, 15]. Although, as discussed in Section 5, somealternatives may lead to potentially better models, we choose topresent in this paper an intuitive construction based on the definition of the HDP and the stick-breaking process [9, 12]. Let
˛ d = βd,1 , βd,2 , . . . , βd,k , . . . be generated from a stick-breaking
process: ˛ d ∼ GEM(γ), that is, βd,m|γ ∼ Beta(1, γ), βd,m =Qm−1βd,m l=1 (1 − βd,l ). Let ı d,k ∼ DP(α, ˛ d ) where ı d,k denotes the state-speciﬁc transition distribution, that is, πd,k,k shows
the state-transition probability from state k to state k . Then, thestates are encouraged to be similar because E[πd,k,k |˛ d ] = βd,k .
The whole model can be expressed by combining Eq. (1) and thefollowing generative process:˛ d ∼ GEM(γ) , ı d,k ∼ DP(α, ˛ d )Zd,t | Zd,t−1 , (πd,k )∞(2)
k=1 ∼ ı d,Zd,t−1
The key property is that only meaningful elements of the unlimited
number of states are active in essence. We thus do not need to ﬁnd
the suitable number of states in advance, as it is assigned in responseto the tested instrumental sounds.
To summarize, we introduced two inﬁnite models into the standard FHMM. We shall refer to this model as the inﬁnite factorial
inﬁnite hidden Markov model (iFiHMM).
4. VARIATIONAL INFERENCE
Inference algorithms for models with an HDP prior are mostly
based on sampling methods. For HDP-HMM, various types ofMarkov chain Monte Carlo methods have been proposed [9, 10, 11].
In the context of Bayesian inference, there are alternatives based on
variational Bayes (VB). Especially, for large-scale problems, VB
may prove to be one of the most convenient methods. In this paper,
we use variational inference for our model because music spectrograms often comprise a large number of parameters.
The VB approach in general assumes a factorized form of theposterior distribution. This implies that the parameters are assumed
to be independent of each other, and such an assumption is likely tooften degrade the performance of the inference when it is not met. It
is the case in iFiHMM, where ı d is likely to have a strong impact on
Z d . In order to better approximate the posterior, we thus use here326

2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

the recently developed collapsed VB (CVB) for HDP, integratingout certain parameters while assuming that other latent variables
are independent [12].
To derive collapsed variational updates, we ﬁrst integrate out ı,leading to a joint distribution over Z, α, ˛ as follows:YnY Γ(αβd,k + nd,j,k ) oΓ(α)
p(Z | α, ˛) =(3)Γ(α + nd,j,· )Γ(αβd,k )d,jkwhere nd,j,k counts the number of times that thePstate transition
from state j to state k is drawn and nd,j,· =k nd,j,k . The
gamma function in Eq. (3) makes the hyperparameter posterior updates cumbersome. We thus introduce auxiliary variables sd,j,k taking integral values, similarly to [12]. The joint distribution over the
expanded system is obtained as:oYn
Y » nd,j,k –Γ(α)
p(Z, s | α, ˛) =(αβd,k )sd,j,kΓ(α + nd,j,· )sd,j,kd,jkˆ ˜where ·· denotes the unsigned Stirling numbers of the ﬁrst kind.
Given the observed spectrogram Y , we want to compute theposterior distribution p(C, „, H, U , Z, ˛, s|α, η, λ, a, b). In thefollowing, we will refer to the parameters of interest as Ξ. CVB approximates the posterior with the variational distribution q(Ξ) =q(C)q(„)q(H)q(U )q(Z)q(˛)q(s|Z). We optimize the factorized distribution of each parameter iteratively while holding the others ﬁxed. Space constraints allow only for a sketch of the update
rules for q(Z):“ X(k)q(Zd,t = k) ∝ exp −E[θd Hω,d Ud,t ]
+

X

ω

(k)

E[Cω,t,d log Hω,d Ud,t ] + E[log(αβd,k + n−td,zd,t−1 ,k )]ω+E[log(αβd,zd,t+1 + n−td,k,zd,t+1 + I(zd,t−1 = k)I(zd,t+1 = k))]”−E[log(α + n−td,k,· + I(zd,t−1 = k))]where n−td,j,k counts the number of times that state transition fromstate j to state k is drawn, excluding the number of transitions
Zd,t−1 → Zd,t or Zd,t → Zd,t+1 .
The update rules for q(C), q(H) and q(U ) are derived bymaking a slight modiﬁcation (adding „) to [3]. The updates of q(˛)
and q(s|Z) are derived similarly to [12]. The updates of q(θ) derived in [5] are generalized to account for the state transitions between spectral bases. It is computationally difﬁcult to obtain someof the expectations. However, we can use the same tractable approximations employed in [12].
5. ALTERNATE CONSTRUCTION
FOR THE HDP-HMM PRIOR
It is worth noting that a different VB approach for HDP has recently
been introduced in [11, 15], based on the Chinese restaurant franchise (CRF) [9]. The whole model can be expressed by combining
Eq. (1) and:˛ d ∼ GEM(γ) , ı d,j ∼ GEM(α) , ςd,j,i ∼ ˛ dηd,j,t ∼ ı d,j , Zd,t = ςd,Zd,t−1 ,ηd,Zd,t−1 ,t(4)

Because this construction does not involve variables having a strong
inﬂuence on others, there is no special incentive to using collapsed
VB. Moreover, owing to the full conjugacy of the construction, VB
can be straightforwardly applied and the update rules obtained inclosed form. This construction thus seems very promising, and weshall report on it in more details in a future publication.

October 16-19, 2011, New Paltz, NY

6. EXPERIMENTS
We present some results on the application of our algorithm to audio signals, for fully unsupervised sound separation. All data were
downmixed to mono and downsampled to 16 kHz. The magnitude
spectrogram was computed using the short time Fourier transformwith 32 ms long Hanning window and with 16 ms overlap.
We ﬁrst generated synthetic data (shown in Fig. 1 (a)) consisting in a mixture of piano (C), violin (E) and ﬂute (G), chosen tohave overlapping harmonic components: ﬁrst, each note is playedalone in turn, then all the combinations of two notes are played andﬁnally all notes are played simultaneously. To clarify the effect ofthe present method, we used as a comparison standard NMF where
the number D of bases is ﬁxed to 3. The result is shown in Fig. 1.
Standard NMF, even when given the proper D, is unable to represent correctly the time-varying spectra, such as the attack part ofthe piano and the vibrato of the violin. It would need more spectral
bases to capture the time-varying spectra, but some post-processing
would then be required to group the bases into single notes. By contrast, the proposed method (truncation levels: D = 10 and K = 30)
automatically factorizes the observed spectrogram into the mainly
active components (shown in Fig. 2). The relative level of the remaining components is −109.2 dB. The proposed method is thus
able to automatically ﬁnd the adequate number of components andgroup the proper number of Markov-chained spectral bases to represent non-stationary spectra.
As an example of real application, we also present results on thetask of note-level manipulation of real-world music signals. Some
audio samples as well as a description of the procedure we used areavailable at http://hil.t.u-tokyo.ac.jp/˜mnakano/.
7. CONCLUSION
This paper presented a Bayesian nonparametric combined NMF and
HMM framework for modeling music spectrograms. We showed
through experiments that our model is capable of automatically determining the number of sources and of modeling the variety of thetime-varying spectra of each instrument sound, via a Bayesian nonparametric approach. In the future, we will apply this model to various tasks related to music signal analysis. We also plan to extend
our model to a semi-supervised setup.
8. REFERENCES[1] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, pp. 788–791, Oct. 1999.
[2] A. Ozerov, C. Févotte, and M. Charbit, “Factorial scaled hidden
Markov model for polyphonic audio representation and source separation,” in Proc. WASPAA, 2009.
[3] M. Nakano, J. Le Roux, H. Kameoka, N. Ono and S. Sagayama,“Inﬁnite-state spectrum model for music signal analysis,” in Proc.
ICASSP, pp. 1972–1975, 2011.
[4] V. Y. F. Tan and C. Févotte, “Automatic relevance determination in
nonnegative matrix factorization,” in Proc. SPARS, 2009.
[5] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametric matrix
factorization for recorded music,” in Proc. ICML, pp. 641–648, 2010.
[6] D. FitzGerald, M. Cranitch and E. Coyle, “On the use of the Beta
divergence for musical source separation,” in Proc. ISSC, 2009.
[7] M. N. Schmidt and M. Mørup, “Inﬁnite non-negative matrix factorization,” in Proc. EUSIPCO, 2010.
[8] J. F. C. Kingman, “Poisson processes,” Oxford University Press, 1993.
[9] Y. Teh, M. Jordan, M. Beal and D. Blei, “Hierarchical Dirichlet processes,” in Proc. Journal of the American Statistical Association, 101,
pp. 1566–1581, 2004.

327

2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

October 16-19, 2011, New Paltz, NY(TGSWGPE[=M*\?

(TGSWGPE[=M*\?























6KOG=U?



6KOG=U?

5VCVG

(a) Input spectrogram, a mixture of piano, violin and ﬂute










6KOG=U?

(a) Estimated model: Piano (C)













(TGSWGPE[=M*\?

(TGSWGPE[=M*\?



6KOG=U?

(b) Estimated model: Piano (C)













6KOG=U?


5VCVG(TGSWGPE[=M*\?














6KOG=U?

(c) Estimated model: Violin (E)

6KOG=U?

(b) Estimated model: Violin (E)(TGSWGPE[=M*\?

(TGSWGPE[=M*\?



















6KOG=U?









6KOG=U?

(d) Estimated model: Flute (G)
5VCVG

Figure 1: Original spectrogram (a), estimated model Hω,d Ud,t (b), (c) and(d) obtained by the standard NMF (D=3). The standard NMF cannot explicitly represent the time-varying nature of the spectra of instrument sounds.
We can see that it does not capture the attack part of piano and ﬂute, and thevibrato of violin.
[10] J. V. Gael, Y. Saatci, Y. W. Teh and Z. Ghahramani, “Beam sampling
for the inﬁnite hidden Markov model,” in Proc. ICML, 2008.
[11] E. B. Fox, E. B. Sudderth, M. I. Jordan and A. S. Willsky, “An HDPHMM for systems with state persistence,” in Proc. ICML, 2008.
[12] Y. W. Teh, K. Kurihara and M. Welling, “Collapsed variational inference for HDP,” in Proc. NIPS, 2008.
[13] A. T. Cemgil and O. Dikmen, “Conjugate Gamma Markov randomﬁelds for modelling nonstationary sources,” in Proc. ICA, pp. 697–705,
2007.
[14] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC music
database: Popular, classical, and jazz music database,” in Proc. ISMIR,
pp. 287–288, 2002.












6KOG=U?

(c) Estimated model: Flute (G)Figure 2: Estimated model obtained by applying the proposed method to(Z)

the input signal (Fig. 1 (a)): E[θd Hω,dd,t Ud,t ] (top) and q(Zd,t = k)(bottom) of each component. We set the truncation levels to D = 10 and
K = 30, and hyperparameters to α = γ = 1, aH = bH = 0.001,
P
aU = 1, η = 0.1 and λ = ΩT / ω,t Yω,t . Three components are mainly
active and the relative level of the remaining components is −109.2 dB. The
adequate numbers of components and states are automatically optimized depending on the observed signal.
[15] C. Wang, J. Paisley and D. M. Blei, “Online variational inference forthe hierarchical Dirichlet process,” in Proc. AISTATS, 2011.

328