MULTIPITCH ANALYSIS WITH HARMONIC NONNEGATIVE MATRIXAPPROXIMATION
Stanisław A. Raczyński Nobutaka Ono Shigeki Sagayama
The University of Tokyo
Graduate School of Information Science and Engineering
E-mail: {raczynski,onono,sagayama}@hil.t.u-tokyo.ac.jp
ABSTRACT

tonal music). It is necessary to distinguish between thefundamental tones and their overtones.

This paper presents a new approach to multipitch analysis
by utilizing the Harmonic Nonnegative Matrix Approximation, a harmonically-constrained and penalized version of the Nonnegative Matrix Approximation (NNMA)
method. It also includes a description of a note onset, offset and amplitude retrieval procedure based on that technique. Compared with the previous NNMA approaches,
specific initialization of the basis matrix is employed – thebasis matrix is initialized with zeros everywhere but at positions corresponding to harmonic frequencies of consequent notes of the equal temperament scale. This results inthe basis containing nothing but harmonically structured
vectors, even after the learning process, and the activity
matrix’s rows containing peaks corresponding to note onset times and amplitudes. Furthermore, additional penalties of mutual uncorrelation and sparseness of rows areplaced upon the activity matrix. The proposed method
is able to uncover the underlying musical structure better
than the previous NNMA approaches and makes the note
detection process very straightforward.

A large variety of methods has been used to tackle themultipitch analysis problem (e.g. [1], [2], [4], [9], [11],[12], [13]; an exhaustive list of methods would be very
long and we are not going to include it here, but for agood summary, see [6]), but so far none of them solving the problem in a satisfactorily precise and universal
way. While our lab has recently developed a powerful
method based on Harmonic Temporal Structured Clustering (HTC) for this purpose [4], the procedure proposed inthis paper is built upon a method from the family of Nonnegative Matrix Approximations (NNMA), which, under
different names and in different varieties, has recently received much attention, also from the music transcription
community. To the best of our knowledge, however, noneof the NNMA-based methods were developed specificallyfor analysis of musical signals. As it will be shown later inthis paper, nature of music can be exploited to increase thetranscription potential of the algorithm. The goal of thispaper was to propose a NNMA variation most suitable formultipitch analysis.

1 INTRODUCTION
The problem of automatic polyphonic music transcription(extracting underlying musical structure from sampled music) has been addressed numerous times, and it still seems
there is a long way to go before arriving at a robust anduniversal technique. This paper tries to lay another bricktowards this goal.
Automatic music transcription of recorded music is usually a two stage process. The first stage is the event detection phase, where music events (note onsets, note offsets,
pitch changes) are detected and identified. In the second
stage, these events are transformed into a musical score.
This paper focuses on the event detection stage, main partbeing multipitch analysis, which aims to uncover the fundamental frequencies of simultaneously played harmonicsounds. It is a difficult task, since each sound, besides
the fundamental tone, consists of many harmonic tones,
some of them having the same frequencies as the fundamental frequencies of other sounds (e.g. in the case of
c 2007 Austrian Computer Society (OCG).

Different variations and extensions of the NNMA algorithm have been used for multipitch analysis: the regular
NNMA [13], its penalized versions, such as the Nonnegative Sparse Coding (NNSC) [2, 1], or NNMA with basis
vectors extended to contain spectrotemporal signatures (anumber of consequent data frames), such as Nonnegative
Matrix Factor 2-D Deconvolution (NMF2D) and Sparse
Nonnegative Matrix Factor 2-D Deconvolution (SNMF2D)[11]. These methods have, however, a few drawbacks.
They do not guarantee to yield basis vectors with harmonic structure. NMF2D and SNMF2D use a single signature for every note (of a single instrument), making useof the shift-similarity of logarithmic frequency scale spectra of notes played on a single instrument. This might
be an oversimplification resulting in an inadequate model.
The note spectra are similar, but not identical, with significant differences for some specific instruments (e.g. flute).
Moreover, spectrotemporal atoms cannot account for different note lengths, which results in multiple activity peaks
when notes are longer than the signature, and lower activity peaks when notes are shorter than the signature. Allof the previous work published on that subject do not propose a complete transcription procedure, simply reporting

good results after visual comparison of the activities andsymbolic data used to generate the analyzed music.
The paper is organized as follows. Section 2 presents
a theoretical introduction to the Nonnegative Matrix Approximation and its extension through the addition of constraints and penalties placed upon both the basis matrix
and the activity matrix. An overview of the proposed
procedure in given in section 3, including description of
the proposed Harmonic Nonnegative Matrix Approximation (HNNMA) technique (3.3) and note detection method(3.4). The procedure is evaluated and compared with theresults of regular NNMA methods in section 4.

the Bregman divergence becomes the I-divergence (generalized Kullback-Leibler divergence):
DKL (P, Q) = P ⊙ log

P−P+Q ,
Q(8)

where the logarithm and the division are calculated elementwise. This situation leads to the simple NNMA, also known
as Nonnegative Matrix Factorization (NMF) [8]. Lee andSeung in [8] has proposed a very fast algorithm for minimizing the I-divergence that can be derived by using auxiliary functions. A function G(P, P′ ) is an auxiliary function for function F (P) if:
1. G(P, P) = F (P),2 THEORETICAL BACKGROUND

2. G(P, P′ ) ≥ F (P).

2.1 Definitions and basic properties
For clarity, the following notation was used in this paper:· is a sum of all the elements of a matrix, ⊙ is the
Hadamard product (calculated element-wise) and 1 is amatrix (of appropriate dimensions) containing nothing but
ones.
A few easy to prove properties were later used. If A ∈RN ×M , then:∇A A ⊙ A = 2A,(1)∇A AT A = 2AT 1,(2)∇A BA = BT 1,(3)

Making use of the convexity of ϕ [14], it can be shown
that:
G(A, A′ ) = ϕ(X) + X̃ − X



A
X
′
− ′ ⊙ A′ ⊙ log ′ S + ϕ(X̃ ) ,
A
X̃

(9)

e ′ = A′ S, is an auxiliary function forwhere XF (A)e
= DKL (X, AS) = DKL (X, X)
X= X ⊙ log − X + X̃
X̃

(10)

∇A f1 (A) + f2 (A) = ∇A f1 (A) + ∇A f2 (A) , (4)∇A B ⊙ (Cυ(A)) = υ ′ (A) ⊙ CT B ,
(5)

and

2.2 Generalized Nonnegative Matrix Approximatione ′ = AS′ , is an auxiliary function forwhere this time X

+

′

+

where υ : R → R is an element-wise function and υ
is its derivative, and f1 , f2 : R+,N ×M → R+,N ×M are
any matrix functions.

Generalized Nonnegative Matrix Approximation (describedin [3] and later developed in [14]), is a method for decomposition of a nonnegative (having only nonnegative elements) matrix X (later referred to as the data matrix) intoa multiplication of two, also nonnegative, matrices A andS (later refereed to as the basis matrix and the activity matrix, respectively):eX∼= AS = X.

(6)

D(P, Q) = ϕ(P) − ϕ(Q) − ϕ′ (Q) ⊙ (P − Q) ,(7)

The Generalized NNMA solves this problem by minimizing a Bregman divergence between the data matrix X
e A Bregman divergence between
and its approximation X.
two matrices is defined aswhere ϕ : S ⊆ R → R is a strictly convex function with
continuous first derivative, calculated here for each element of a matrix separately. If ϕ(p) = p log p − p, thene −X
G(S, S′ ) = ϕ(X) + X
 


S
X
′
′
e
⊙ A S ⊙ log ′ + ϕ(X ) , (11)−e′S
X
F (S) = DKL (X, AS).

(12)

It can also be shown [14] that F (S′ ) is non-increasing under the update
S′ ← arg min G(S, S′ ).

(13)

S

To solve this optimization problem, we calculate thegradient of the auxiliary function and force it to zero:e −X∇S G(S, S′ ) = ∇S ϕ(X) + X


 
X
S
e ′ ) = 0. (14)−⊙ A S′ ⊙ log ′ + ϕ(Xe′SXUsing properties (3), (4) and (5), we can easily calculatethat gradient as:

 
S
X
⊙ A S′ ⊙ log ′
= 0,
(15)
∇S AS −
e′
S
X

 

S′
X
T
= 0,
A 1−
⊙ A
e′
S
X
 
X
AT X
e′
.
S = S′ ⊙AT 1This suggest a multiplicative update rule:

X
AT AS
S←S⊙
.
AT 1
T

(16)

3.1 Overview of the procedure
(17)

(18)

We can come up with a similar update rule for the basis
matrix A:∇A G(A, A′ ) = 0,
(19)


′
X T
A
⊙
= 0,
(20)
S
1ST −
e′
A
X

X T
S
(21)
A ← A ⊙ AS T .
1SSo, by using simple multiplicative rules from (21) and(18), we can find matrices A and S that minimize the Idivergence from equation (8).

2.3 Penalized NNMA
Because H(P, P′ ) = G(P, P′ )+α(P) is an auxiliary function for F (P) + α(P) (the proof is straightforward), the
NNMA algorithm can be extended by placing additional
penalties upon both estimated decomposition matrices, expressed as additional element in the objective function.
∇A H(A, A′ ) = ∇A G(A, A′ ) + ∇A α(A),
(22)


A′
X T
∇A H(A, A′ ) = 1ST −+ ∇A α(A).
⊙
S
A
X′
(23)
While it is very difficult to solve this non-linear equation with respect to A, the following approximation canbe used [14]:∇A α(A) ∼= ∇A α(A)

,

(24)

A=A′

which is asymptotically true, as difference between A inconsequent iterations tends to 0. Now:


X T
A′
+ ∇A′ α(A′ ) = 0,
(25)
S
⊙
1ST −
A
X′
which yields an update rule:
A←A⊙

X T
AS S

1ST + ∇A α(A).

3 MULTIPITCH ANALYSIS PROCEDURE(26)

Similarly:

The NNMA algorithm decomposes the data matrix X,
which does not contain musical data, but rather some midlevel representation of it. In most cases its columns are
power spectra of consecutive frames of time-domain musical data. In the proposed procedure a constant-Q transform is used. The central frequencies of the constant-Q
filters can be set to correspond to the frequencies of thenotes of the most common twelve-tone equal temperament (12-TET) scale or can further divide each semitone,
which is a very useful property for analyzing musical signals. What is more, the dimensionality of a constant-Qtransformed data is much lower than the dimensionalityof a Fourier-transformed data, which makes the computation of the NNMA faster. After calculating the constant-Q
transform, the resulting data is fed through the HNNMA
algorithm, which decomposes it to a product of the basis
matrix and activity matrix. Activity matrix is analyzed in
the last part of the procedure – the note detector, describedin section 3.4.
3.2 Matrix initialization in HNNMA
Because zero-valued elements of basis vectors will remain
zero-valued throughout the learning process, we could initialize them to have zeros everywhere but at the positionsof fundamentals of notes from a specific range of the 12TET scale and their harmonics. Furthermore, that would
guarantee that the basis vectors are sorted by their fundamental frequencies, and that corresponding rows in theactivity matrix contain activities of consequent notes fromthat range, resulting in a harmonicaly-constrained NNMA.
This would make analysis of the results of the algorithm
straightforward – one would only have to analyze the note
activities and find peaks corresponding to instances of thesenotes.
In the proposed procedure, after initialization, each basis vector is multiplied by the normalized mean value ofthe constant-Q transform of the data at the bin corresponding to this note. This should discourage the HNNMA fromlearning these notes and using them to reconstruct the analyzed data. During the learning process, each row of theactivity matrix is, as it is usually done in the learning process of the NNMA methods, normalized to unit squaredsum, while the basis matrix is simply normalized by its
maximal value to let the basis vectors, that correspond tonotes not existing in analyzed music, freely decrease.
3.3 Additional penalties in HNNMA∇S HS (S, S ) ∼= ∇S G(S, S′ ) + ∇S′ β(S′ ),
′

(27)
(28)

HNNMA extends the regular NNMA to include additional
penalties on the activity matrix S. We would like to findsuch an activity matrix that would:

It must be noted that the new update rules may result inthe matrices A and S becoming negative, so caution must
be taken while constructing the objective function.

1. be sparse, i.e. each row should contain only very
few non-zero elements (to reduce the low-valued
noise),

S←S⊙

X
AT AS

AT 1 + ∇S β(S)

.

d

e

f#

g#

a#

c’

d’

e’

f#’

g#’

Note activities

0

5

10

15

20

Time [s]Figure 1. Basis vectors after analysis of the Ode to Joy.
Its harmonic structure is clearly visible. Vertical dottedlines indicate expected harmonic peaks positions.
2. contain mutually uncorrelated rows (to reduce theinter-row crosstalk, like e.g. octave errors).
The above can be reformulated, accordingly, in termsof a objective function β:β(S)= −µ1 log(1 + S ⊙ S)+µ2

ST S − S ⊙ S .

(29)

The first element, log (1 + S ⊙ S) , is one of the often
used sparseness measures [5]. The second one is a measure of correlation between every pair of different matrix
rows:
XX
ST S − S ⊙ S =sTi sj ,
(30)
i

j6=i

where sTi is the i-th row and si is its transposition. Using
properties (1), (2) and (5) from section 2.1, we can easily
calculate the gradient:∇S β(S) = −2µ1 S/ (1 + S ⊙ S) + 2µ2 S (1 − I) . (31)
Similar penalties could be used for the basis matrix, but
our experiments with sparsity, column uncorrelation andcolumn shift-similarity showed that these constraints do
not improve procedure’s accuracy when the basis matrix
was initialized in the way described in the next subsection.
When the matrix was initialized with traditional noise, theconstraints would result in basis vectors containing peaks,
although the structure was not always purely harmonic.
Thus, either further, much more complex constraint are
required to enforce this structure, or we could take the advantage of the multiplicative nature of HNNMA algorithm
update rules.
3.4 Note detector
Before analysis, each row of the activity matrix is multiplied by the height of the peak at the fundamental frequency in corresponding basis vector (all rows are being
normalized to unit squared sum during the learning process). This should make the activities of notes that do

Figure 2. Note activities (S) after analysis of the Ode to
Joy. Peaks correspond to notes detected in the signal
not exist in the analyzed music significantly smaller than
the activities of notes that occur in the music. After thateach row with values higher than some arbitrarily chosen
threshold are normalized to the maximal value in each ofthem.
It turns out that in preliminary experiments the resulting activities clearly correspond to notes in the analyzed
musical piece. However, if the notes were played shortly
one after another, their peaks blend to form a single peakwith multiple sub-peaks. Because of that, a simple thresholding is not enough. A still simple, but much more robust thresholding method was used. First, the activities
are thresholded to detect peaks and blended peaks (e.g.
two blended peaks depicted on Figure 3). Then, for eachdetection all local maxima and local minima are found.
Some of the maxima correspond to actual sub-peaks, whilesome are just fluctuations in the note activity. Two thresholds are set between the highest local maximum and thelowest local minimum. All maxima that are above thehigher threshold (upper ligth-gray range on Figure 3) andhas at least one minimum lower than the lower threshold(lower light-gray range on picture 3) are marked as subpeaks and are assumed to correspond to individual notes.
In similar fashion, all minima under the lower thresholdthat lay between two sub-peaks are assumed to be the offset time of the note corresponding to the left sub-peak andonset time of the note corresponding to the right sub-peak.
The beginning and the end of the blended group of subpeaks are assumed to be the onset time of the first note in
the group and the offset time of the last note in the group.
The last step of the note detection process is acceptancedecision for each of the detected peaks and sub-peaks. A
peak is accepted and regarded a note only if its width multiplied by its height is greater than some threshold.
4 EXPERIMENTAL RESULTS
4.1 Experiment conditions
To validate our approach, we tested our procedure on a
few recordings. All analyzed recordings were played on

Figure 3. Results of note detection for a quarter note andan eighth blended together (example of real data). Darkgray circles mark the detected sub-peaks and the offset
time of the first note.
Composer
Title
Notes Acc. Corr.
L. Beethoven Symphony in D minor, 101 96% 86%
Op. 125, No. 9 (lastmovement, Ode to joy)F. Chopin
Nocturne in E# major, 328 87% 72%
Op. 9, No. 2 (part)F. Chopin
Nocturne in Bb minor, 358 70% 74%
Op. 9, No. 1 (part)J. S. Bach
Minuet No 4 in G
102 97% 100%

Table 1. Piano pieces used for algorithm evaluation
piano, as listed in Table 1, but experiments show equally
good results for acoustic guitar and violin (though lower
detection accuracy for violin).
During the experiments, all the parameters of the algorithm were kept constant in order to evaluate its robustness, though by fine-tuning the parameters for each musical piece separately, much better results can be obtained.
It has been noted that the best results were achieved byapplying this method for shorter blocks of data (30-60 s),
instead of the whole song at once.
The input data was first mixed down to a monaural signal and resampled to 11025 kHz. The constant-Q transform was calculated for frames shifted 12 ms. During
learning, µ1 = 1, µ2 = 10 were used. During note detection phase, before normalization of activities, rows having
maximal values lower than 0.125 of the maximal value ofthe activity matrix S were set to zero. The main detection
threshold was set at 0.25, the lower threshold to 0.25 ofdifference between the lowest minimum and the highest
maximum and the higher threshold to 0.75 of that difference.
After learning the basis matrix contained very well structured vectors, each one having a stronger peak for the fundamental tone and weaker peaks for the harmonics (Figure

Figure 4. Three bars from the middle of Chopin’s Nocturne in E# major, Op. 9, No. 2(a) Note activities obtained by NNMA(b) The first 4 bars of the played score, as a reference
Figure 5. Note activities after note detection for Ode to
Joy with grey squares being the original notes
1). The results of note detection for few example pieces
of music are presented in Table 1. Correctness of transcription is the ratio of the difference between the numberof notes in analyzed music and the number of deletions, tothe number of notes in analyzed music. Accuracy is the ratio of the difference between the number of detected notes
and the number of insertions, to the number of detected
notes. The first and the last musical piece were relatively
easy to analyze, containing notes from a rather short range(about 2 octaves). The two Chopin’s nocturnes were, on
the other hand, very difficult – played with a big dynamic
and wide range of note lengths, and containing notes fromwithin 5-6 octaves (see e.g. Figure 4).
4.2 Comparison with previous methods
Figure 6 depicts the basis matrix obtained using standard
NNMA (NMF) method after fundamental frequency estimation and basis vector sorting. It does not contain clear
harmonic structure – many of the vectors have two (ormore) dominant peaks, sometimes with highest peaks being the overtones instead of the fundamental, sometimes
having the same fundamental frequency as different basis
vectors (Figure 6). Slightly better results are obtained byutilizing different penalized NNMA methods proposed in
the literature (e.g NNSC [1] or Local Nonnegative Matrix Factorization [10]), however, the basis matrix nevercontains as highly harmonically structured vectors as theones obtained with the proposed method. The activities
obtained with these methods contain a lot of fluctuations
and assigning note names to the activities depends on thehighly dubious operation of fundamental frequency estimation of the basis vectors.

19Trans. on Neural Networks, vol. 17, no. 1, pp. 179–
196, 2006.

13
11
9

[4] Kameoka, H., Nishimoto, T., Sagayama, S. “Multipitch Analyzer Based on Harmonic Temporal Structured Clustering,” IEEE Trans. on Audio, Speech andLanguage Processing, vol. 15, no. 3, pp. 982–994,
Mar, 2007.

1

3

5

7

Basis vector nr

1517[3] Dhillon, I.S. and Sra, S. “Generalized Nonnegative
Matrix Approximations with Bregman Divergences,”
Proc. Neural Information Processing Systems, Vancouver, USA 2005.

0

50

100

150

Frequency bin

Figure 6. Basis obtained with standard NMF. Note detection relies on pitch estimation of these, often multipeaked, basis vectors.
It seems that the activity matrix in the proposed procedure contains very easy to analyze and at the same time
almost complete information about the underlying musical structure of the analyzed signal. This method is agood compromise between full basis estimation methods(such as NMF and other NNMA-based approaches) andmethods that use pre-learned basis vectors (e.g. [12] or[9]). The achieved results are similar to the results of different recently developed music transcription techniques(e.g. [4]), but by fine-tuning the method’s parameters,even greater accuracy could be achieved. The proposed
procedure uses a relatively simple method of analyzing
the activity matrix, making room for future research inmore advanced techniques, such as modeling the temporal envelopes of notes or using models of musical rhythm
and harmony.
5 CONCLUSION
In this paper, we discussed the use of Harmonic Nonnegative Matrix Approximation for multipitch analysis ofpolyphonic music signals. By initializing the basis matrix with harmonic structure and using new penalties ofsparsity and uncorrelation of rows of the activity matrix,
this approach yielded higher note detection accuracy compared with previous extensions of the Nonnegative Matrix
Approximation algorithm. The future work includes improving the post-processing of the HNNMA results by incorporating models of musical rhythm and harmonicity.
6 REFERENCES[1] Abdallah, S.A. and Plumbley, M.D. “Polyphonic music transcription by non-negative sparse coding of
power spectra,” Proc. 5th International Conference onMusic Information Retrieval, pp. 318–325, Barcelona,
Spain, 2004.
[2] Abdallah, S.A. and Plumbley, M.D. “Unsupervised
analysis of polyphonic music by sparse coding,” IEEE

[5] Karvanen, J. and Cichocki, A. “Measuring sparsenessof noisy signals,” Proc. 4th International Symposiumon Independent Component Analysis and Blind Signal
Separation, 2003.
[6] Klapuri, A.P. “Automatic Music Transcription as We
Know it Today,” Journal of New Music Research, vol.
33, no. 3, pp. 269–282, 2004.
[7] Lee, D.D. and Seung, H.S. “Algorithms for Nonnegative Matrix Factorization,” Advances in Neural
Information Processing Systems, vol. 13, pp. 556–562,
2001.
[8] Lee, D.D. and Seung, H.S. “Learning the parts ofobjects by nonnegative matrix factorization,” Nature,
vol. 401, no. 6755, pp. 788–791, 1999.
[9] Lepain, P. “Polyphonic Pitch Extraction from MusicalSignals,” Journal of New Music Research, vol. 28, no.
4, pp. 296–309, 1999.
[10] Li, S.Z. and Hou, X.W. and Zhang, H.J. and Cheng,
Q.S. “Learning spatially localized, parts-based representation,” Proc. IEEE Conf. Computer Vision and
Pattern Recognition, pp. 1–6, 2001.
[11] Schmidt, M.N. and Mørup M. “Sparse Non-negative
Matrix Factor 2-D Deconvolution for Automatic Transcription of Polyphonic Music,” Proc. 6th International Symposium on Independent Component Analysis and Blind Signal Separation, Charleston, USA,2006.
[12] Sha, F. and Saul, L.K. “Real-Time Pitch Determination of One or More Voices by Nonnegative Matrix
Factorization,” Advances in Neural Information Processing Systems, vol. 17, 2005.
[13] Smaragdis, P. and Brown. J.C. “Non-Negative Matrix Factorization for Polyphonic Music Transcription,” Proc. 2003 IEEE Workshop on Applications of
Signal Processing to Audio and Acoustics, New York,USA, 2003.
[14] Sra, S. and Dhillon, I.S. “Nonnegative Matrix Approximations: Algorithms and Application,” Technical Report Tr-06-27, Computer Sciences, Universityof Texas, Austin, USA 2006.