SEPARATION OF A MONAURAL AUDIO SIGNAL INTO HARMONIC/PERCUSSIVECOMPONENTS BY COMPLEMENTARY DIFFUSION ON SPECTROGRAM
Nobutaka Ono, Kenichi Miyamoto, Jonathan Le Roux, Hirokazu Kameoka, and Shigeki Sagayama
Department of Information Physics and Computing,
Graduate School of Information Science and Technology, The University of Tokyo
7-3-1 Hongo Bunkyo-ku, Tokyo, 113-8656, Japan
Email: {onono, miyamoto, leroux, kameoka, sagayama}@hil.t.u-tokyo.ac.jp
web: http://hil.t.u-tokyo.ac.jp/index-e.html

ABSTRACTIn this paper, we present a simple and fast method to separate a monaural audio signal into harmonic and percussivecomponents, which is much useful for multi-pitch analysis,
automatic music transcription, drum detection, modificationof music, and so on. Exploiting the differences in the spectrograms of harmonic and percussive components, the objectivefunction is defined in a quadrature form of the spectrogram
gradients. Applying the auxiliary function approach to that,
simple and fast update equations are derived, which guarantee the decrease of the objective function at each iteration.
We show some experimental results by applying our method
to popular and jazz music songs.
1. INTRODUCTION
Recently, music signal has become a significant target in thesignal processing field and various tasks have been discussed
like audio music retrieval, audio onset detection, multiple
fundamental frequency estimation, and so on [1]. The music signal often consists of two different components: harmonic one and percussive one. The simultaneous presence ofthem makes some tasks harder because of their much different spectral structures. For instance, most of the multipitch
analysis are interfered by percussive tones, while suppressionof harmonic components will facilitate the drum detection or
the rhythm analysis.
In this paper, aiming to efficient pre-processing of musicsignal analysis, we present a simple and fast method to separate a monaural audio signal into harmonic and percussive
components. This kind of separation problem has been discussed in several pilot research. Uhle et al. applied Independent Component Analysis (ICA) to the power spectrogram,and classified the extracted independent components into aharmonic and a percussive groups based on the several features like percussiveness, noise-likeness, etc [2]. Helen andVirtanen utilized Non-negative Matrix Factorization (NMF)
for decomposing the spectrogram into elementary patterns
and classified them by pre-trained Support Vector Machine(SVM) [3]. Through modeling harmonic and inharmonic
tones on spectrogram, Itoyama et al. proposed separationof an audio signal to each track part based on the MIDI information synchronized to the input audio signal [4]. Daudet
reviewed recent separation algorithms of music signals intosteady and transient components [5]. Other kinds of singlechannel separation problem have been recently developed in[6, 7, 8, 9]. The contribution of this paper is to derive a simple
and fast algorithm specifically for the harmonic/percussive
separation based on the anisotropy of them on spectrogram.

We present the formulation of the separation as an optimization problem, derive the fast iterative solution to it by auxiliary function approach, and examine the performance byexperiments to popular and jazz music songs.
2. FORMULATION OF HARMONIC/PERCUSSIVESEPARATION
Let Fh,i be a Short Time Fourier Transform (STFT) of amonaural audio signal f(t) and Wh,i = |Fh,i|2 be its power
spectrogram, where h and i represent indices of frequency
and time bins. A typical spectrogram of an audio signal isshown in Fig. 1. In it, the vertical and horizontal structures
are clearly observed. The harmonic component usually has astable pitch and forms parallel ridges with smooth temporalenvelopes on the spectrogram, while the energy of a percussive tone is concentrated in a short time frame, which forms
a vertical ridge with a wideband spectral envelope. Exploiting the anisotropy, we decompose the original power spectrogram Wh,i into the harmonic component Hh,i and the percussive component Ph,i on the spectrogram. For evaluating
their anisotropic smoothness, L2 norm of the power spectrogram gradients is examined here, that is, Hh,i and Ph,i arefound by minimizing
J(H, P ) =1 (Hh,i−1 − Hh,i )22
2σHh,i1 + 2(Ph−1,i − Ph,i )22σP

(1)

h,i

under the constraint as
Hh,i + Ph,i = Wh,iHh,i ≥ 0, Ph,i ≥ 0,(2)(3)

where H and P are sets of Hh,i and Ph,i , respectively,and σH and σP are parameters to control the weights of
the horizontal and vertical smoothness. Minimizing eq. (1)is equivalent to maximum likelihood estimation under theassumption that the spectrogram gradients (Hh,i−1 − Hh,i)and (Ph−1,i − Ph,i ) follow independent Gaussian distributions. Although the actual distributions of the spectrogram
gradients are different from them, the assumption leads us
to simple and comprehensive formulation and solution. As
confirmed later, replacing the power spectrogram Wh,i by
the range-compressed version: W̃h,i = |Fh,i |2γ (0 < γ ≤ 1)
partially bridges a gap between the assumption and the realsituation.

Then, the following updates:{U (k+1), V (k+1) } = min Q(H (k), P (k), U, V ),U ,V(11){H (k+1) , P (k+1) } = min Q(H, P , U (k+1), V (k+1) ), (12)H ,Pdecrease J monotonically, where k represents the number of
iterations, and U and V are sets of Uh,i and Vh,i , respectively.
4. DERIVATION OF UPDATE RULESFirst, we here derive H (k+1) and P (k+1) satisfying eq. (12)under the constraint of eq. (2). With introducing Lagrange
multipliers λh,i , consider
Figure 1: A spectrogram of a popular music song

Q̃(H, P ) = Q(H , P , U (k+1) , V (k+1))+λh,i (Hh,i + Ph,i − Wh,i ).

3. DESIGN OF AUXILIARY FUNCTION
Since the objective function defined by eq. (1) is quadratic
for all variables, it is monomodal and has a single global
minimum, which can be directly obtained by solving∂J/∂Hh,i = 0 and ∂J/∂Ph,i = 0. But they yield much largesize simultaneous equations of Hh,i and Ph,i (the order of theequations is equal to the number of all time-frequency binsof (h, i)). To avoid it and derive a simple iterative solution,
we adopt an auxiliary function approach, which is used inother signal processing methods as NMF [10] or HarmonicTemporal Clustering (HTC) [11].
In order to design the auxiliary function of our problem,
note that(A − B)2 ≤ 2(A − X)2 + 2(B − X)2(4)for any A, B, and X since2(A − X)2 + 2(B − X)2 − (A − B)22A+B= 4 X−2(13)h,i(5)

is obviously nonnegative and equal to zero where X = (A +B)/2. Applying the inequality to eq. (1), we have

Differentiating it by Hh,i , Ph,i , and λh,i yields, respectively2(k+1)(k+1)) + λh,i = 0,
2 (2Hh,i − Uh,i+1 − Uh,iσH2(k+1)(k+1)(2Ph,i − Vh+1,i − Vh,i ) + λh,i = 0,
σP2Hh,i + Ph,i − Wh,i = 0.

(14)
(15)
(16)

Solving them, we obtainα (k+1)
(k+1)
(k+1)
Hh,i=+ Uh,i )(U2 h,i+1(1 − α)(k+1)(k+1)+(2Wh,i − Vh+1,i − Vh,i ), (17)2
(1 − α) (k+1)
(k+1)
(k+1)
=
(Vh+1,i + Vh,i )Ph,i2α
(k+1)
(k+1)
(18)
+ (2Wh,i − Uh,i+1 − Uh,i ),2whereσ2α= 2 P 2 .
(19)
σH + σP(Hh,i−1 − Hh,i)2 ≤ 2(Hh,i−1 − Uh,i )2 + 2(Hh,i − Uh,i )2 ,(6)222
(Ph−1,i − Ph,i) ≤ 2(Ph−1,i − Vh,i ) + 2(Ph,i − Vh,i ) ,(7)While, the auxiliary parameters U (k+1) and V (k+1) satisfying eq. (11) are easily given byfor any Uh,i and Vh,i . The equalities are valid for Uh,i =(Hh,i−1 + Hh,i)/2 and Vh,i = (Ph−1,i + Ph,i )/2. Hence, theauxiliary function:(20)
By substituting eq. (20) into eq. (17) and eq. (18), we can remove the auxiliary parameters Uh,i and Vh,i from the update
rules as

Q(H, P , U, V )1 (Hh,i−1 − Uh,i )2 + (Hh,i − Uh,i )2= 2
σH(k)
(k)
Hh,i−1 + Hh,i(k+1),Uh,i =

2h,ih,isatisfies
J(H, P )J(H, P )≤=Q(H, P , U, V ),min Q(H, P , U, V ).
U ,V(9)(10)

whereΔ(k)= Hh,i + Δ(k),

(21)

(k+1)
Ph,i(k)Hh,i − Δ(k),(22)= α2(k+1)

Hh,i
1 + 2
(Ph−1,i − Vh,i )2 + (Ph,i − Vh,i )2 (8)
σP(k)
(k)
Ph−1,i + Ph,i(k+1)=Vh,i.

(k)

(k)

=(k)

(k)

Hh,i−1 − 2Hh,i + Hh,i+1−(1 − α)

4(k)
(k)
(k)
Ph−1,i − 2Ph,i + Ph+1,i4


.

(23)

Introducing a process to restrict the obtained solution tosatisfy eq. (3) and binarizing a separation result, which isexperimentally confirmed to yield better separation performance, the separation algorithm is consequently summarized
as follows.
1. Calculate Fh,i , the STFT of an input signal f(t).
2. Calculate a range-compressed version of the power spectrogram by
Wh,i = |Fh,i|2γ(0 < γ ≤ 1).

(24)

3. Set initial values as1(0)(0)
Hh,i = Ph,i = Wh,i ,2(25)for all h and i and set k = 0.
4. Calculate the update variables Δ(k) defined as eq. (23).
5. Update Hh,i and Ph,i as(k+1)= min(max(H h,i + Δ(k), 0), Wh,i), (26)

(k)

(k+1)

= Wh,i − Hw,iHh,i(k+1)

Ph,i.

(27)

6. Increment k. If k < kmax −1 (kmax : the maximum number of iterations), then, go to step 4, else, go to step 7.
7. Binarize the separation result as(k)(k)(Hh,imax , Ph,imax )(k=−1)(k−1)(0, Wh,i) (Hh,imax< Ph,imax ).
(k−1)(k−1) (28)
≥ Ph,imax )(Wh,i , 0) (Hh,imax(k)(k)

8. Convert Hh,imax and Ph,imax into waveforms by(k)h(t) =ISTFT[(Hh,imax )1/2γ ej Fh,i ],(29)p(t) =

(k)ISTFT[(Ph,imax )1/2γ ej Fh,i ],(30)

where ISTFT represents the inverse STFT.
Since the update variable Δ consists of the second orderderivative of Hh,i and Ph,i , the update in step 4 has basically
dfd2 fthe same form as the diffusion equation:= C 2 , whichdtdx
represents dynamics of diffusion phenomena like heat. Butunlike the physical diffusion process, each update of Hh,iand Ph,i includes a negative diffusion term derived fromthe other. As iterations, the energy distribution of Hh,i on
the spectrogram diffuses horizontally and concentrates vertically. Ph,i follows the inverse way. Both of which hold that
Hh,i + Ph,i = Wh,i , Hh,i ≥ 0, and Ph,i ≥ 0. We denote thediffusion-like process of two energy distributions with a balance by complementary diffusion. The balance parameter α(0 < α < 1) controls the strength of the diffusion along thevertical and the horizontal directions.
5. EXPERIMENTAL EVALUATIONS
We here show several results by experiments. The target signals were chosen from the RWC Music Database (PopularMusic Database) [12] and the sampling frequency was converted into 16kHz. The frame length of STFT is 1024 andFigure 2: The spectrograms of the iteratively-updated har(k)monic component Hh,i (left) and the percussive component(k)Ph,i (right) at k = 0, k = 3, k = 10, k = 50, and after binarization, from top to bottom, respectively.
the frame shift is 512. The balance parameter α was experimentally determined to be 0.3. The resultant spectrograms of(k)the harmonic component Hh,i and the percussive component(k)Ph,i to 6.25s fragment of RWC-MDB-P-2001 No.7 is shownin Fig. 2, where γ = 0.3 was used. We can see that the energyof spectrogram is splitting to two components as iterations,each of which is forming horizontal and vertical ridges, respectively. The computational time for a 6.25s-length signal
with 50 iterations is about 2.3s at a laptop-PC with 1.20GHz
Pentium in our implementation, which is nearly three times
faster than real-time processing. Indeed, we have developed
the real-time processing version of it.
In auditory evaluation, pitched instrument tracks anddrum tracks are well separated into h(t) and p(t), respectively. Nevertheless, the duration of some percussion, typically bass drum was almost separated into h(t) since it hasa smooth temporal envelope, too. Inversely, the attack ofpitched tone had a tendency to belong to p(t). A singing

energy ratio of harmonic component

1
0.9
0.8
0.7
0.6

piano
bass
synthesizer
bell 1snare drum
bass drum
high-hat
melody
brass

0.5
0.4
0.3
0.2
0.1
00

0.1

0.2

0.3

0.4

0.5
gamma

0.6

0.7

0.8

0.9

1

Figure 3: The energy ratio of the separated harmonic component in each track of RWC-MDB-P-2001 No.18
1
energy ratio of harmonic component

voice is also difficult to be perfectly classified to h(t) due tothe nature of time-varying pitch. In the condition, most ofthem except pitch-varying component belonged to h(t). We
think that utilizing wavelet transform instead of STFT hassome potential to improve the performance for pitch-varying
components. The resultant sound will be presented in theconference.
In order to quantitatively evaluate the separation performance of our algorithm, we prepared each track of audio
signals in two songs (RWC-MDB-P-2001 No.18 and RWCMDB-J-2001 No.16) by MIDI-to-WAV conversion. Let fi (t)be one audio track signal where i is an index of the track.
and inputed the summation of all tracks to our algorithm.
Then, the energy ratio of each track included in h(t) and p(t)was calculated as rh = Eh /(Eh + Ep ), rp = Ep /(Eh + Ep ),where Eh =< fi (t), h(t) >2 , Ep =< fi (t), p(t) >2 , and <>represents the cross correlation operation. The energy ratio of the harmonic component, rh in each track is shown in
Fig. 3 and Fig. 4, where several γs were examined. The
tendency that pitched instrumental tracks and the bass drum
track belong to h(t) and other percussion tracks to p(t) is
consistent to the auditory evaluation. We can see that theseparation performance gets worse as γ is larger, especiallyfor drum tracks. Thus, the range compression facilitates theseparation and γ  0.3 seems to give the best separation performance for these two songs. Probably, it is related to theactual distribution of spectrogram gradients of harmonic andpercussive components.

0.9
0.8
0.7
piano
bass
snare drum
high-hat 1
high-hat 20.60.5
0.4
0.3
0.2
0.1

6. CONCLUSION
In this paper, we derived a simple and fast method for amonaural audio signal into harmonic and percussive components, which is performed by the complementary diffusion
on spectrogram. Despite of the simplicity of the algorithm,
pitched instruments and drums are well separated. Applying the proposed algorithm as a pre-process to multi-pitch
analysis or other tasks of music signal analysis is our current
concern.
REFERENCES[1] http://www.music-ir.org/mirex2007/index.php[2] C. Uhle, C. Dittmar and T. Sporer, “Extraction of drum
tracks from polyphonic music using independent subspace analysis,” Proc. ICA, pp. 843–847, Apr. 2003.
[3] M. Helen and T. Virtanen, “Separation of drums frompolyphonic music using non-negative matrix factorization and support vector machine,” Proc. EUSIPCO,
Sep. 2005.
[4] K. Itoyama, M. Goto, K. Komatani, T. Ogata, and H.
Okuno, “Integration and Adaptation of Harmonic andInharmonic Models for Separating Polyphonic MusicalSignals,” Proc, ICASSP, pp. 57–60, Apr. 2007.
[5] L. Daudet, “A Review on Techniques for the Extraction of Transients in Musical Signals,” Proc. CMMR,
pp. 219-232, 2005.
[6] D. P. W. Ellis and R. J. Weiss, “Model-based monaural source separation using a vector-quantized phasevocoder representation,” Proc. ICASSP, pp. 957-960,
May 2006.

00

0.1

0.2

0.3

0.4

0.5
gamma

0.6

0.7

0.8

0.9

1

Figure 4: The energy ratio of the separated harmonic component in each track of RWC-MDB-J-2001 No.16[7] P. Smaragdis, B. Raj, and M. Shashanka, “Supervised and Semi-supervised Separation of Sounds fromSingle-Channel Mixtures,” Proc. ICA pp. 414-421, Sep.
2007.
[8] J. Le Roux, H. Kameoka, N. Ono, A. de Cheveigne,
S. Sagayama, “Single Channel Speech and Background
Segregation through Harmonic-Temporal Clustering,”
Proc. WASPAA, Oct, 2007.
[9] K. Miyamoto, H. Kameoka, T. Nishimoto, N. Ono, andS. Sagayama, “Harmonic-Temporal-Timbral Clustering (HTTC) for the Analysis of Multi-Instrument Polyphonic Music Signals, Proc. ICASSP2008, pp. 113116, Apr. 2008.
[10] D. D. Lee and H. S. Seung, “Algorithms for NonNegative Matrix Factorization” Proc. NIPS, pp. 556–562, 2000.
[11] H. Kameoka, T. Nishimoto, S. Sagayama, “A Multipitch Analyzer Based on Harmonic Temporal Structured Clustering,” IEEE Trans. on Audio, Speech andLanguage Processing, vol. 15, no. 3, pp.982-994, Mar.
2007.
[12] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka,“RWC music database: Popular, classical, and jazz music databases,” Proc. ISMIR, pp. 287-288, Oct. 2002.