INFINITE-STATE SPECTRUM MODEL FOR MUSIC SIGNAL ANALYSIS
Masahiro Nakano† , Jonathan Le Roux‡ , Hirokazu Kameoka‡ , Nobutaka Ono† , 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 nonparametric Bayesian extension of nonnegative matrix factorization (NMF) for music signal analysis.
Instrument sounds often exhibit non-stationary spectral characteristics. We introduce infinite-state spectral bases into NMF to represent
time-varying spectra in polyphonic music signals. We describe our
extension of NMF with infinite-state spectral bases generated by
the Dirichlet process in a statistical framework, derive an efficient
optimization algorithm based on collapsed variational inference, andvalidate the framework on audio data.
Index Terms— Nonnegative matrix factorization (NMF),
Dirichlet process, Collapsed variational Bayes, NonparametricBayes
1. INTRODUCTION
Nonnegative matrix factorization (NMF) [1] is an unsupervised decomposition technique allowing the representation of two-dimensional nonnegative data as a linear combination of a few meaningful
elementary bases. In particular, NMF has been applied successfully
to music spectrograms in audio signal processing, with such applications as automatic music transcription or sound source separation.
NMF is able to project all signals that have a similar spectral shapeon a single basis, allowing one to represent a variety of phenomena
efficiently using a very compact set of spectrum bases.
While real world sounds typically exhibit non-stationary spectral characteristics, the standard NMF implicitly assumes that eachelementary component (which is expected to correspond to a single
note activation) of a signal under analysis is represented by a ”rank
1” spectrogram. This means that the spectrum of each note is assumed to be constant over time up to a scale factor. Thus, learning
an important spectral variability with standard NMF would require
to use a large number of bases, and some post-processing to groupthe bases into single events.
Several approaches have been proposed to overcome this problem. For example, nonnegative matrix factor deconvolution [2]introduces temporal-spectral bases, while [3] considers an extension of NMF where temporal activations are replaced with timefrequency activations based on a source/filter model. Another kindof approaches based on state variations of the spectral patterns has
been an active area of research to model sounds whose spectral
characteristics evolve other time, such as the piano with attack andrelease parts. Factorial scaled hidden Markov model [4] or NMF
with Markov-chained bases [5] intend to represent time-varying
spectra as state transitions through a limited and fixed number ofspectral patterns. However, the important problem of the automatic

978-1-4577-0539-7/11/$26.00 ©2011 IEEE1972

determination of the number of states remained to be solved. For
example, while a piano note is likely to be well expressed with a few
spectral patterns, one would expect a singing voice to require many.
The number of states is thus an important factor to characterize
the diversity of each instrument. Hopefully, the adequate numberof states should be assigned in response to the tested instrumental
sounds.
This paper proposes to extend NMF to represent time-varying
spectral patterns of instrument sounds by introducing infinite-state
spectral bases. The number of states of basis is adaptively optimized depending on the instrument sounds. Bayesian NMF with
deformable bases is explained in Section 2. Section 3 presents theinfinite-state spectrum model, an extension of NMF with deformablebases. Section 4 derives an efficient optimization algorithm basedon collapsed variational inference. Validation experiments are presented in Section 5.
2. BAYESIAN NONNEGATIVE MATRIX FACTORIZATIONWITH DEFORMABLE BASESNMF applied to audio signal analysis is based on a signal model
where the magnitude or power spectrogram Y = (Yω,t )Ω×T ∈R≥0,Ω×T , where ω = 1, . . . , Ω is a frequency bin index, andt = 1, . . . , T is a time frame index, is factorized into nonnegative
parameters, H = (Hω,d )Ω×D ∈ R≥0,Ω×D and U = (Ud,t )D×T ∈R≥0,D×T . This can be written as
XHω,d Ud,t ,(1)
Yω,t ≈dwhere D is the number of bases hd = [H1,d , . . . , HΩ,d ]. The termcomponent is used to refer to one basis hd and its time-varying gain
Ud,t . The bases can be considered as spectral patterns which arefrequently observed.
Hopefully, one factorized component should represent a single event, but audio events actually often have varying spectral patterns. We thus propose to consider deformable bases
(k)
(k)
(k)
H = {(Hω,1 )Ω×K , . . . , (Hω,D )Ω×K }, where (Hω,d )Ω denotes
the k-th possible state for the spectral basis of the d-th component.
If we let Z = (Zd,t )D×T ∈ N denote which spectral basis stateof the d-th component is activated at time t, NMF with deformablebases can be written as
X (Zd,t )Yω,t ≈
Hω,d Ud,t .
(2)dAssuming that the generalized Kullback-Leibler divergence is usedas the divergence measure, the model can be expressed as the following generative model, similarly to [6]:X(Z)
Yω,t =
Cω,t,d , Cω,t,d ∼ Poisson(Cω,t,d | Hω,dd,t Ud,t ) .
d

ICASSP 2011

It is possible to introduce various prior structures for H and U
according to the data and the requirements of the considered applications. For the spectral bases H, we here propose to use a Gamma
distribution, which is the conjugate prior to the Poisson distribution:(k)

Hω,d ∼ Gamma(aω , bω ) ,(3)

the primary motivation for this choice being computational convenience. To promote temporal continuity, often encountered in realworld sounds, we introduce the following prior distribution on thetime-varying gains U , similarly to [7]:Wd,t | Ud,t
Ud,t+1 | Wd,t∼∼Gamma(α, αUd,t ) ,Gamma(β, βWd,t ) ,

(4)
(5)

where W = (Wd,t )D×T are auxiliary variables. If we set α and β
to large values, the prior will enforce the temporal smoothness of U .
3. INFINITE-STATE SPECTRUM MODEL
It is important to determine the number of states for each component.
For example, a piano note would be often characterized by several
spectral patterns such as “attack”, “decay”, “sustain” and “release”.
Three or four bases may thus be needed to represent the spectrogramof the piano. As another example, singing voices and stringed instruments feature a particular musical effect, vibrato. Vibrato may need
to be expressed using many spectral patterns. Thus, it is desirable
to automatically determine the appropriate number of deformablebases for fitting to each instrument sound. To achieve this, we introduce the Dirichlet process (DP) into the deformable bases. DP
is popularly used as a nonparametric prior in hierarchical Bayesian
specification [8, 9]. Several practical methods for the constructionof DP have been derived [8, 10].
Let us first show how to construct it here based on a symmetric
Dirichlet prior. We consider the sequence {Zd,1 , . . . , Zd,T } of state
indices for the d-th component as a sequence of discrete indicatorvariable, where each Zd,t can take on values 1, . . . , K with proportions given by ıd = {πd,1 , . . . , πd,K }. The joint distribution of thesequence is multinomial
K
(k)
Y
nd
p(Zd,1 , . . . , Zd,T | ıd ) =
πd,k
,
(6)
k=1
P
(k)
where nd = Tt =1 δ(Zd,t − k) and δ(·) denotes the Kronecker(k)delta function to count the number of times nd that Zd,t = k (t =1, . . . , T ) has been drawn. Let us give the mixing proportions asymmetric Dirichlet prior, which is a conjugate prior of the multinomial distribution, with positive concentration hyperparameter ‚ ={γ1 , . . . , γD }:p(πd | γd ) ∼ Dirichlet(γd /K, . . . , γd /K) .
(7)

The critical property of the DP is shown in the conditional probability of Zd,t given the setting of all other indices in the sequence{Zd,1 , . . . , Zd,t−1 , Zd,t+1 , . . . , Zd,T } (denoted Zd,−t ) under K →∞:p(Zd,t = k | Zd,−t , γd )8(k)>< nd,−t(k ∈ {1, . . . , K} i.e. represented)= T − 1 + γdγd>:(for all unrepresented k, combined)T − 1 + γd(k)where nd,−t counts the number of times that Zd,t = k has beendrawn in Z d,−t . As we can see, Zd,t tends to choose an alreadypopular state. The concentration parameter γd controls the tendency

1973

to populate a previously unrepresented state. Each component thus
tends to keep an adequate number of states depending on the observed signals.
Another formulation of the DP has also been proposed in termsof a stick-breaking construction [10]:k−1YVd,k ∼ Beta(1, γd ) , πd,k (Vd ) = Vd,k
(1 − Vd,j ) ,(8)
j=1

where Vd = {Vd,1 , Vd,2 , . . . }. This construction is regarded as
breaking pieces off a unit-length stick successively with size determined by variable Vd,k ∈ [0, 1] drawn from Beta(1, γd ). The k-th
broken piece shows the proportion πd,k .
Note that the present model can be considered as an HMM withuniform transition probabilities. Future work will include the introduction of the hierarchical DP [11] to incorporate transition probabilities in order to more accurately model the succession of spectral patterns. Also note that although the number of components D
needs here to be given by hand, it should ideally be determined adaptively according to the input data for example using the techniques
described in [12, 13].
4. VARIATIONAL INFERENCE ALGORITHM
Inference algorithms for models with DP prior are mostly based on
sampling methods. However, a variational Bayesian (VB) approach
is often preferable for large-scale problems. Even though we target ashort sequence of music, its spectrogram comprises a large numberof parameters. We thus prefer here the framework of VB inference.
Some approaches have been proposed for variational inference
involving DP [9, 14]. In the context of VB, DP is usually approximated by finite symmetric Dirichlet prior (FSD) or truncated stickbreaking construction (TSB) [15]. In this paper, we will discuss only
the TSB for the sake of brevity. TSB is obtained by using a fixedvalue K and setting ∀d, p(Vd,K = 1) = 1; this implies that the
mixture proportions πd,k are equal to zero for k > K. Note that the
proposed model does have “infinite-state” spectral bases, and that thefinite approximation is only linked to the use of variational inference.
As a matter of fact, we do not have to make such an approximation
if we use sampling methods, such as Gibbs sampling.
The VB approach in general assumes a factorized form for theposterior distribution. That can be regarded as assuming that the parameters are independent of each other. Here, because of the strong
impact of ı d on Zd , it seems difficult to make such an assumption without impeding the performance of the inference. As shownby [14], a collapsed approach based on integrating out ı d can beapplied to overcome this problem, and we thus use collapsed variational Bayes instead of standard VB. Here, we integrate out ı d andthen the joint collapsed model is given by
P (Y , C, H, U , Z, W )= p(Y | C)p(C | H, U , Z)p(Z)p(H)p(U , W )(9)where p(Z) is given by (marginalizing out ı d from Eq. (8))
P(k )Y Y Γ(1 + n(k))
d )Γ(γd +k >k nd.
p(Z) =P(k )Γ(1 + γd + k ≥k nd )
d
k

(10)

We want to maximize the log marginal likelihood L(Y ). As it isdifficult to do it directly, VB uses a lower bound B(Y ) for L(Y )based on the assumption that parameters are independent of eachother. The lower bound is given by
XZ
P (Y , C, H, U , Z, W )B(Y ) =
Q(Ξ) log
dCdHdU dW
Q(Ξ)
Z

(TGSWGPE[=M*\?

(TGSWGPE[=M*\?



























6KOG=U?

(b) Estimated component: F(TGSWGPE[=M*\?

(TGSWGPE[=M*\?

(a) Original spectrogram











6KOG=U?

















6KOG=U?



6KOG=U?

(c) Estimated component: A(c) Estimated component: DP(k)k E[Hω,d ]q(Zd,t = k)E[Ud,t ] (not E[Cω,t,d ]) of each component (b),Fig. 1. Original spectrogram of vocal signals (a), estimated model(c) and (d). We set the number of bases to D = 3 and truncated level to K = 30. Owing to parsimony and temporal smoothness, the proposed
method was able to decompose a note with vibrato in an unsupervised way.
where Q(Ξ) = q(C)q(Z)q(H)q(U )q(W ). In the following, forthe sake of brievity, we only sketch the update rules. Let us first
focus on the update rules for q(Z). We find the update:“”q(Zd,t ) ∝ exp EQm=t q(Zd,m ) [log p(Zd,t | Z d,−t )]„«× exp Eq(ΘZ )q(C ) [log p(C t,d | ΘZd,t )](11)

We can then find the following update rule for q(H):(k)

q(Hω,d )

∝

(k)
μω,d

=

(k)

(Z

)

=

bω +

“
”
(k)
≈ log 1 + E[nd,−t ] −+X(

2(1 + E[nkd,−t ])2

j<k

−

X

(

(>j)

V[nd,−t ]

j≥k

(≥j)

V[nd,−t ]

(13)

ω,k

τd,t =X(k)

E[Hω,d ]q(Zd,t = k) + αE[Wd,t ] + βE[Wd,t−1 ] ,ω,kand finally for q(W ):(14)

ϕd,t = αE[Ud,t ] + βE[Ud,t+1 ] .

5. EXPERIMENTS)

2
2(1 + γd + E[n≥jd,−t ]).

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 16kHz. The magnitude
spectrogram was computed using the short time Fourier transformwith 32 ms long Hanning window and with 16 ms overlap.
The treatment of the hyperparameters is important. One possible
approach is introducing hyperpriors (such as non-informative priors)on the hyperparameters, which we expect would lead to automatic
estimation of the hyperparameters. However, we could instead indicate our wish to fix them so that the model is enforced by sparsenessand temporal smoothness. Here, the hyperparameters were set toaω = bω = 0, α = β = 5, γ = 1.
At first, we generated synthetic data, a mixture of vocal signals
taken from RWC database (RWC-MDB-I-2001 No.45) [16]. The
sequence is composed of 3 notes (Major chord: D，F and A have

As for the second term of Eq. (11), we findEq(ΘZ )q(C ) [log p(C t,d | ΘZd,t )]
d,t
X
X
(Z
)
(Z
)
=−
E[Hω,dd,t Ud,t ] +E[Cω,t,d log Hω,dd,t Ud,t ] .
ωsimilarly for q(U ):q(Ud,t ) ∝ Gamma(ηd,t , τd,t ) ,X
E[Cω,t,d ]q(Zd,t = k) + α + β ,ηd,t =

To avoid unexpected local solutions, we propose introducing aweight parameter Υ into p(C | H, U , Z)Υ , and gradually increasing the value of Υ to 1.

)

2
2(γd + E[n>jd,−t ])”“(≥j)
log 1 + γd + E[nd,−t ] −q(Zd,t = k)E[Ud,t ] ,tφd,t = α + β ,(k)V[nd,−t ]”“(>j)
log γd + E[nd,−t ] −Xq(Wd,t ) ∝ Gamma(φd,t , ϕd,t ) ,EQm=t q(Zd,m ) [log p(Zd,t | Z d,−t )]

(12)

t

(k)
νω,d

d,t

where ΘZd,t = {hd d,t , Ud,t }. Integrating out parameters often results in expensive computational cost. This can however be avoided
by using a Gaussian approximation [14]. This approximation can(k)(≥k)
be applied to random variables nd,−t , nd,−t which are sums overBernoulli variables. Thanks to the central limit theorem, these sums
can be effectively approximated using Gaussian distributions. In order to apply this framework, we use the following second order Taylor expansion: E(log ρm ) ≈ log(E[ρm ]) − V[ρm ]/2E[ρm ]2 for any
probabilistic variable ρm . We have as for the first term of Eq. (11)(k)

Gamma(μω,d , νω,d ) ,Xaω +E[Cω,t,d ]q(Zd,t = k) ,ωNext, we find the update for q(C):q(Cω,t,d ) ∝ Multinomial(Cω,t,d | Yω,t , λω,t,d ) ,
”
“P
(k)
exp
k q(Zd,t = k)E[log Hω,d Ud,t ]“P” .
λω,t,d = P
(k)
d exp
k q(Zd,t = k)E[log Hω,d Ud,t ]1974

(TGSWGPE[=M*\?

gorithm based on collapsed variational inference, and presented experimental results showing that our model is fitted for the modelingof non-stationary audio signals. In the future, we will extend thismodel to an infinite-state Markov chain spectrum model with hierarchical DP [11]. A remaining important problem for audio signal
analysis based on NMF is that of the automatic determination of thenumber of components [12, 13]: future work will include the introduction of an Indian buffet process prior into our model to overcome
this problem.







6KOG=U?

(TGSWGPE[=M*\?

(a) Original spectrogram, a mixture of piano (D) and violin (A)


7. REFERENCES



5VCVG

6KOG=U?








6KOG=U?

(TGSWGPE[=M*\?

(b) Estimated model: Violin (A)







5VCVG

6KOG=U?








6KOG=U?

(c) Estimated model: Piano (D)Fig. 2. Original spectrogram (a), estimated model (b) and (c). EstiP(k)mated model: k E[Hω,d ]q(Zd,t = k)E[Ud,t ] (top) and q(Zd,t =k)E[Ud,t ] (bottom) of each component. We set the number of bases
to D = 2 and truncated level to K = 30. The number of states was
automatically optimized depending on the instrumental sounds.

overlapping harmonic components): first, each note is played alone
in turn, then all the combinations of two notes are played and finally all notes are played simultaneously. “synthetic” means thateach note is played four times and they are synthesized from thesame signal. The result is shown in Fig. 1. The proposed method
was confirmed to appropriately decompose a note with vibrato with
deformable bases.
Next, we generated as test data (Fig. 2), a mixture of the piano(D: RWC-MDB-I-2001 No. 1) and the violin (A: RWC-MDBI-2001 No. 15). Each note is played alone in turn, then two notes
are played simultaneously. Each note is played twice and each time
synthesized from a different signal. The results are shown in Fig. 2.
The number of states is determined automatically depending on the
instruments.
6. CONCLUSION
We presented a nonparametric Bayesian extension of NMF, which
we call infinite-state spectrum model. The proposed model applied
to audio signals represents the time-varying spectrum of each note
event in polyphonic music. We derived an efficient optimization al-

1975[1] D. D. Lee and H. S. Seung, “Learning the parts of objects by
non-negative matrix factorization,” Nature, vol. 401, pp. 788–
791, Oct. 1999.
[2] Paris Smaragdis, “Non-negative matrix factor deconvolution;
extraction of multiple sound sources from monophonic inputs,”
in Proc. ICA, 2004, pp. 494–499.
[3] R. Hennequin, R. Badeau, and B. David, “NMF with timefrequency activations to model non stationary audio events,”
IEEE Trans. on Audio, Speech, and Language Process., 2010.
[4] A. Ozerov, C. Févotte, and M. Charbit, “Factorial scaled hidden Markov model for polyphonic audio representation andsource separation,” in Proc. WASPAA, 2009.
[5] M. Nakano, J. Le Roux, H. Kameoka, Y. Kitano, N. Ono, and
S. Sagayama, “Nonnegative matrix factorization with markovchained bases for modeling time-varying patterns in music
spectrograms,” in Proc. LVA/ICA, 2010, pp. 149–156.
[6] C. Févotte and A. T. Cemgil, “Nonnegative matrix factorizations as probabilistic inference in composite models,” in Proc.
EUSIPCO, 2009, vol. 47, pp. 1913–1917.
[7] O. Dikman and A. T. Cemgil, “Unsupervised single-channel
source separation using Bayesian NMF,” in Proc. WASPAA,
2009, pp. 93–96.
[8] C. Antoniak, “Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems,” Annals of Statistics, vol. 2, no. 6, pp. 1152–1174, 1974.
[9] D. M. Blei and M. I. Jordan, “Variational inference for Dirichlet process mixtures,” Journal of Bayesian Analysis, vol. 1, no.
1, pp. 121–144, 2005.
[10] J. Ishwaran and L. James, “Gibbs sampling methods for stickbreaking priors,” Journal of the American Statistical Association, pp. 161–174, 2001.
[11] Y. Teh, M. Jordan, M. Beal, and D. Blei, “Hierarchical Dirichlet processes,” Advances in neural information processing systems, vol. 17, pp. 1385–1392, 2004.
[12] V. Y. F. Tan and C. Fevotte, “Automatic relevance determination in nonnegative matrix factorization,” in Proc. SPARS,
2009.
[13] M. Hoffman, D. Blei, and P. Cook, “Bayesian nonparametricmatrix factorization for recorded music,” in Proc. ICML, 2010.
[14] K. Kenichi, MaxWelling, and Y. Whye Teh, “Collapsed variational dirichlet process mixture models,” in Proc. IJCAI, 2007.
[15] J. Sethuraman, “A constructive definition of Dirichlet priors,”
Statistica Sinica, vol. 4, pp. 639–650, 1994.
[16] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC
music database: Popular, classical, and jazz music database,”in Proc. ISMIR, 2002, pp. 287–288.