EXPLICIT BEAT STRUCTURE MODELING FOR NON-NEGATIVE MATRIXFACTORIZATION-BASED MULTIPITCH ANALYSIS
Kazuki Ochiai, Hirokazu Kameoka, Shigeki Sagayama
Graduate School of Information Science and Technology, The University of Tokyo
Hongo 7-3-1, Bunkyo, Tokyo, 113-8656, Japan
ABSTRACTThis paper proposes model-based non-negative matrix factorization(NMF) for estimating basis spectra and activations, detecting note
onsets and offsets, and determining beat locations, simultaneously.
Multipitch analysis is a process of detecting the pitch and onsetof each note from a musical signal. Conventional NMF-based approaches often lead to unsatisfactory results very possibly due to thelack of musically meaningful constraints. As music is highly structured in terms of the temporal regularity underlying the onset occurrences of notes, we use this rhythmic structure to constrain NMF byparametrically modeling each note activation with a Gaussian mixture and derive an algorithm for iteratively updating model parameters. It is experimentally shown that the proposed model outperforms
the standard NMF algorithms as regards onset detection rate.
Index Terms— Polyphonic pitch transcription, Non-negative
matrix factorization, Rhythmic/Beat structure, Onset detection
1. INTRODUCTION
Music transcription is a process of obtaining a symbolic representation (such as a set of MIDI messages or a score) of an audio signal.
While there are a number of viable solutions for transcribing monophonic music, polyphonic music still poses a formidable challenge.
Most existing methods for polyphonic pitch transcription rely
on prior knowledge about the sources contained in the polyphonicdata being analyzed. The main weakness of these kinds of methods
is that they lack the capacity to adapt to signals that do not comply
with the assumption they make about the sources. On the other hand,relatively recent techniques based on sparse representations use as
few hypotheses as possible about the audio content to separate thenotes. The goal of this kind of approach is to ﬁnd a set of basis spectra such that any observed spectrum can be concisely represented
as a linear combination of a small number of ‘active’ basis spectra.
One successful approach involves applying the Non-negative Matrix Factorization (NMF) to a power spectrogram (a time-frequency
representation) interpreted as a non-negative matrix [1]. In this approach, a spectrogram Y is factorized into a product of two factors
with non-negative entries, one being a basis matrix H, consistingof basis spectra, and the other an activation matrix U , consisting oftime-varying amplitudes associated with the basis spectra.
One way of obtaining a symbolic representation from a polyphonic audio signal is to apply NMF to its spectrogram and thenperform onset/offset detection in the activations associated with thebasis spectra [2]. This transcription method is based on the assumption that each basis spectrum obtained with NMF corresponds to asingle pitch. We should thus be able to determine note onsets andoffsets by simply thresholding the lines of the activation matrix U .
In practice, however, the lack of constraints in the NMF model often
leads to unsatisfactory results. One reason is that the envelope of

978-1-4673-0046-9/12/$26.00 ©2012 IEEE133

each basis activation typically contains many noisy peaks and dipsdue to the mismatch between the actual source spectra and the basis spectra. These kinds of errors occurring in the NMF phase willpropagate through to the onset/offset detection phase, thus making
it difﬁcult to reach correct decisions. To obtain the basis activity information with as few errors as possible, it is necessary to incorporate an appropriate constraint into the NMF model. Many attempts
have already been made to develop constrained variants of the NMF
model: “Sparse NMF” promotes the sparsity of the basis activations
while performing decomposition [3]. “Non-negative matrix factor
deconvolution (NMFD)” [4] and “non-negative matrix factor 2D deconvolution (NMF2D)” [5] use spectro-temporal signatures rather
than the basis spectra to represent spectrograms. “Smooth NMF”
enforces temporal continuity on the basis activations [6]. “Harmonic
NMF” imposes harmonicity constraints on the basis spectra [7]. The“Non-negative hidden Markov model” uses Markov-chained basis
spectra to represent time-varying patterns in spectrograms [8–11].
In our view, transcription methods that consist of performing
NMF or one of its variants mentioned above for the preprocessing,and then performing onset/offset detection for the postprocessing
will at a certain point face limitations in terms of pitch transcription
performance. This is because the estimation of basis activations andthe detection of note onsets and offsets are each a prerequisite of theother. If we were able to obtain the correct basis activations, then it
would be a relatively simple matter to detect the onsets and offsets of
the underlying notes. On the other hand, if we were provided withthe onsets and offsets of the notes, they could constitute very useful information for estimating the basis activations. Therefore, thisleads to a “chicken and egg” situation. Furthermore, as the onsets
of notes are governed by the rhythmic structure of a piece of music,
the “chicken and egg” situation also applies to the detection of note
onsets and the determination of beat locations. If we knew the beatlocations of a piece of music, then it would be much easier to detect
note onsets from the basis activations, and vice versa. Because on
this, we consider it necessary to introduce a uniﬁed model for multipitch analysis, which could be used to jointly solve the problemsof determining the basis spectra and activations, detecting the note
onsets and offsets, and determining the beat locations.
2. PROPOSED MODEL2.1. Non-negative Matrix Factorization model
Let us start our discussion with the standard NMF modelXω,t =

K

Hω,k Uk,t ,

(1)

k=1

where k is the index of each basis spectrum, ω and t are frequency and time indices, respectively. The set consisting of
H1,k , . . . , HΩ,k represents the k-th basis spectrum and Uk,t is

ICASSP 2012

Power density

where vk,n is the total energy of object n, and τk,n is the center of theﬁrst Gaussian, which can be considered an estimate of the onset time.
The centers of the Gaussians are constrained to be equally spacedwith the distance φk,n , which is equal to the “standard deviation” of
all the Gaussians. This speciﬁc constraint allows the entire object tobe stretched or shrunk linearly in the time direction according to thevalue of φk,n . Thus, φk,n can be regarded as a parameter related tothe duration of a note. wk,n,1 , . . . , wk,n,M are weights associatedwith the M Gaussians, which determine the shape of the object. To
avoid an indeterminacy in the scaling of vk,n and wk,n,m , we assume

Time

Pitch

∀k,n :

M

wk,n,m = 1.

(6)

m=1

Time

Fig. 1. An illustration of the parametric activation model, Vk,n,t ,consisting of M = 7 Gaussians (top) and along with the beat locations (bottom). The vertical dashed lines and the dotted lines represent multiples and fractions of the beat periods.
the activity of the k-th basis spectrum at time t. Given an observed spectrogram Y = (Yω,t )Ω×T , the goal of NMF is to ﬁnd
H = (Hω,k )Ω×K ∈ R≥0,Ω×K and U = (Uk,t )K×T ∈ R≥0,K×T
such that Y  X where X = (Xω,t )Ω,T = HU . We assumeHω,k = 1 (k = 1, . . . , K),(2)
ω

in order to avoid an indeterminacy in the scaling of H and U .
2.2. Explicit beat structure constraint
Music is highly structured in terms of the temporal regularity underlying the onset occurrences of notes. In general, the time between
consecutive onsets corresponds to multiples and fractions of the beatperiod, with small deviations in timing and tempo.
Based on the rhythmic structure of music, we make the following assumptions to constrain the activity function Uk,t :1. Each activity function consists of local activity patterns,
called “objects”, each of which we expect to correspond to asingle note activation.
2. Each object is characterized by a fast/slow rise at the onset
time followed by a continuous contour.
3. The onset of each object is likely to be located at multiples or
fractions of the beat period.
4. The beat period varies gradually over time.
First, from assumption 1, Uk,t can be expressed as
Uk,t =Nk

Vk,n,t ,(3)

n=1

where Vk,n,t denotes the n-th object and Nk is the number of objects
in the k-th activity function. To incorporate assumption 2 into Vk,n,t ,we introduce a parametric model from [12], which is expressed as asum of Gaussians (Fig. 1):
Vk,n,t =M

Gk,n,m,t ,

(4)

m=1

vk,n wk,n,m −(t−(m−1)φk,n −τk,n )2 /2φ2k,ne,
Gk,n,m,t = √2πφk,n(5)Now, recall that our goal is to achieve Y  X. We thus
need to deﬁne a discrepancy measure between these two quantities.
Let D(yx) be a discrepancy measure between y and x such thatD(yx) ≥ 0 and D(yx) = 0 only if y = x. We can then deﬁnea goodness-of-ﬁt measure between the observed spectrogram Y andthe current NMF model X

D(Yω,t Xω,t ).
(7)
J (Y X) =ω,tAs in [13], we choose to use the I-divergence as the distortion measure D(yx)D(yx) = y log(8)

Minimizing J (Y X) with respect to X is then known to be equivalent to maximizing the Poisson likelihood. The likelihood of theunknown parameters, H, v, w, τ , φ, given Y is thereforep(Y |H, v, w, τ , φ) =Poisson(Yω,t ; Xω,t ),

(9)

ω,t

where Poisson(y; x) = xy e−x /y!.
To avoid overﬁtting the shape of each object, we place a
Dirichlet prior over wk,n = {wk,n,m }1≤m≤M , namely p(w) = xi −1
k,n Dirichlet(w k,n ; α), where Dirichlet(y; x) ∝i yiMand α = (α1 , . . . , αM ). αm / m αm =1 is the expected valueof wk,n,m . To enforce the sparseness of the activity function,
we place a generalized Gaussian prior over vk,n , as in [14],pp(v) = k,n GN (vk,n ; 0, λ, p), where GN (y; 0, λ, p) ∝ e−λ|y| .
When 0 < p < 2 this distribution becomes super-Gaussianand promotes the sparsity of active objects. Henceforth, we
assume 0 < p ≤ 1. Furthermore, we would like to ensurethat each basis spectrum has a harmonic structure of a particular pitch. To do so, we shall place a Gamma prior over H,
namely p(H) =ω,k Gamma(Hω,k ; γ H̄ω,k + 1, β), whereGamma(x; a, b) ∝ xa−1 e−bx , and its mode is given by H̄ω,k .
H̄ω,k is thus the most likely value for Hω,k and β determines thepeakiness of the density around the mode.
Next, to impose Assumption 3, we ﬁrst introduce a set of hyperparameters, ψ = {ψd }0≤d≤D , where ψd corresponds to the timeinterval between the (d + 1)-th and d-th beat locations. With these
hyperparameters, we can design a Gaussian prior distribution over
the onset parameter τk,np(τ |ψ) =


k,n

134

y
− (y − x).
x

N (τk,n ; ρn , ν 2 ),

(10)

ρn =

d−1


ψl +

l=0
2

i−1
ψd ,
I

(11)

2

where N (x; μ, σ 2 ) ∝ e−(x−μ) /2σ , ρn denotes the most expected
location of the onset of the n-th object, ν 2 is the variance of theGaussian indicating how much τk,n is allowed to deviate from ρn , I
is the number of divisions per beat, and the indices d, i are such thatn = (d − 1)I + i. To impose Assumption 4, we place a Gaussian
chain hyperprior over ψp(ψ) = p(ψ0 )Dp(ψd |ψd−1 ),

(12)

d=1

p(ψd |ψd−1 ) = N (ψd ; ψd−1 , σ 2 ).

(13)

The next step is to derive update equations for the unknown parameters using the auxiliary function. Here, we must take account of thenormalization conditions for Hω,k and wk,n,m . As for Hω,k , below
we consider a simple updating procedure, which consists of solving
the unconstrained maximization of the auxiliary function and projecting its solution onto the constrained space. The update equations
are as follows:
p−1),
vk,n ← Ak,n /(1 + λpηk,nω,t Wk,ω,t Gk,n,m,t + αm − 1
wk,n,m ←,Ak,n + m (αm − 1)−bk,n + bk,n 2 + 4Ak,n ck,n,φk,n ←2Ak,nν 2 Bk,n + φ2k,n {Ψd + (i − 1)ψd /I},ν 2 Ak,n + φ2k,n
Iσ 2 i,k (i−1)(τk,n−Ψd )+I 2ν 2 (ψd−1 +ψd+1 )ψd ←,Kσ 2 i (i − 1)2 + 2I 2 ν 2Yω,t Uk,t /Xω,t + β H̄ω,k
Hω,k ← Hω,k t,t Uk,t + β

τk,n ←

3. PARAMETER ESTIMATION
3.1. Maximum a posteriori (MAP) estimation problem
Given an observed spectrogram Y , we would like to ﬁnd the estimates of H, v, w, τ , φ, ψ that maximize the log posterior densityp(H, v, w, τ , φ, ψ|Y ) ∝ p(Y |H, v, w, τ , φ)p(H)p(v)p(w)p(τ |ψ)p(ψ). We therefore consider the problem of maximizing
L(θ) = log p(Y |H, v, w, τ , φ) + log p(H)+log p(v) + log p(w) + log p(τ |ψ) + log p(ψ)

where Wk,ω,t = Yω,t Hω,k /Xω,t , Ψd =Ak,n =

subject to∀k :Hω,k = 1, ∀k,n :
uk,n,m = 1, ∀d : ψd ≥ 0,
ω3.3. Parameter updates

Bk,n =

ω,tmbk,n =∀ω,k : Hω,k ≥ 0, ∀k,n,m : vk,n , wk,n,m , φk,n , τk,n ≥ 0,

3.2. Designing an auxiliary function
When applying an auxiliary function approach to a maximization
problem, the ﬁrst step is to deﬁne a lower bound function for the objective function. The difﬁculty with the current maximization problem lies in the terms log p(Y |H, v, w, τ , φ) and log p(v). We use
the following inequalities to bound these terms from below:
log p(Yω,t |H, v, w, τ , φ)Hω,k Gk,n,m,t
γk,n,m,ω,t log− Xω,t − log(Yω,t !),≥ Yω,tγk,n,m,ω,t
k,n,m

p−1
log p(vk,n ) ≥ −λ(pηk,n(vk,n − ηk,n ) + ηk,n ) + log

pλ1/p,2Γ(1/p)for vk,n > 0, in which the exact bounds are achieved when
γk,n,m,ω,t = 

Hω,k Gk,n,m,tk ,n ,m Hω,k Gk ,n ,m ,t

ηk,n = vk,n .









,

(14)
(15)

The auxiliary function for our objective function will thus be deﬁned
by replacing the terms log p(Yω,t |H, v, w, τ , φ) and log p(vk,n )with the lower bounds shown above.

135

ck,n =



(18)
(19)
(20)
(21)

l=0 ψl and(22){t − (m − 1)φk,n } Gk,n,m,t ,

(23)

m



(t − τk,n )Wk,ω,t(m − 1)Gk,n,m,t ,ω,twhere θ denotes a set consisting of H, v, w, τ , φ, ψ. Although it
is difﬁcult to solve this optimization problem analytically, we candevelop a computationally efﬁcient algorithm for ﬁnding a locally
optimal solution based on the auxiliary function concept, similar tothe one used in [12, 13].

Wk,ω,t(17)

d−1

Wk,ω,t Vk,n,t ,

ω,t

(16)

(24)

m


(t − τk,n )2 Wk,ω,t Vk,n,t .

(25)

ω,t

Based on the auxiliary function principle, the objective function is
non-decreasing under the updates (14), (15) and (16)–(21), and so
the convergence of the iterative procedure is guaranteed if each iteration involves these updates.
4. EXPERIMENTAL RESULTS
To verify the performance of our method, we evaluated the onset detection rate using a few recordings found in the RWC music database[15]. The data for each evaluation were the ﬁrst 23 s, mixed down toa monaural signal and resampled to 16 kHz. The constant-Q transform was used to compute spectrograms where the time resolution,
the lower bound of the frequency range, and the frequency resolution were set at 16 ms, 60 Hz and 12 cents, respectively. D and theinitial value of ψd were set at the values obtained with [16]. K,
I, N , and M were set at 74, 4, D × I, and 20, respectively. ν,σ, λ were all set at 1. p, β, and αm were set at 0.5, 1000, and1000 × e−m/3 / m e−m /3 , respectively. The initial values ofmmvk,n , wk,n,m , φk,n , and τk,n were set at 20, e− 3 / m e− 3 , 5,d−1
and l=0 ψl + (i−1)ψd /I, respectively. The initial values of Hω,kand H̄ω,k were set at the value obtained with the standard NMF applied to the piano excerpts from the RWC musical instrument sounddatabase [17]. The algorithm was run for 300 iterations, after which
any object whose energy vn,k /M φn,k was less than a certain threshold was considered to be silent.

C6

Table 1. Onset detection performance of the standard NMF and proposedmodels. F, P and R correspond to F -measure, precision, and recall, respectively (%).

Note

C5
C4

Data and Instruments

Notes

C3

C-2001, 30, Piano
C-2001, 35, Piano
J-2001, 1, Piano
J-2001, 2, Piano
J-2001, 6, Guitar
J-2001, 7, Guitar
J-2001, 8, Guitar
J-2001, 9, Guitar

119
50
155
69
161
99
79
94

C2
3

6

9

Time [s]1215

standard NMF
F
P
R
18.9
13.7
30.2
10.3
5.8
44.0
15.0
8.5
59.4
11.8
6.8
43.5
14.8
9.0
41.6
12.3
8.2
24.2
7.5
4.1
41.8
11.7
7.1
33.0

proposed model
F
P
R
68.8
66.4
71.4
64.4
72.5
58.0
71.3
62.7
82.6
64.9
60.8
69.6
59.3
58.9
59.6
44.0
48.2
40.4
30.7
25.2
39.2
51.8
49.5
54.3

this approach properly yields the onset times and the beat locations,and also provides higher pitch estimation accuracy than the standard
NMF. Our future work will include employing time-varying spectral
basis patterns to deal with the timbre change of the instruments [10],improving the estimation of the beat locations by incorporating models of musical rhythm and implementing a transcription application
using this model.

Fig. 2. An estimated piano roll (top) and the reference pitch data (bottom)
for Chopin’s Nocturne No. 2 in E-maj, Op. 9. The vertical dashed lines represent the estimated beats.

6. REFERENCES

with the present method applied to Chopin’s Nocturne No. 2 inE-maj, Op. 9 (RWC-MDB-C-2001 No. 30). We compared the proposed method with the standard NMF by employing the F-measure.
We used the same reference pitch data as in [12]. The number ofnotes for the proposed method was simply the number of objects.
As for the standard NMF, Hω,k was initialized in the same manneras the proposed method, and after convergence the k-th basis spectrum was considered to be turned on at time t when Uk,t exceededa threshold and turned off when it was below this threshold. The
threshold for each method varied from zero to the maximum valueof the activations, and the maximum F -measure was selected. The
results are shown in Tab. 1. It is worth noting that it is harder toobtain a high detection rate for an onset detection task than a frameby-frame pitch detection task. To illustrate this, consider a task thatinvolves detecting the pitch of a note with a duration of 100 frames.
Suppose that method A was able to detect all the frames correctly
except for one in the middle of the duration. The frame-by-frame
pitch detection rate in this case will be 99%. On the other hand, theonset detection rate will be only 50% (one correct onset at the starting frame and one inserted onset in the middle of the duration). As
can be seen from Tab. 1, the standard NMF model provided poor
results, implying that the basis activations obtained with the standard NMF had redundant peaks and dips. On the other hand, theproposed method obtained a signiﬁcantly higher detection rate especially for the piano recordings. The results for the guitar recordings
were not as high as for the piano recordings because the initial basis
H and the mode of the Gamma prior were set at the spectra of piano
notes. If we were also provided with a set of basis spectra learned
from other instruments, we would be able to obtain a higher detection rate for polyphonic music played on many kinds of instruments.
Furthermore, the initial setting of the estimates of the beat locations
is an important issue, that must be solved in the future, since it became clear that the estimates of the beat locations were sensitive toinitialization.

[1] P. Smaragdis and J.C. Brown, “Non-negative matrix factorization for
polyphonic music transcription,” in Proc. WASPAA, 2003, pp. 177–180.
[2] N. Bertin, R. Badeau, and G. Richard, “Blind signal decompositions
for automatic transcription of polyphonic music: NMF and K-SVD on
the benchmark,” in Proc. ICASSP, 2011, vol. 1, pp. 65–68.
[3] S.A. Abdallah and M.D. Plumbley, “Unsupervised analysis of polyphonic music by sparse coding,” IEEE Trans. on Neural Networks, vol.
17, no. 1, pp. 179–196, 2006.
[4] P. Smaragdis, “Non-negative matrix factor deconvolution; extraction ofmultiple sound sources from monophonic inputs,” in Proc. ICA, 2004,
pp. 494–499.
[5] M. N. Schmidt and M Mørup, “Sparse non-negative matrix factor 2D deconvolution for automatic transcription of polyphonic music,” inProc. ICA, 2006, pp. 700–707.
[6] T. Virtanen, “Monaural sound source separation by nonnegative matrix
factorization with temporal continuity and sparseness criteria,” IEEE
Trans. on Audio, Speech, and Language Processing, vol. 15, no. 3, pp.
1066–1074, 2007.
[7] S. A. Raczyński, N. Ono, and S. Sagayama, “Multipitch analysis withharmonic nonnegative matrix approximation,” in Proc. ISMIR, 2007,
pp. 381–386.
[8] A. Ozerov, C. Févotte, and M. Charbit, “Factorial scaled hidden
Markov model for polyphonic audio representation and source separation,” in Proc. WASPAA, 2009, pp. 121–124.
[9] M. Nakano, J. Le Roux, H. Kameoka, Y. Kitano, N. Ono, and
S. Sagayama, “Nonnegative matrix factorization with Markov-Chainedbases for modeling time-varying patterns in music spectrograms,” inProc. LVA/ICA, 2010, pp. 149–156.
[10] M. Nakano, J. Le Roux, H. Kameoka, T. Nakamura, N. Ono, and
S. Sagayama, “Bayesian nonparametric spectrogram modeling basedon inﬁnite factorial inﬁnite hidden Markov model,” in Proc. WASPAA,
2011, pp. 325–328.
[11] G. J. Mysore, P. Smaragdis, and B. Raj, “Non-negative hidden Markov
modeling of audio with application to source separation,” in Proc.
LVA/ICA, 2010, pp. 140–148.
[12] H. Kameoka, T. Nishimoto, and S. Sagayama, “A multipitch analyzer
based on harmonic temporal structured clustering,” IEEE Trans. on
Audio, Speech, and Language Processing, vol. 15, no. 3, pp. 982–994,
2007.
[13] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in Proc. NIPS, 2000, pp. 556–562.
[14] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF:
A new sparse representation for acoustic signals,” in Proc. ICASSP,
2009, pp. 3437–3440.
[15] 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.
[16] D. P. W. Ellis, “Beat tracking by dynamic programming,” Journal of
New Music Research, vol. 36, no. 1, pp. 51–60, 2007.
[17] M. Goto, “Development of the RWC music database,” in Proc. the 18th
International Congress on Acoustics (ICA 2004), 2004, pp. I–553–556.

5. CONCLUSION
In this paper, we proposed a new framework for the multipitch analysis of polyphonic music signals based on NMF. By incorporating
the object model into the activation matrix, the onsets and the beatstructure, which are very important for a proper audio-to-score music transcription, can be obtained simultaneously, together with thespectral basis matrix. We performed experiments that showed that136