15th International Society for Music Information Retrieval Conference (ISMIR 2014)

HARMONIC-TEMPORAL FACTOR DECOMPOSITIONINCORPORATING MUSIC PRIOR INFORMATION FOR
INFORMED MONAURAL SOURCE SEPARATION
Tomohiko Nakamura† , Kotaro Shikata† , Norihiro Takamune† , Hirokazu Kameoka†‡†
Graduate School of Information Science and Technology, The University of Tokyo.
‡
NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation.
{nakamura,k-shikata,takamune,kameoka}@hil.t.u-tokyo.ac.jp

monicity and the coherences and continuities of amplitude and frequency modulations). This method is called“harmonic-temporal clustering (HTC).”
The other approach involves applying non-negative matrix factorization (NMF) to an observed magnitude spectrogram (time-frequency representation) interpreted as a
non-negative matrix [19]. The idea behind this approachis that the spectrum at each frame is assumed to be represented as a weighted sum of a limited number of common
spectral templates. Since the spectral templates and themixing weights should both be non-negative, this implies
that an observed spectrogram is modeled as the product of
two non-negative matrices. Thus, factorizing an observed
spectrogram into the product of two non-negative matrices allows us to estimate the unknown spectral templates
constituting the observed spectra and decompose the observed spectra into components associated with the estimated spectral templates.
The two approaches described above rely on diﬀerent
clues for making separation possible. Roughly speaking,the former approach focuses on the local time-frequency
structure of each source, while the latter approach focuses on a relatively global structure of music spectrograms (such a property that a music signal typically consists of a limited number of recurring note events). Rather
than discussing which clues are more useful, we believethat both of these clues can be useful for achieving a reliable monaural source separation algorithm. This belief has
led us to develop a new model and method for monauralsource separation that combine the features of both HTC
and NMF. We call the present method “harmonic-temporalfactor decomposition (HTFD).”
The present model is formulated as a probabilistic generative model in such a way that musically relevant information can be flexibly incorporated into the prior distributions of the model parameters. Given the recent progressof state-of-the-art methods for a variety of music information retrieval (MIR)-related tasks such as audio key detection, audio chord detection, and audio beat tracking, information such as key, chord and beat extracted from thegiven signal can potentially be utilized as reliable and useful prior information for source separation. The inclusionof auxiliary information in the separation scheme is referred to as informed source separation and is gaining increasing momentum in recent years (see e.g., among others, [5,15,18,20]). This paper further describes some ideas
how to design the prior distributions for the present model
to incorporate musically relevant information.
We henceforth denote the normal, Dirichlet and Poisson

ABSTRACT
For monaural source separation two main approaches havethus far been adopted. One approach involves applying
non-negative matrix factorization (NMF) to an observed
magnitude spectrogram, interpreted as a non-negative matrix. The other approach is based on the concept ofcomputational auditory scene analysis (CASA). A CASAbased approach called the “harmonic-temporal clustering(HTC)” aims to cluster the time-frequency components
of an observed signal based on a constraint designed according to the local time-frequency structure common in
many sound sources (such as harmonicity and the continuity of frequency and amplitude modulations). This paper proposes a new approach for monaural source separation called the “Harmonic-Temporal Factor Decomposition (HTFD)” by introducing a spectrogram model thatcombines the features of the models employed in the NMF
and HTC approaches. We further describe some ideas how
to design the prior distributions for the present model toincorporate musically relevant information into the separation scheme.
1. INTRODUCTION
Monaural source separation is a process in which the signals of concurrent sources are estimated from a monauralpolyphonic signal and is one of fundamental objectives offering a wide range of applications such as music information retrieval, music transcription and audio editing.
While we can use spatial cues for blind source separation with multichannel inputs, for monaural source separation we need other cues instead of the spatial cues.
For monaural source separation two main approaches havethus far been adopted. One approach is based on the concept of computational auditory scene analysis (e.g., [7]).
The auditory scene analysis process described by Bregman [1] involves grouping elements that are likely to haveoriginated from the same source into a perceptual structure called an auditory stream. In [8, 10], an attempt
has been made to imitate this process by clustering timefrequency components based on a constraint designed according to the auditory grouping cues (such as the harc Tomohiko Nakamura† , Kotaro Shikata† , Norihiro⃝
Takamune† , Hirokazu Kameoka†‡ .
Licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0). Attribution:
Tomohiko Nakamura† , Kotaro
Shikata† , Norihiro Takamune† , Hirokazu Kameoka†‡ . “Harmonictemporal factor decomposition incorporating music prior information for
informed monaural source separation”, 15th International Society forMusic Information Retrieval Conference, 2014.

623

15th International Society for Music Information Retrieval Conference (ISMIR 2014)
distributions by N, Dir and Pois, respectively.

at time t is expressed as a harmonically-spaced Gaussian
mixture function. If we assume the additivity of power
spectra, the power spectrogram of a superposition of K
pitched sounds is given by the sum of Eq. (8) over k. It
should be noted that this model is identical to the one employed in the HTC approach [8].
Although we have defined the spectrogram model abovein continuous time and continuous log-frequency, we actually obtain observed spectrograms as a discrete timefrequency representation through computer implementations. Thus, we henceforth use Yl,m := Y(xl , tm ) to denote an observed spectrogram where xl (l = 1, . . . , L) andtm (m = 1, . . . , M) stand for the uniformly-quantized logfrequency points and time points, respectively. We will
also use the notation Ωk,m and ak,n,m to indicate Ωk (tm ) andak,n (tm ).

2. SPECTROGRAM MODEL OF MUSIC SIGNAL
2.1 Wavelet transform of source signal model
As in [8], this section derives the continuous wavelet transform of a source signal. Let us first consider as a signal
model for the sound of the kth pitch the analytic signal
representation of a pseudo-periodic signal given byN∑
fk (u) =ak,n (u)ej(nθk (u)+φk,n ) ,(1)
n=1

where u denotes the time, nθk (u) + φk,n the instantaneous
phase of the n-th harmonic and ak,n (u) the instantaneous
amplitude. This signal model implicitly ensures not to violate the ‘harmonicity’ and ‘coherent frequency modulation’ constraints of the auditory grouping cues. Now, let
the wavelet basis function be defined by(u − t)1
ψα,t (u) = √
ψ
,
(2)
α
2πα
where α is the scale parameter such that α > 0, t the shiftparameter and ψ(u) the mother wavelet with the center frequency of 1 satisfying the admissibility condition. ψα,t (u)can thus be used to measure the component of period α attime t. The continuous wavelet transform of fk (u) is thendefined by∫ ∞∑N
Wk (log α1 , t) =ak,n (u)ej(nθk (u)+φk,n ) ψ∗α,t (u)du. (3)2.2 Incorporating source-filter model
The generating processes of many sound sources in realworld can be explained fairly well by the source-filter theory. In this section, we follow the idea described in [12] toincorporate the source-filter model into the above model.
Let us assume that each signal fk (u) within a short-time
segment is an output of an all-pole system. That is, if
we use fk,m [i] to denote the discrete-time representationof fk (u) within a short-time segment centered at time tm ,fk,m [i] can be described asP∑
βk,m [0] fk,m [i] =βk,m [p] fk,m [i − p] + ϵk,m [i], (9)−∞ n=1Since the dominant part of ψ∗α,t (u) is typically localizedaround time t, the result of the integral in Eq. (3) shalldepend only on the values of θk (u) and ak,n (u) near t. By
taking this into account, we replace θk (t) and ak,n (t) withzero- and first-order approximations around time t:ak,n (u) ≃ ak,n (t), θk (u) ≃ θk (t) + θ̇k (t)(u − t). (4)p=1

where i, ϵk,m [i], and βk,m [p] (p = 0, . . . , P) denote thediscrete-time index, an excitation signal, and the autoregressive (AR) coeﬃcients, respectively. As we have already assumed in 2.1 that the F0 of fk,m [i] is eΩk,m , to makethe assumption consistent, the F0 of the excitation signal
ϵk,m [i] must also be eΩk,m . We thus define ϵk,m [i] asN∑Ωk,mϵk,m [i] =vk,n,m e jne iu0 ,(10)Note that the variable θ̇k (u) corresponds to the instantaneous fundamental frequency (F0 ). By undertaking theabove approximations, applying the Parseval’s theorem,
and putting x = log(1/α) and Ωk (t) = log θ̇k (t), we canfurther write Eq. (3) asN∑
Wk (x, t) =ak,n (t)Ψ∗ (ne−x+Ωk (t) )e j(nθk (t)+φk,n ) , (5)n=1where u0 denotes the sampling period of the discrete-time
representation and vk,n,m denotes the complex amplitude of
the nth partial. By applying the discrete-time Fourier transform (DTFT) to Eq. (9) and putting Bk,m (z) := βk,m [0] −βk,m [1]z−1 · · · − βk,m [P]z−P , we obtain√N
2π ∑vk,n,m δ(ω − neΩk,m u0 ),(11)
Fk,m (ω) =Bk,m (e jω ) n=1

n=1

where x denotes log-frequency and Ψ the Fourier transformof ψ. Since the function Ψ can be chosen arbitrarily, aswith [8], we employ the following unimodal real function
whose maximum is taken at ω = 1:
 (log ω)2e− 4σ2(ω > 0) .
Ψ(ω) = (6)0(ω ≤ 0)Eq. (5) can then be written asN∑(x−Ωk (t)−log n)2
Wk (x, t) =ak,n (t)e− 4σ2e j(nθk (t)+φk,n ) .

where Fk,m denotes the DTFT of fk,m , ω the normalized angular frequency, and δ the Dirac delta function. The inverse
DTFT of Eq. (11) gives us another expression of fk,m [i]:N∑
Ωk,mvk,n,mfk,m [i] =e jne iu0 .
(12)
Ωk,m
jne
u0 )
n=1 Bk,m (eBy comparing Eq. (12) and the discrete-time representation of Eq. (1), we can associate the parameters of thesource filter model defined above with the parameters introduced in 2.1 through the explicit relationship:
vk,n,m|ak,n,m | =.
(13)Ω
Bk,m (e jne k,m u0 )

(7)

n=1

If we now assume that the time-frequency components aresparsely distributed so that the partials rarely overlap eachother, |Wk (x, t)|2 is given approximately asN∑(x−Ωk (t)−log n)2|Wk (x, t)|2 ≃|ak,n (t)|2 e− 2σ2
.
(8)
n=1

2.3 Constraining model parameters
The key assumption behind the NMF model is that the
spectra of the sound of a particular pitch is expressed as

This assumption means that the power spectra of the partials can approximately be considered additive. Note thata cutting plane of the spectrogram model given by Eq. (8)62415th International Society for Music Information Retrieval Conference (ISMIR 2014)a multiplication of time-independent and time-dependent
factors. In order to extend the NMF model to a more reasonable one, we consider it important to clarify which factors involved in the spectra should be assumed to be timedependent and which factors should not. For example, theF0 must be assumed to vary in time during vibrato or portamento. Of course, the scale of the spectrum should also be
assumed to be time-varying (as with the NMF model). On
the other hand, the timbre of an instrument can be considered relatively static throughout an entire piece of music.
We can reflect these assumptions in the present model inthe following way. For convenience of the following analysis, we factorize |ak,n,m | into the product of two variables,
wk,n,m and Uk,m√|ak,n,m | = wk,n,m Uk,m .
(14)
wk,n,m can be interpreted as the relative power of the nthharmonic and Uk,m as the time-varying normalized
ampli∑
tude of the sound of the kth pitch such that k,m Uk,m = 1.
In the same way, let us put vk,n,m as√
vk,n,m = w̃k,n,m Uk,m .
(15)

7000

Frequency [Hz]

3929

|n=1

{z

1238
695
390
0

0.37

0.73

1.1

1.46

1.83

Time [s]Figure 1. Power spectrogram of a violin vibrato sound.
p(w|w̃, β, Ω) expressed by the Dirac delta function)∏ (w̃k,n,m.
p(w|w̃, β, Ω) =δ wk,n,m −Ωu
Bk (e jne k,m 0 )
k,n,m

(20)

The conditional distribution p(w|β, Ω) can thus be obtainedby defining the distribution p(w̃) and marginalizing overw̃. If we now assume that the complex amplitude w̃k,n,mfollows a circular complex normal distributionw̃k,n,m ∼ NC (w̃k,n,m ; 0, ν2 ), n = 1, . . . , N,
(21)

Since the all-pole spectrum 1/|Bk,m (e jω )|2 is related to thetimbre of the sound of the kth pitch, we want to constrain
it to be time-invariant. This can be done simply by eliminating the subscript m. Eq. (13) can thus be rewritten asw̃k,n,mwk,n,m =.
(16)Ω
Bk (e jne k,m u0 )We can use Ωk,m as is, since it is already dependent on m.
To sum up, we obtain a spectrogram model Xl,m as
 N

K
∑
(x −Ω
−log n)2 ∑ 2− l k,m 22σXl,m =
Ck,l,m , Ck,l,m =  wk,n,m e Uk,m ,
k=1

2205

where NC (z; 0, ξ 2 ) = e−|z| /ξ /(πξ2 ), we can show, as in [12],that wk,n,m follows a Rayleigh distribution:
2

2

wk,n,m ∼ Rayleigh(wk,n,m ; ν/|Bk (e jne

Ωk,m
u

0

)|),

(22)

−z2 /(2ξ2 )where Rayleigh(z; ξ) = (z/ξ )e. This defines theconditional distribution p(w|β, Ω).
The F0 of stringed and wind instruments often varies
continuously over time with musical expressions such as
vibrato. For example, the F0 of a violin sound varies periodically around the note frequency during vibrato, as depicted in Fig. 1. Let us denote the standard log-F0 corresponding to the kth note by µk . To appropriately describe the variability of an F0 contour in both the global
and local time scales, we design a prior distribution forΩk := (Ωk,1 , Ωk,2 , . . . , Ωk,M )T by employing the productof-experts (PoE) [6] concept using two probability distributions. First, we design a distribution qg (Ωk ) describing how likely Ωk,1 , . . . , Ωk,L stay near µk . Second, we
design another distribution ql (Ωk ) describing how likelyΩk,1 , . . . , Ωk,L are locally continuous along time. Here we
define qg (Ωk ) and ql (Ωk ) as2}

Hk,l,m(17)
where Ck,l,m stands for the spectrogram of the kth pitch. If
we denote the term inside∑the parenthesis by Hk,l,m , Xl,mcan be rewritten as Xl,m = k Hk,l,m Uk,m and so the relation
to the NMF model may become much clearer.
2.4 Formulating probabilistic model
Since the assumptions and approximations we made so far
do not always hold exactly in reality, an observed spectrogram Yl,m may diverge from Xl,m even though the parameters are optimally determined. One way to simplify theprocess by which this kind of deviation occurs would be to
assume a probability distribution of Yl,m with the expectedvalue of Xl,m . Here, we assume that Yl,m follows a Poisson
distribution with mean Xl,m
Yl,m ∼ Pois(Yl,m ; Xl,m ),(18)

qg (Ωk ) = N(Ωk ; µk 1 M , υ2k I M ),ql (Ωk ) = N(Ω k ; 0 M , τ2k D−1 ),

(23)

(24)


 1 −1 0 0 · · · 0 −1 2 −1 0 · · · 0  0 −1 2 −1 · · · 0 

(25)
D =  ..
. . . . . . ..  , .
. . . .  0 · · · 0 −1 2 −1
0 · · · 0 0 −1 1
where I M denotes an M × M identity matrix, D an M × M
band matrix, 1 M an M-dimensional all-one vector, and 0 M
an M-dimensional all-zero vector, respectively. υk denotes
the standard deviation from mean µk , and τk the standard
deviation of the F0 jumps between adjacent frames. The
prior distribution of Ωk is then derived asp(Ωk ) ∝ qg (Ωk )αg ql (Ωk )αl(26)
where αg and αl are the hyperparameters that weigh thewhere Pois(z; ξ) = ξz e−ξ /Γ(z). This defines our likelihood
function∏
Pois(Yl,m ; Xl,m ),(19)p(Y|θ) =l,mwhere Y denotes the set consisting of Yl,m and Θ the entire
set consisting of the unknown model parameters. It should
be noted that the maximization of the Poisson likelihood
with respect to Xl,m amounts to optimally fitting Xl,m to Yl,mby using the I-divergence as the fitting criterion.
Eq. (16) implicitly defines the conditional distribution

625

15th International Society for Music Information Retrieval Conference (ISMIR 2014)
Ak = (Ak,1 , . . . , Ak,M )T . Here we introduce Dirichlet distributions:Ak ∼ Dir( Ak ; γ(A)R ∼ Dir(R; γ(R) ),
(28)
k ),
∏ ξi (A)(A)(A) T
where Dir(z; ξ) ∝ i zi , γk := (γk,1 , . . . , γk,M ) , and(R)
T
γ(R) := (γ1(R) , . . . , γ(R)K ) . For p(R), we set γk at a reasonably high value if the kth pitch is contained in the scale and(A)vice versa. For p( Ak ), we set γk,m
< 1 so that the Dirichlet
distribution becomes a sparsity inducing distribution.

contributions of qg (Ωk ) and ql (Ωk ) to the prior distribution.
2.5 Relation to other models
It should be noted that the present model is related to othermodels proposed previously.
If we do not assume a parametric model for Hk,l,m andtreat each Hk,l,m itself as the parameter, the spectrogram
model Xl,m can be seen as an NMF model with timevarying basis spectra, as in [14]. In addition to this assumption, if we assume that Hk,l,m is time-invariant (i.e.,
Hk,l,m = Hk,l ), Xl,m reduces to the regular NMF model [19].
Furthermore, if we assume each basis spectrum to havea harmonic structure, Xl,m becomes equivalent to the harmonic NMF model [16, 21].
If we assume that Ωk,m is equal over time m, Xl,m reduces
to a model similar to the ones described in [17, 22]. Furthermore, if we describe Uk,m using a parametric functionof m, Xl,m becomes equivalent to the HTC model [8, 10].
With a similar motivation, Hennequin et al. developed
an extension to the NMF model defined in the short-time
Fourier transform domain to allow the F0 of each basis
spectrum to be time-varying [4].

4. PARAMETER ESTIMATION ALGORITHM
Given an observed power spectrogram Y := {Yl,m }l,m ,we would like to find the estimates of Θ :={Ω, w, β, V, R, A} that maximizes the posterior densityp(Θ|Y) ∝ p(Y|Θ)p(Θ). We therefore consider the problem of maximizing
L(Θ) := ln p(Y|Θ) + ln p(Θ),(29)
with respect to Θ where∑()
ln p(Y|Θ) =Yl,m ln Xl,m − Xl,m
(30)
c

l,m

ln p(Θ) = ln p(w|β, Ω) +3. INCORPORATION OF AUXILIARYINFORMATION
3.1 Use of musically relevant information
We consider using side-information obtained with thestate-of-the-art methods for MIR-related tasks including
key detection, chord detection and beat tracking to assist
source separation.
When multiple types of side-information are obtained
for a specific parameter, we can combine the use of the
mixture-of-experts and PoE [6] concepts according to the“AND” and “OR” conditions we design. For example,
pitch occurrences typically depend on both the chord andkey of a piece of music. Thus, when the chord and key information are obtained, we may use the product-of-experts
concept to define a prior distribution for the parameters
governing the likeliness of the occurrences of the pitches.
In the next subsection, we describe specifically how to design the prior distributions.

+ ln p(R) +∑∑ln p(Ωk )kln p( Ak ).

(31)

k

=c denotes equality up to constant terms. Since the firstterm of Eq. (30) involves summation over k and n, analytically solving the current maximization problem is intractable. However, we can develop a computationally eﬃcient algorithm for finding a locally optimal solution basedon the auxiliary function concept, by using a similar idea
described in [8, 12].
When applying an auxiliary function approach to a certain maximization problem, the first step is to define alower bound function for the objective function. As mentioned earlier, the diﬃculty with the current maximization
problem lies in the first term in Eq. (30) . By using the factthat the logarithm function is a concave function, we can
invoke the Jensen’s inequality

3.2 Designing prior distributions
The likeliness of the pitch occurrences in popular and classical western music usually depend on the key or the chordused in that piece. The likeliness of the pitch occurrences
can be described as a probability distribution over the relative energies of the sounds of the individual pitches.
Since the number of times each note is activated is usually limited, inducing sparsity to the temporal activation ofeach note event would facilitate the source separation. The
likeliness of the number of times each note is activated canbe described as well as a probability distribution over the
temporal activations of the sound of each pitch.
To allow for designing such prior distributions, we decompose Uk,m as the productthe pitch∑∑ of two variables:
wise relative energy Rk = m Uk,m (i.e. k Rk = 1), andthe pitch-wise normalized amplitude Ak,m = Uk,m /Rk (i.e.
∑m Ak,m = 1). Hence, we can write
Uk,m = Rk Ak,m .
(27)This decomposition allows us to incorporate diﬀerent
kinds of prior information into our model by separately
defining prior distributions over R = (R1 , . . . , RK )T and

Yl,m ln Xl,m ≥ Yl,m∑

λk,n,l,m ln

w2k,n,m e−k,n(xl −Ωk,m −log n)22σ2

λk,n,l,m

Uk,m

,

(32)
to obtain a lower bound function,∑ where λk,n,l,m is a positive
variable that sums to unity: k,n λk,n,l,m = 1. Equality of(32) holds if and only if
λk,n,l,m =w2k,n,m e−(xl −Ωk,m −log n)2
2σ2

Xl,m

Uk,m

.

(33)

Although one may notice that the second term in
Eq. (30) is nonlinear in Ωk,m , the summation of Xl,mover
fairly well using the integral∫ ∞ l can be approximated∑Xisthe sum of the values at the
X(x,
t
)dx,
since
l,m
m
l
−∞
sampled points X(x1 , tm ), . . . , X(xL , tm ) with an equal interval, say ∆ x . Hence,∫ ∞∑1
Xl,m ≃X(x, tm )dx∆ x −∞l∫ ∞ (x−Ω −log n)2
k,m
1 ∑ 2
=
wk,n,m Uk,me− 2σ2dx∆ x k,n−∞626

15th International Society for Music Information Retrieval Conference (ISMIR 2014)
√

∑
2πσ ∑
Uk,m
w2k,n,m .
∆x
n
k

(34)

5555

Frequency [Hz]=This approximation implies that the second term in Eq.
(30) depends little on Ωk,m .
An auxiliary function can thus be written as+

L (Θ, λ) =

∑

Yl,m

∑

λk,n,l,m ln

w2k,n,m e−−ln n)2(xl −Ωk,m2σ2

λk,n,l,m
l,m
k,n
√
∑
2πσ ∑ ∑−
Uk,mw2k,n,m + ln p(Θ).
∆x
n
m
k
c

Uk,m

1750
551
174

2

4

6

8

Time [s](35)

Figure 2. Power spectrogram of a mixed audio signal ofthree violin vibrato sounds (D♭4, F4 and A♭4).

We can derive update equations for the model parameters,
using the above auxiliary function. By setting at zero thepartial derivative of L+ (Θ, λ) with respect to each of themodel parameters, we obtain∑l Yl,m λk,n,l,m + 1/2, (36)
w2k,n,m ← √Ω
2πRk Ak,m σ/∆ x + ν2 /(2|Bk (e jne k,m u0 )|2 )−1∑ αlαgΩk ←  2 D + 2 I M +diag(pk,n,l )
τ
υk
n,l


∑
 αg(37)× µk 2 1 M +(xl − ln n)pk,n,l  ,
υk
n,l
∑
∑
(R)
l,m Yl,m n λk,n,l,m + γk − 1
Rk ∝
,
(38)
∑
2
m,n Ak,m wk,m,n∑∑(A)
l Yl,m n λk,n,l,m + γk,m − 1
Ak,m ∝,(39)∑
Rk n w2k,m,n]⊤
1 [
pk,n,l := 2 Yl,1 λk,n,l,1 , Yl,2 λk,n,l,2 , · · · , Yl,M λk,n,l,M , (40)σwere artificially made by mixing D♭4, F4 and A♭4 violin
vibrato sounds from the RWC instrument database [3]. In
this paper, the F0 of the pitch name A4 was set at 440
Hz. The power spectrogram of the mixed signal is shownin Fig. 2. To convert the signal into a spectrogram, we
employed the fast approximate continuous wavelet transform [9] with a 16 ms time-shift interval. {xl }l ranged 55
to 7040 Hz per 10 cent. The parameters of HTFD wereset at γ(A)= (1 − 3.96 × 10−6 )1I , (τk , vk ) = (0.83, 1.25)
kfor all k, (N, K, σ, αg , α s ) = (8, 73, 0.02, 1, 1), and γ(R) =(1−2.4×10−3 )1K . {µk }k ranged A1 to A♯7 with a chromaticinterval, i.e. µk = ln(55) + ln(2) × (k − 1)/12. The number
of NMF bases were set at three. The parameter updates ofboth HTFD and NMF were stopped at 100 iterations.
While the estimates of spectrograms obtained with
NMF were flat and the vibrato spectra seemed to be averaged (Fig. 3 (a)), those obtained with HTFD tracked the F0
contours of the vibrato sounds appropriately (Fig. 3 (b)),
and clear vibrato sounds were contained in the separated
audio signals by HTFD.

where diag( p) converts a vector p into a diagonal matrix
with the elements of p on the main diagonal.
As for the update equations for the AR coeﬃcients β,
we can invoke the method described in [23] with a slight
modification, since the terms in the auxiliary function that
depend on β has the similar form as the objective function
defined in [23]. It can be shown that L+ can be increasedby the following updates (the details are omitted owing tospace limitations):hk ← Ĉk (βk )βk , βk ← Ck−1 hk ,(41)

5.2 Separation using key information
We next examined whether the prior information of asound improve source separation accuracy. The key of thesound used in 5.1, was assumed as D♭ major. The key information was incorporated in the estimation scheme bysetting γk(R) = 1 − 2.4 × 10−3 for the pitch indices that are
not contained in the D♭ major scale and γk(R) = 1−3.0×10−3
for the pitch indices contained in that scale. The other conditions were the same as 5.1.
With HTFD without using the key information, the estimated activations of the pitch indices that were not contained in the scale, in particular D4, were high as illustrated in Fig. 4 (a). In contrast, those estimated activations
with HTFD using the key information were suppressed asshown in Fig. 4 (b). These results thus support strongly thatincorporating prior information improve the source separation accuracy.

where Ck and Ĉk (βk ) are (P + 1) × (P + 1) Toeplitz matrices,
whose (p, q)-th elements are2
1 ∑ wk,m,n
Ck,p,q =cos[(p − q)neΩk,m u0 ],MN m,n 2ν
1 ∑1
Ĉk,p,q (βk ) =cos[(p − q)neΩk,m u0 ].
jne
MN m,n |Bk (e Ωk,m u0 )|2(42)
5. EXPERIMENTS

5.3 Transposing from one key to another
Here we show some results of an experiment on automatic
key transposition [11] using HTFD. The aim of key transposition is to change the key of a musical piece to anotherkey. We separated the spectrogram of a polyphonic soundinto spectrograms of individual pitches using HFTD, transposed the pitches of the subset of the separated components, added all the spectrograms together to construct apitch-modified polyphonic spectrogram, and constructed a

In the following preliminary experiments, we simplified
HTFD by omitting the source filter model and assuming
the time-invariance of wk,m,n .
5.1 F0 tracking of violin soundTo confirm whether HTFD can track the F0 contour ofa sound, we compared HTFD with NMF with the Idivergence, by using a 16 kHz-sampled audio signal which627

15th International Society for Music Information Retrieval Conference (ISMIR 2014)
A♭4

311
0

2

4

Time [s]

6

554

311
0

2

4

Time [s]

6

F4
Pitch

554

988

Frequency [Hz]988

Frequency [Hz]Frequency [Hz]

988

554

311
Time

0

2

4

Time [s]6(a) Without key information(a) Estimates of spectrograms and F0 contours (orange lines) obtained with HTFD

311
0

2

4

Time [s]

6

554

311
0

2

4

Time [s]

6

F4
Pitch

554

A♭4

988

Frequency [Hz]988

Frequency [Hz]Frequency [Hz]

988

D4
D♭4

554

D4
D♭4

311
Time

0

2

4

Time [s]6(b) With key information(b) Estimates of spectrograms obtained with NMF

Figure 4. Temporal activations of
Figure 3. Estimated spectrogram models by harmonic-temporal factor decomposi- A3–A♭4 estimated with HTFD using
tion (HTFD) and non-negative matrix factorization (NMF). In left-to-right fashion, and without using prior informationof the key. The red curves represent
the spectrogram models are for D♭4, F4 and A♭4.
the temporal activations of D4.
time-domain signal from the modified spectrogram using
the method described in [13]. For the key transposition,
we adopted a simple way: To transpose, for example, from
A major scale to A natural minor scale, we changed thepitches of the separated spectrograms corresponding to C♯,F♯ and G♯ to C, F and G, respectively.
Some results are demonstrated in http://hil.t.
u-tokyo.ac.jp/˜nakamura/demo/HTFD.html.
6. CONCLUSION
This paper proposed a new approach for monaural source
separation called the “Harmonic-Temporal Factor Decomposition (HTFD)” by introducing a spectrogram model thatcombines the features of the models employed in the NMF
and HTC approaches. We further described some ideas
how to design the prior distributions for the present model
to incorporate musically relevant information into the separation scheme.
7. ACKNOWLEDGEMENTS
This work was supported by JSPS Grant-in-Aid for Young
Scientists B Grant Number 26730100.
8. REFERENCES[1] A. S. Bregman: Auditory Scene Analysis, MIT Press, Cambridge, 1990.
[2] J. S. Downie, D. Byrd, and T. Crawford: “Ten years of ISMIR: Reflections on challenges and opportunities,” Proc. ISMIR, pp. 13–18, 2009.
[3] M. Goto: “Development of the RWC Music Database,” Proc.
ICA, pp. l–553–556, 2004.
[4] R. Hennequin, R. Badeau, and B. David: “Time-dependent
parametric and harmonic templates in non-negative matrix
factorization,” Proc. DAFx, pp. 246–253, 2010.
[5] R. Hennequin, B. David, and R. Badeau: “Score informed
audio source separation using a parametric model of nonnegative spectrogram,” Proc. ICASSP, pp. 45–48, 2011.
[6] G. E. Hinton: “Training products of experts by minimizing contrastive divergence,” Neural Comput., vol. 14, no. 8,
pp. 1771–1800, 2002.
[7] G. Hu, and D. L. Wang: “An auditory scene analysis approach
to monaural speech segregation,” Topics in Acoust. Echo andNoise Contr., pp. 485–515, 2006.
[8] H. Kameoka: Statistical Approach to Multipitch Analysis,
PhD thesis, The University of Tokyo, Mar. 2007.

[9] H. Kameoka, T. Tabaru, T. Nishimoto, and S. Sagayama:(Patent) Signal processing method and unit, in Japanese, Nov.
2008.
[10] H. Kameoka, T. Nishimoto, and S. Sagayama: “A multipitch
analyzer based on harmonic temporal structured clustering,”
IEEE Trans. ASLP, vol. 15, no. 3, pp. 982–994, 2007.
[11] H. Kameoka, J. Le Roux, Y. Ohishi, and K. Kashino: “Music
Factorizer: A note-by-note editing interface for music waveforms,” IPSJ SIG Tech. Rep., 2009-MUS-81-9, in Japanese,
Jul. 2009.
[12] H. Kameoka: “Statistical speech spectrum model incorporating all-pole vocal tract model and F0 contour generating process model,” IEICE Tech. Rep., vol. 110, no. 297, SP2010-74,
pp. 29–34, in Japanese, Nov. 2010.
[13] T. Nakamura and H. Kameoka: “Fast signal reconstruction
from magnitude spectrogram of continuous wavelet transform based on spectrogram consistency,” Proc. DAFx, 40, toappear, 2014.
[14] M. Nakano, J. Le Roux, H. Kameoka, Y. Kitano, N. Ono,and S. Sagayama: “Nonnegative matrix factorization with
Markov-chained bases for modeling time-varying patterns inmusic spectrograms,” Proc. LVA/ICA, pp. 149–156, 2010.
[15] A. Ozerov, C. Févotte, R. Blouet, and J. L. Durrieu: “Multichannel nonnegative tensor factorization with structured
constraints for user-guided audio source separation,” Proc.
ICASSP., pp. 257–260, 2011.
[16] S. A. Raczyński, N. Ono, and S. Sagayama: “Multipitch analysis with harmonic nonnegative matrix approximation,” Proc.
ISMIR, pp. 381–386, 2007.
[17] D. Sakaue, T. Otsuka, K. Itoyama, and H. G. Okuno:“Bayesian nonnegative harmonic-temporal factorization andits application to multipitch analysis,” Proc. ISMIR, pp. 91–
96, 2012.
[18] U. Simsekli and A. T. Cemgil: “Score guided musical source
separation using generalized coupled tensor factorization,”
Proc. EUSIPCO, pp. 2639–2643, 2012.
[19] P. Smaragdis and J. C. Brown: “Non-negative matrix factorization for polyphonic music transcription,” Proc. WASPAA,
pp. 177–180, 2003.
[20] P. Smaragdis and G. J. Mysore: “Separation by ”humming”:
User-guided sound extraction from monophonic mixtures,”Proc. WASPAA, pp. 69–72, 2009.
[21] E. Vincent, N. Bertin, and R. Badeau: “Harmonic and inharmonic nonnegative matrix factorization for polyphonic pitchtranscription,” Proc. ICASSP, pp. 109–112, 2008.
[22] K. Yoshii and M. Goto: “Infinite latent harmonic allocation:
A nonparametric Bayesian approach to multipitch analysis,”
Proc. ISMIR, pp. 309–314, 2010.
[23] A. El-Jaroudi, J. Makhoul: “Discrete all-pole modeling,”
IEEE Trans. SP, vol. 39, no. 2, pp. 411–423, 1991.

628