SHIFTED AND CONVOLUTIVE SOURCE-FILTER NON-NEGATIVE MATRIXFACTORIZATION FOR MONAURAL AUDIO SOURCE SEPARATION
Tomohiko Nakamura† and Hirokazu Kameoka†,‡†

Graduate School of Information Science and Technology, The University of Tokyo.
7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan.
‡
NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation.
3-1, Morinosato Wakamiya, Atsugi, Kanagawa, 243-0198, Japan

This paper proposes an extension of non-negative matrix factorization (NMF), which combines the shifted NMF model with thesource-filter model. Shifted NMF was proposed as a powerful approach for monaural source separation and multiple fundamental frequency (F0 ) estimation, which is particularly unique in that it takes
account of the constant inter-harmonic spacings of a harmonic structure in log-frequency representations and uses a shifted copy of aspectrum template to represent the spectra of diﬀerent F0 s. However,for those sounds that follow the source-filter model, this assumption
does not hold in reality, since the filter spectra are usually invariant
under F0 changes. A more reasonable way to represent the spectrum of a diﬀerent F0 is to use a shifted copy of a harmonic structure template as the excitation spectrum and keep the filter spectrum
fixed. Thus, we can describe the spectrogram of a mixture signal
as the sum of the products between the shifted copies of excitation
spectrum templates and filter spectrum templates. Furthermore, thetime course of filter spectra represents the dynamics of the timbre,
which is important for characterizing the feature of an instrument
sound. Thus, we further incorporate the non-negative matrix factor deconvolution (NMFD) model into the above model to describe
the filter spectrogram. We derive a computationally eﬃcient andconvergence-guaranteed algorithm for estimating the unknown parameters of the constructed model based on the auxiliary function
approach. Experimental results revealed that the proposed method
outperformed shifted NMF in terms of the source separation accuracy.
Index Terms— Audio source separation, Shifted non-negative
matrix factorization, Shift-invariant probabilistic latent component
analysis, Source-filter theory
1. INTRODUCTION
One major approach to monaural source separation involves applying non-negative matrix factorization (NMF) to an observed magnitude (or power) spectrogram interpreted as a non-negative matrix [1].
While many variants and extensions of NMF have been applied tospectrograms with linear frequency resolution given by the shorttime Fourier transform (STFT), spectrograms with log-frequency
resolution obtained with the continuous wavelet transform (CWT)
This study was supported by JSPS KAKENHI Grant Number 15J09992
and JSPS Grant-in-Aid for Young Scientists B Grant Number 26730100.

978-1-4799-9988-0/16/$31.00 ©2016 IEEE

489

Magnitude

ABSTRACT

48

85

151

269

480

855 1523 2714 4836

Log−Frequency [Hz]Fig. 1: Two spectra of clarinet sounds at diﬀerent pitches.

are particularly well suited to music data in the sense that the fundamental frequencies (F0 s) of the tones in music are geometricallyspaced [2–7]. In particular, shifted NMF [2], also known as theshift-invariant probabilistic latent component analysis (PLCA) [3],has proved successful for monaural source separation and multiple
F0 estimation tasks [8–10]. The uniqueness of this method lies inthat it takes account of the constant inter-harmonic spacings of aharmonic structure in log-frequency representations and uses shifted
copies of a spectrum template to represent the spectra of diﬀerent
F0 s where the shape of each spectrum template is assumed to be associated with an individual instrument. It should be noted that thisidea was first introduced in [11, 12] to develop a method for multiple F0 estimation. However, the above assumption does not hold in
real situations. Fig. 1 shows two examples of the spectra producedby the same instrument with diﬀerent F0 s. As Fig. 1 shows, the relative energies of the harmonic partials of the two spectra appear to
be completely diﬀerent, so that it is diﬃcult to represent them onlywith a single template. To address this mismatch, previous work proposed using multiple templates to represent the spectra of the same
instrument [2, 8–10].
To construct a more reasonable model, we focus on the fact thatthe generating processes of many instrument sounds can be explained fairly well by the source-filter theory. Specifically, we assume that the spectrum of an instrument sound is given by the product of the excitation and filter spectra, where the “excitation” source
represents a vibrating object such as a violin string and “filter” refers
to the resonance structure of the instrument body, which is usuallyinvariant under F0 changes. While most of the NMF variants incorporating the source-filter model such as [13–17] were designed
to model STFT spectrograms, we consider modeling CWT spectrograms so as to utilize the shift-invariant property of the excitation

ICASSP 2016

spectrum. We can thus describe the spectrum of each instrument at
a particular time as the product of a shifted copy of an excitation
spectrum template and a filter spectrum.
As regards the filter spectrum, if it can be assumed to be constantover time, we can use a single spectrum template for each instrument. However, for some instruments such as a singing voice, thetime course of the filter spectrum represents the dynamics of thetimbre, which is important for characterizing the feature of an instrument. We thus consider it reasonable to represent the filter spectrogram using a short-range spectrogram rather than a single-frame
spectrum as the template. This can be achieved by introducing theconvolutive NMF model [18] to express the filter spectrogram ofeach instrument.
In summary, this paper proposes modeling the CWT spectrogramof a polyphonic music as the sum of components each described by
the product of the excitation and filter spectrogram, where the excitation and filter spectrograms are respectively expressed using theshifted NMF model and the convolutive NMF model. We further derive a convergence-guaranteed iterative algorithm for minimizing thediﬀerence between an observed spectrogram and the present model
based on the auxiliary function approach [19, 20].
2. SHIFTED AND CONVOLUTIVE SOURCE-FILTERNON-NEGATIVE MATRIX FACTORIZATION
2.1. Spectrogram Model
Let us denote the indices of log-frequency and time by l = 0, . . . , L −1 and m = 0, . . . , M − 1, respectively. A spectrogram of an audio
signal that follows the source-filter model can be described as the
product of an excitation spectrogram and a filter spectrogram. By using the fact that the inter-harmonic spacings of a harmonic structure
are constant in the log-frequency domain, we use a shifted copy ofan excitation spectrum template to represent the excitation spectrum(ex)of each F0 . The excitation spectrogram X̃l,m,k≥ 0 of source excitation k(= 0, · · · , K − 1) is modeled as the convolution of an excitation(ex)
spectrum template S k,l ≥ 0 and time-varying gains Uk,p,m≥ 0, i.e.
∑
(ex)
(ex)
, where p is the frequency shift index andX̃l,m,k= p∈P S k,l−p Uk,p,m
P is the set of possible frequency shifts. By abuse of notation, we
understand that S k,l−p = 0 unless 0 ≤ l − p ≤ L − 1.
On the other hand, we describe the filter spectrogram in a similar manner to NMFD [18] and NMF-2D [4] to capture the dynamics of the timbre, which is important for characterizing the fea(filt)
ture of an instrument sound. The spectrogram X̃l,m,r≥ 0 of filterr(= 0, · · · , R − 1) is represented by a time convolution of a short(filt)≥ 0,
range spectrogram Fr,l,τ ≥ 0 with time-varying gains Ur,m
∑ M(tap) −1(filt)(filt)(tap)− 1 is thei.e. X̃l,m,r = τ=0Fr,l,τ Ur,m−τ , where τ = 0, · · · , M
time shift index and M (tap) is the tap size of the short-range spectro(filt)
grams. By abuse of notation, we understand that Ur,m−τ= 0 unless0 ≤ m − τ ≤ M − 1.
As we want the magnitude spectra of filters to be smooth andnon-negative in the log-frequency domain, we parameterize Fr,l,τ by
N envelope∑ kernels Gl,n ≥ 0 and their mixture weights Wr,n,τ ≥ 0 thatsatisfies n,τ Wr,n,τ = 1:∑
Fr,l,τ =Wr,n,τ Gl,n ,(1)n(ω −ρ )21
− l n
Gl,n := √
(2)
e 2ν2
2πν2
where n = 0, · · · , N − 1 is the index of envelope kernel and ωl ∈(0, π] is the normalized angular frequency corresponding to the lth490

log-frequency. The kernel Gl,n for n ≥ 1 is identical to a normal
distribution of a normalized angular frequency with mean ρn andvariance ν2 .
Multiple excitation and filter spectra can be used for an instrument to describe complex spectral changes, but we hereafter assign an excitation spectrum and a filter spectrum to each instrument(ex)(filt)
for the simplicity. By putting Uk,r,p,m−τ = Uk,p,m
Uk,r,p,m−τ
and treating Uk,r,p,m−τ itself as a parameter, the spectrogram of an instrument
sound associated with source excitation k and filter r can be written
as∑Fr,l,τ S k,l−p Uk,r,p,m−τ .
(3)X̃l,m,k,r =p,τAssuming the additivity of magnitude spectrograms as with conventional NMFs, the observed spectrogram can be represented as∑
Xl,m =
X̃l,m,k,r .
(4)
k,r

∑
To avoid the indeterminacy in scaling, we put l S k,l = 1 for all k.
Although a model similar to the above has been mentioned in[21], the temporal dynamics was not incorporated into the sourcefilter model in the literature. Any experimental evaluation was not
given in the literature and the incorporation of the source-filter model
into shifted NMF has yet been validated. We will thus confirm theeﬃcacy of the incorporation of the source-filter model in Sec. 4.
2.2. Formulation
For a given magnitude spectrogram Y := {Yl,m }l,m , we would liketo find the parameters S := {S k,p }k,p , W := {Wr,n,τ }r,n,τ and U :={Uk,r,l,m }k,r,l,m of the proposed model such that minimizes∑
L∗ (S , W, U) =D∗ (Yl,m ||Xl,m ) + R∗ (U).
(5)
l,mThe first term of Eq. (5) is a goodness-of-fit measure between Y andX := {Xl,m }l,m . How to define the measure is very important since
it corresponds to an assumption to the statistical nature of observed
data. If we define DI as the generalized Kullback-Leibler divergence(a.k.a I divergence), it implicitly assumes that Yl,m follows a Poisson
distribution with mean Xl,m :Yl,m
− Yl,m + Xl,m .
(6)
Xl,m
∑
From this fact, it is known that minimizing l,m DI (Yl,m ||Xl,m ) with respect to Xl,m amounts to the maximum likelihood estimation of Xl,m .
This measure is frequently used in conventional NMF algorithms
and has been confirmed to work well for audio source separation
empirically. Another commonly-used measure is the Itakura-Saito
divergence DIS :DI (Yl,m ||Xl,m ) = Yl,m ln

2||Xl,m ) =DIS (Yl,m

2
Yl,m

Xl,m

− ln

2
Yl,m

Xl,m

− 1.

(7)

This corresponds to the assumption that an observed complex spectrogram follows a circularly-symmetric complex normal distribution
with mean zero and variance Xl,m , in which Xl,m can be interpreted asa model of a power spectral density of the observed signal.
The second term R∗ (U) is a regularizer for U. In popular andclassical western music, the number of pitches occurred in a musicalpiece and the number of times each note is performed are usuallylimited, and so inducing the sparsity of U would facilitate the source
separation. To reflect it, we can design the regularizer in analogy to
the Bayesian modeling. The conjugate prior of the Poisson distribution is a gamma distribution Gam(x; α, β) ∝ xα−1 e−βx and thus wedesign the regularizer RI (U) for DI as∑ {}RI (U) =−(α(I) − 1) ln Uk,r,p,m + β(I) Uk,r,p,m ,(8)k,r,p,mwhere α(I) > 0 and β(I) > 0 are associated with the shape and rate
parameters of a gamma distribution, respectively. Similarly, the conjugate prior of a circularly-symmetric complex normal distribution
with known mean and unknown variance is an inverse gamma distribution InvGam(x; α, β) ∝ x−α−1 e−β/x , and we design the regularizer
for DIS as
}
∑ {
β(IS)
RIS (U) =(α(IS) + 1) ln Uk,r,p,m +,(9)
Uk,r,p,mk,r,p,mwhere α(IS) > 0 and β(IS) > 0 are associated with the shape and scaleparameters of an inverse gamma distribution, respectively. The lessα(I) (α(IS) ), the more sparse U tends to become.
3. PARAMETER ESTIMATION ALGORITHMS BASED ON
AUXILIARY FUNCTION APPROACH
We first derive a parameter estimation algorithm for the I divergence. Since LI (S , W, U) involves summations over k, r, p, τ andn in the logarithmic function, the current minimization problem isdiﬃcult to solve analytically. However, we can develop a computationally eﬃcient algorithm for finding a locally optimal solution based on the auxiliary function approach [19, 20, 22]. The firststep to apply the auxiliary function approach, is to define an upper bound function for the objective function L(S , W, U), arranged
as L+ (S , W, U, Λ), such that L(S , W, U) = minΛ L+ (S , W, U, Λ).
We call Λ an auxiliary variable and L+ (S , W, U, Λ) an auxiliary
function. If we can construct L+ (S , W, U, Λ), L(S , W, U) is nonincreasing under the updates {S , W, U} ← argmin L+ (S , W, U, Λ) and
S ,W,UΛ ← argmin L+ (S , W, U, Λ).
Λ

Since the logarithmic function is a concave function, we can obtain an upper bound function by invoking the Jensen’s inequality:∑(−Yl,m ln Xl,m ≤ − Yl,m
λl,m,k,r,p,τ,n ln S k,l−p + ln Wr,n,τ
k,r,p,τ,n+ ln Gl,n + ln Uk,r,p,m−τ − ln λl,m,k,r,p,τ,n

)

(10)

where
λl,m,k,r,p,τ,n ≥ 0 is an auxiliary variable such that∑k,r,p,τ,n λl,m,k,r,p,τ,n = 1 for all l and m. The equality holds if andonly if
S k,l−p Wr,n,τ Gl,n Uk,r,p,m−τ
.
λl,m,k,r,p,τ,n =
Xl,m

(11)

The auxiliary function can thus be written as∑∑(
L+I (S , W, U, Λ) = −Yl,mλl,m,k,r,p,τ,n ln S k,l−p + ln Wr,n,τ
c

l,m

∑ {

c

{λl,m,k,r,p,τ,n }l,m,k,r,p,τ,n .
By setting the partial derivatives of
L+ (S , W, U, Λ) with respect to S , W and U at zeros and substituting Eq. (11) into Λ, we can derive the following update equations:
Yl,m ∑
Fr,l,τ Uk,r,l−l′ ,m−τXl,m r,τ
∑,S k,l′ ←S k,l′l,m,r,τ F r,l,τ U k,r,l−l′ ,m−τ
∑ Yl,m ∑Gl,n S k,l−p Uk,r,p,m−τl,mXl,m p,k∑
Wr,n,τ ←Wr,n,τ
,l,m,p,k G l,n S k,l−p U k,r,p,m−τ∑ Yl,m ∑Fr,l,m−m′ S k,l−p Uk,r,p,m′ + α(I) − 1l,m
Xl,m p,k∑
Uk,r,p,m′ ←,(I)l,m,p,k F r,l,m−m′ S k,l−p + β
∑

l,m

(13)

(14)

(15)

The
W are followed by normalization such that
∑ update rules of S and∑l S k,l = 1 for all k andn,τ Wr,n,τ = 1 for all r. It is important tonote that once the initial values of W and S are set to be non-negative,
the multiplicative update equations ensures the non-negativity of theentries of W and S . Since the non-negativity of U does not hold,we can ensure it by simply performing Uk,r,p,m ← max{0, Uk,r,p,m }at each update. One may think that the update equations containtime-consuming convolutions and correlations and would require along computation time. However, we can invoke the fast Fouriertransform (FFT) to calculate the convolutions and correlations, andthey are not time-consuming in practice.
Similarly to the above, we can construct an auxiliary function for
the IS divergence as with [23]. We here omit details of the derivation of an auxiliary function due to limitations of space. Update
equations can be derived as
v
u
u
2
u
u
Yl,m
u
∑
u
u
u
Fr,l,τ Uk,r,p,ml,m,r,τ
u
2
u
u
Xl,m
t
′
′
,
(16)
S k,l ←S k,l
Fr,l,τ Uk,r,p,m∑l,m,r,τ
Xl,m
v
u
u
2
u
u
Y
u
∑
l,m
u
u
u
S k,l−p Gl,n Uk,r,p,ml,m,k,p
u
2
u
u
Xl,m
t
,
(17)
Wr,n,τ ←Wr,n,τ
S k,l−p Gl,n Uk,r,p,m∑l,m,k,pXl,mUk,r,p,m′ = √Ak,r,p,m′
(

α

(IS)

+1

(18)

)2

2

+B

k,r,p,m′A

k,r,p,m′α(IS) + 1
+
2

where
2
∑ Yl,m

Fr,l,m−m′ S k,l−p Uk,r,p,m′2 + β(IS) .
2
X
l,m
l,m
∑ S k,l−p Fr,l,m−m′
Bk,r,p,m′ =
.
Xl,m
l,m

Ak,r,p,m′ =

(19)
(20)

k,r,p,τ,n) ∑+ ln Uk,r,p,m−τ − ln λl,m,k,r,p,τ,n +Xl,m−where = denotes the equality up to constant terms and Λ :=l,m(α(I) − 1) ln Uk,r,p,m − β(I) Uk,r,p,m4. EXPERIMENTS}(12)k,r,p,m491

To evaluate the proposed algorithms, we conducted a supervised
source separation experiment. For the convenience, we call the proposed algorithm with the I divergence criterion (the IS divergence

Violin
Saxophone

Algorithm
I-SNMFwSF(0.6, 1.0 × 10−10 )IS-SNMFwSF(1.0, 0.6)
I-SNMF(1.0, 0.4)
IS-SNMF(0.4, 1.0)

SDR
SIR
improvements improvements

SDR improvement [dB]

Table 1: Average SDR improvements, SIR improvements andSARs with standard errors obtained with the proposed algorithms (ISNMFwSF and IS-SNMFwSF) and the shifted NMFs (I-SNMF and
IS-SNMF).
SARs

6.01 ± 0.58

11.06 ± 1.09 3.91 ± 0.64

4.77 ± 0.29

8.98 ± 0.40

3.46 ± 0.50

3.96 ± 0.44

9.58 ± 0.69

1.79 ± 0.85

2.81 ± 0.51

6.00 ± 0.61

4.09 ± 0.94

11
10
9
8
7
6
5
4
3
2
1
0
SNMF

criterion) I-SNMFwSF (IS-SNMFwSF, respectively). For comparison, we employed shifted NMF with the I divergence criterion (ISNMF) and that with the IS divergence criterion (IS-SNMF). While
the original shifted NMF [2] does not contain any terms inducing thesparsity of parameters, we found that R∗ (U) improved SDRs in theexperiment and we here used R∗ (U).
The experimental data was the Bach10 dataset [24], which consists of audio recordings of ten four-part chorales by J. S. Bach. Each
recording is a mixture of violin, clarinet saxophone and bassoon performances, which correspond to the soprano, alto, tenor and bass
parts of each musical piece, respectively. Audio recordings of individual parts are also contained in the dataset. All recordings weremonaural and downsampled to 16 kHz. Magnitude spectrograms
were computed with the fast approximate CWT algorithm [25, 26].
The center frequencies ranged from 27.5 to 7902 Hz with 100/3 centinterval and the log-normal wavelet [27] was used as an analyzing
wavelet. The wavelet has a Gaussian shape with a common variancein the log-frequency domain, and we set a parameter corresponding
to the standard deviation of the Gaussian as a one fifth of a semitone
interval.
We first trained S and W of the proposed models and basis spectra of the shifted NMFs with the audio recordings of individual
parts of the five musical pieces (training data), and then performed
source separation on the audio recordings of the other five musical pieces (test data). With the proposed algorithms, a pair of asource excitation and a filter was trained for each instrument, anda total of four pairs of a source and a filter were used for the separation. With the shifted NMFs, one basis spectrum was assigned toeach instrument and a total of four basis spectra were used for the
separation. For each test data, we designed a soft time-frequency
mask as X̃l,m,k,r /Xl,m to obtain separated audio signals of the sources.
The proposed methods and the shifted NMFs ran for 100 iterations both in the training and test stages. As α(I) or α(IS) , we useα(train) = 1.0 × 10−10 , 0.2, 0.4, 0.6, 0.8, 1.0 for the training data andα(test) = 1.0 × 10−10 , 0.2, 0.4, 0.6, 0.8, 1.0 for the test data. The otherparameters were set as follows: β(I) = β(IS) = 1.0 × 10−10 , M (tap) = 1,
N = 140, ν = π/(2N − 2) and ρn = πn/(N − 1) for n = 0, · · · , N − 1.
Table 1 shows SDR improvements, SIR improvements and SARs
obtained with all algorithms for all data. SDRs, SIRs and SARs werecomputed with the BSSEval toolbox [28]. The pairs of two values
below the algorithm names are (α(train) , α(test) ), which provided thehighest average SDR improvement for each algorithm. These results
show that the incorporation of the source-filter model improves the492

Clarinet
Bassoon

20 40 60 80 100 120 140 160 180 200

Fig. 2: Average SDR improvements and standard errors obtainedwith the proposed algorithm (I-SNMFwSF) and I-SNMF for eachmusical instrument. “SNMF” corresponds to I-SNMF.

source separation accuracy in the CWT domain.
We also conducted an experiment to examine the eﬀect of N of
the proposed algorithms. Fig. 2 displays average SDR improvements
and standard errors obtained with I-SNMFwSF and I-SNMF for individual musical instruments, where each algorithm ran with the same(α(train) , α(test) ) as in the above. I-SNMFwSF with 100 ≤ N ≤ 160
provided significantly higher SDR improvements for all musical instruments compared to I-SNMF. We found that the best N was different for each musical instrument, and so exploring the best Ns forother musical instruments and classifying them may be required forpractical use.
5. CONCLUSION
This paper has developed a new source separation method by incorporating the source-filter model into shifted NMF. With the proposed model, the observed spectrogram is represented by a productof excitation and filter spectrograms. The excitation spectrogram is
described with shifted NMF to exploit the constant inter-harmonic
spacings of a harmonic structure in the log-frequency domain, andthe filter spectrogram is modeled by NMFD to represent temporal
dynamics of timbre. We have derived iterative algorithms of estimating parameters for the objective functions using the I divergence and
IS divergence criterions based on the auxiliary function approach.
We have experimentally confirmed that the proposed algorithms outperformed shifted NMFs in the accuracy of source separation. In future, we will examine the eﬀect of setting M (tap) > 1 to the source
separation accuracy.

6. REFERENCES[1] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” in Proc. IEEE Workshop Applications Signal Process. Audio Acoust. IEEE, 2003,
pp. 177–180.
[2] D. FitzGerald, M. Cranitch, and E. Coyle, “Shifted nonnegative matrix factorisation for sound source separation,”in IEEE/SP 13th Workshop on Statistical Signal Processing.
IEEE, 2005, pp. 1132–1137.
[3] P. Smaragdis, B. Raj, and M. Shashanka, “Sparse and shiftinvariant feature extraction from non-negative data,” in Proc.
Int. Conf. Acoust. Speech Signal Process., Mar. 2008, pp.
2069–2072.
[4] M. N. Schmidt and M. Mørup, “Nonnegative matrix factor 2D deconvolution for blind single channel source separation,”in Proc. Int. Conf. Independent Component Analysis and Blind
Signal Separation, 2006, pp. 700–707.

[15] R. Hennequin, B. David, and R. Badeau, “Score informed audio source separation using a parametric model of non-negative
spectrogram,” in Proc. Int. Conf. Acoust. Speech Signal Process., 2011, pp. 45–48.
[16] J. J. Carabias-Orti, T. Virtanen, P. Vera-Candeas, N. RuizReyes, and F. J. Canadas-Quesada, “Musical instrument soundmulti-excitation model for non-negative spectrogram,” IEEE J.
Sel. Topics. Signal Process., vol. 5, no. 6, pp. 1144–1158, Oct.
2011.
[17] H. Kirchhoﬀ, S. Dixon, and A. Klapuri, “Missing templateestimation for user-assisted music transcription,” in Proc. Int.
Conf. Acoust. Speech Signal Process., 2013, pp. 26–30.
[18] P. Smaragdis, “Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs,”
in Proc. Int. Conf. Independent Component Analysis and Blind
Signal Separation, 2004, pp. 494–499.
[19] J. M. Ortega and W. C. Rheinboldt, Iterative solution of nonlinear equations in several variables, Number 30. SIAM, 1970.

[5] S. A. Raczyński, N. Ono, and S. Sagayama, “Multipitch analysis with harmonic nonnegative matrix approximation,” in Proc.
Int. Conf. Music Info. Retrieval, 2007, pp. 381–386.

[20] D. R. Hunter and K. Lange, “Quantile regression via an mm algorithm,” Journal of Computational and Graphical Statistics,
vol. 9, no. 1, pp. 60–77, 2000.

[6] E. Vincent, N. Bertin, and R. Badeau, “Harmonic and inharmonic nonnegative matrix factorization for polyphonic pitchtranscription,” in Proc. Int. Conf. Acoust. Speech Signal Process., Mar. 2008, pp. 109–112.

[21] T. Virtanen, “Unsupervised learning methods for source separation in monaural music signals,” in Signal Processing Methods for Music Transcription, A. Klapuri and M. Davy, Eds.,
2006, pp. 267–296.

[7] T. Nakamura, K. Shikata, N. Takamune, and H. Kameoka,“Harmonic-temporal factor decomposition incorporating music prior information for informed monaural source separation,” in Proc. Int. Symposium Music Info. Retrieval, 2014,
pp. 623–628.

[22] H. Kameoka, “Statistical speech spectrum model incorporating
all-pole vocal tract model and f0 contour generating processmodel,” in IEICE Technical Report, Nov. 2010, vol. 110, pp.
29–34, in Japanese.

[8] “2014: Multiple Fundamental Frequency Estimation & Tracking Results - MIREX Wiki,” http://www.music-ir.org/
mirex/wiki/2014:Multiple_Fundamental_Frequency_
Estimation_%26_Tracking_Results.
[9] B. Fuentes, R. Badeau, and G. Richard, “Harmonic adaptive
latent component analysis of audio and application to musictranscription,” IEEE Trans. Acoust., Speech, and Language
Process., vol. 21, no. 9, Sept. 2013.
[10] R. Jaiswal, D. FitzGerald, E. Coyle, and S. Rickard, “Towards
shifted NMF for improved monaural separation,” in Proc. 24th
IET Irish Signals and Systems Conference, July 2013.
[11] S. Sagayama, H. Kameoka, and T. Nishimoto, “Specmurt analysis: A piano-roll-visualization of polyphonic music signal bydeconvolution of log-frequency spectrum,” in Proc. ISCA Tutorial and Research Workshops on Statistical and Perceptual
Audition.
[12] S. Saito, H. Kameoka, K. Takahashi, T. Nishimoto, andS Sagayama, “Specmurt analysis of polyphonic music signals,”IEEE Trans. Acoust., Speech, and Language Process., vol. 16,
no. 3, pp. 639–650, Mar. 2008.
[13] A. Klapuri, “Analysis of musical instrument sounds by sourcefilter-decay model,” in Proc. Int. Conf. Acoust. Speech Signal
Process., 2007, vol. 1, pp. I–53–I–56.
[14] H. Kameoka and K. Kashino, “Composite autoregressive system for sparse source-filter representation of speech,” in IEEE
International Symposium on Circuits and Systems, 2009, pp.
2477–2480.

493[23] C. Févotte, N. Bertin, and J. L. Durrieu, “Nonnegative matrix
factorization with the itakura-saito divergence: With application to music analysis,” Neural computation, vol. 21, no. 3, pp.
793–830, 2009.
[24] Z. Duan and B. Pardo, “Soundprism: An online system forscore-informed source separation of music audio,” IEEE J. Sel.
Topics. Signal Process., vol. 5, no. 6, pp. 1205–1215, 2011.
[25] H. Kameoka, T. Tahara, T. Nishimoto, and S. Shigeki, “Signal processing method and device,” Nov. 2008, Japan Patent
JP2008-281898.
[26] T. Nakamura and H. Kameoka, “Fast signal reconstruction
from magnitude spectroggram of continuous wavelet transform
based on spectrogram consistency,” in Proc. Int. Conf. Digital
Audio Eﬀects, 2014, pp. 129–135.
[27] H. Kameoka, Statistical Approach to Multipitch Analysis,
Ph.D. thesis, The University of Tokyo, Mar. 2007.
[28] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans.
Acoust., Speech, and Language Process., vol. 14, no. 4, pp.
1462–1469, 2006.