2015 23rd European Signal Processing Conference (EUSIPCO)UNIFIED APPROACH FOR AUDIO SOURCE SEPARATION WITH MULTICHANNEL
FACTORIAL HMM AND DOA MIXTURE MODEL
Takuya Higuchi1) and Hirokazu Kameoka1),2)
1)

Graduate School of Information Science and Technology, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
2)
NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation
3-1 Morinosatowakamiya, Atsugi, Kanagawa 243-0198, Japan{higuchi,kameoka}@hil.t.u-tokyo.ac.jp
ABSTRACTWe deal with the problems of blind source separation, dereverberation, audio event detection and direction-of-arrival(DOA) estimation. We previously proposed a generativemodel of multichannel signals called the multichannel factorial hidden Markov model, which allows us to simultaneously
solve these problems through a joint optimization problem
formulation. In this approach, we modeled the spatial correlation matrix of each source as a weighted sum of thespatial correlation matrices corresponding to all possible
DOAs. However, it became clear through real environment
experiments that the estimate of the spatial correlation matrix
tended to deviate from the actual correlation matrix since
the plane wave assumption does not hold due to reverberation and noise components. To handle such deviations,
we propose introducing a prior distribution over the spatial
correlation matrices called the DOA mixture model insteadof using the weighted sum model. The experiment showedthat the proposed method provided 1.94 [dB] improvementcompared with our previous method in terms of the the signalto-distortion ratios of separated signals.
Index Terms— Blind source separation, voice activity
detection, dereverberation, DOA estimation
1. INTRODUCTION
Blind source separation (BSS) refers to a technique for separating out individual source signals from microphone array
inputs when the transfer characteristics between the sources
and microphones are unknown. To solve the BSS problem,
it is generally necessary to make some assumptions about thesources, and formulate an appropriate optimization problem
based on criteria designed according to those assumptions.
For example, if the observed signals outnumber the sources,
we can employ independent component analysis (ICA) [1]by assuming that the sources are statistically independent ofeach other. For monaural source separation, one successful
approach involves applying non-negative matrix factorization(NMF) to the magnitude spectrogram of a mixture signal, interpreted as a non-negative matrix [2, 3]. Up to now, several attempts have been made to extend this approach to a
multichannel case in order to allow us to use the spatial information as an additional clue to source separation [4, 5].
Moreover, we previously proposed the multichannel factorial
hidden Markov model [6], where we used the information ofaudio events for source separation and simultaneously solved
the problems of BSS and audio event detection. Furthermore,
This work was supported by JSPS KAKENHI Grant Number 26730100.

978-0-9928626-3-3/15/$31.00 ©2015 IEEE

2043

we extended the model for dereverberation [7], by approximating the multichannel observed signal recorded in a reverberant condition as a form of a convolution of the frequency
array response and the source signal in the time-frequency
domain. Thus, we modeled the impulse response out of theframe of STFT by introducing a time sequence of frequency
response arrays in the time-frequency domain. Furthermore,
we used direction-of-arrivals (DOAs) of sources as a clue toestimation of spatial correlation matrices [8]. Generally, it is
known that the spatial correlation matrix of the direct sound of
a point source has a certain structure described by the DOA.
In the model proposed in [8], we used a weighted sum of matrices parametrized by all possible DOAs, same as in [9], andestimates of the weight variables corresponded to estimatesof spatial correlation matrices at all frequencies. However,in practical situations, a spatial correlation matrix often deviates from a theoretical structure because of reverberation andnoises within the frame of the STFT. Therefore, in this paper, we apply the DOA mixture model [10] for modeling thegenerative process of a spatial correlation matrix, and design
a prior distribution over a spatial correlation matrix based on
DOAs. In the present model, we estimate a spatial correlation matrix at each frequency bin based on DOAs. Throughthe parameter inference of our new generative model, we cansimultaneously perform source separation, source activity detection, dereverberation and DOA estimation based on a uniﬁed maximum likelihood criterion.
2. MULTICHANNEL FACTORIAL HIDDENMARKOV MODEL2.1. Mixing model with convolutive mixture approximation [7]First we consider a situation where I source signals arerecorded by M microphones. In a reverberant condition, thelength of the impulse responses are relatively long and soan instantaneous mixture approximation is not always true.
Therefore we approximately express the observed signals as
a form of a convolution of the frequency array response andthe source signal in the time-frequency domain.
TI y(ωk , tl ) ≈ai (ωk , tτ )si (ωk , tl − tτ ).
(1)i=1 τ =0

Here, let ym (ωk , tl ) ∈ C be the short-time Fourier transform (STFT) component observed at the m-th microphone,and si (ωk , tl ) ∈ C be the STFT component of the i-th source
signal. 1 ≤ k ≤ K and 1 ≤ l ≤ L are the frequency

2015 23rd European Signal Processing Conference (EUSIPCO)

and time indices, respectively. ai (ωk , tτ ) denotes the frequency array response for source i at frequency ωk and timetτ . 0 ≤ τ ≤ T is the time index of the frequency array response in the time-frequency domain. Note that ai (ωk , t1:T )denote the frequency array responses which correspond to theimpulse responses out of the frames of the STFT. For convenience of notation, we hereafter use subscripts k, l and τ toindicate ωk , tl and tτ respectively.
2.2. Generative modeling of source signals based on audio
events [6]
We now describe the generative process of the source signal
si,k,l based on its audio event. First, we assume each signal
2, which ishas a speciﬁc spectrum and utilize NMF for σi,k,l2
an expected value of power of si,k,l . σi,k,l is factorized as2= wi,k hi,l ,(2)
σi,k,lwhere wi,k and hi,l are non-negative variables. The generative model of si,k,l is written as(3)
si,k,l |wi,k , hi,l ∼ NC (si,k,l ; 0, wi,k hi,l ),conditioned on wi,k and hi,l . Under the condition of Eq. (3),the generative model of the multichannel observed signal isequivalent to the method proposed in [4, 5].
In many cases, a sound signal consists of several spec(q)tra. Therefore, ﬁrst we introduce latent variable zi,l ∈{1, . . . , Q} to denote a state of i-th source’s spectrum at(q)(q)
time l, and the state sequence zi,1 , . . . , zi,L follows a Markov
chain:(q) (q)(q)(q)
zi,l |zi,l−1 ∼ Categorical(zi,l ; ρ (q) ),(4)i,zi,l−1(q)yx , ρi,q

(q)

(q)

= (ρi,q,1 , . . . , ρi,q,Q )where Categorical(x; y) =denotes the transition probability of state q to each state(q)(q)
1, . . . , Q, and ρi = (ρi,q,q )Q×Q denotes the transition matrix. Then the power spectrum of the i-th source at time l is(q)assumed to be determined according to zi,l and wi,k,q denotes a spectrum of i-th source at state q. Thus, the generativemodel of si,k,l is rewritten as(q)si,k,l |wi,k,1:Q , hi,l , zi,l ∼ NC (si,k,l ; 0, wi,k,z(q) hi,l ).
i,l

(5)

Next, the source’s power also tend to differ according to thesource’s state such as “active” or “inactive.” Therefore we(j)ﬁrst introduce latent variable zi,l ∈ 1, . . . , J to denote astate of i-th source’s power at time l, and the state sequence(j)(j)
zi,1 , . . . , zi,L follows a Markov chain:
(j)

(j)

(j)

zi,l |zi,l−1 ∼ Categorical(zi,l ; ρ

(j)

(j)

i,zi,l−1).

(6)Then we assume hi,l follows a gamma distribution which has(j)different parameters according to zi,l ,(j)hi,l |zi,l ∼ Gamma(hi,l ; αz(j) βz(j) ),
i,l
i,l

(7)

where α1:J and β1:J are the shape and scale parameters of
α−1 −x/βe.
a gamma distribution, and Gamma(x; α, β) = x Γ(α)βα(j)
As we want hi,l to take a small value when zi,l is the “in-active” (i.e., silent) state, we set the hyperparameters of the2044

gamma distribution of that state so that it becomes a sparsityinducing distribution. As regards the gamma distributions of
the remaining states, we consider setting the hyperparameters
so that they are like uniform distributions.
2.3. Generative process of a spatial correlation matrix
based on DOA mixture model
If we assume that a source is far from the microphones, thefrequency array response has a certain structure depending on
DOA. Therefore, we can express a spatial correlation matrix
by using the DOAs. Since a spatial correlation matrix C i,k,τ
is deﬁned as C i,k,τ = ai,k,τ aH
i,k,τ by using ai,k,τ , speciﬁcally, with M = 2 microphones, the spatial correlation matrix
for a source at direction θ such that 0 ≤ θ ≤ π is deﬁned as afunction of ω depending on θ



∗
1
J (θ, ω) = jωB cos θ/C 1 ejωB cos θ/C ,(8)ewhere j is the imaginary unit, B [m] is the distance between
the two microphones, and C [m/s] is the speed of sound. If the
DOA θi of source i is known, the spatial correlation matrix forthe direct wave should be equal to J (θi , ωk ). However, C i,k,0
would not be equal to J (θi , ωk ) because of reverberant components and noises within the frame of the STFT. One way
is modeling a spatial correlation matrix by using a weighted
sum of J (θi , ωk ) of all possible DOAs, same as in [8, 9]. In
that model, estimates of the weight variables corresponded
to estimates of spatial correlation matrices at all frequencies,
and so it would not be enough for modeling the deviation of
spatial correlation matrices. Therefore, in this paper, we incorporate the DOA mixture model [10] into the multichannel
factorial hidden Markov model, and estimate a spatial correlation matrix at each frequency bin using the information of
DOAs. First, we now introduce a discrete set of O possible directions, ϑ1 , . . . , ϑO , which are all assumed to be constants. We then assume that each source signal propagates(o)from one of these directions. For each source i, an index zi
of direction is drawn according to a categorical distribution(o)ρi = (ρi,1 , . . . , ρi,O )
(o)

(o)

(o)

(o)

zi |ρi ∼ Categorical(zi ; ρi ),θi = ϑz(o) .

(9)
(10)

i

Then, a spatial correlation matrix follows a whishart distribution:(o)C i,k,0 |zi ∼ WC (C i,k,0 ; γ, J ϑ (o) ,k + I),
(11)
z

i

exp(− 12 tr(XY −1 ))where WC (X; γ, Y ) ∝ |X|
and γ(≥ M + 1) is a hyperparameter of the whishart distribution. Here I is an identity matrix and is a small value,
and they enable inverse operation.
Our overall generative model is given by Eqs. (4), Eqs.
(6), Eqs. (7) Eqs. (9), Eqs. (11) and(γ−M )/2(q)

y k,l |a1:I,k,0:T , w1:I,k,1:Q , h1:I,l−T :l , z1:I,l−T :l∼ NC (y k,l ; 0,C i,k,τ wi,k,z(q) hi,l−τ ),
i,τ

i,l−τ

(q)

(12)

conditioned on a1:I,k,0:T , w1:I,k,1:Q , h1:I,l−T :l , z1:I,l−T :l .
Parameter estimation of the generative model allows us tosolve the problems of source separation, source activity detec-

2015 23rd European Signal Processing Conference (EUSIPCO)+

tion, dereverberation and DOA estimation based on a uniﬁed
maximum likelihood criterion.

(j)

λj,i,l log p(Z (j) ) + CL+ ,

j,i,l

(q)

(13)

(j)

where Ri,k,l,τ,q , U k,l , Ũ i,k , λq,i,l , λj,i,l and di,o are auxiliary
variables. Ri,k,l,τ,q , U k,l and Ũ i,k satisfy Hermitian positive(q)(j)
deﬁniteness and i,τ Ri,k,l,τ,q = I. λq,i,l , λj,i,l and di,o

3. PARAMETER ESTIMATION ALGORITHM
BASED ON AUXILIARY FUNCTION METHOD
In this section, we describe a parameter estimation algorithm for our generative model based on an auxiliary function
method. The random variables of interest in our model are
W = w1:I,1:K,1:Q , H = h1:I,1:L and C = C 1:I,1:K,0:T. We denote the entire set of the above parameters as Θ.
(q)
(j)
(o)
Z (q) = z1:I,1:L , Z (j) = z1:I,1:L and Z (o) = z1:I are latent variables. In the following, ρ(q) , ρ(j) , ρ(o) , α and β areconstants that are determined experimentally.
The objective function L(Θ) is deﬁned as L(Θ) =log p(Θ|Y ), where Y = y 1:K,1:L is a set consisting of
the time-frequency components of observed multichannel
signals. By using the conditional distributions deﬁned inSec. 2, we can write down log p(Θ|Y ). However, L(Θ) is
difﬁcult to maximize analytically, and so we derive an efﬁcient parameter estimation algorithm based on an auxiliary
function method.
The optimization problem of maximizing L(Θ) with respect to Θ is difﬁcult to solve analytically. However, we can
invoke an auxiliary function approach to derive an iterativealgorithm that searches for the estimate of Θ, as with [5]. To
apply an auxiliary function approach to the current optimization problem, the ﬁrst step is to construct an auxiliary function
L+ (Θ, Λ) satisfying L(Θ) = maxΛ L+ (Θ, Λ). We refer toΛ as an auxiliary variable. It can then be shown that L(Θ) is
non-decreasing under the updates Θ ← argmaxΘ L+ (Θ, Λ)and Λ ← argmaxΛ L+ (Θ, Λ). The proof of this shall beomitted owing to space limitations. Thus, L+ (Θ, Λ) should
be designed as a function that can be maximized analyticallywith respect to Θ and Λ. Such a function can be constructed
as follows.
L(Θ) ≥L+ (Θ, Λ)
1   (q)λq,i,l−τ=−2k,l i,q,τ

tr(y k,l y k,l H Ri,k,l,q,τ C −1
i,k Ri,k,l,q,τ )

(q)

−1

Ri,k,l,τ,q = C i,k,τ wi,k,q hi,l−τ X̂ k,l ,(14)

U k,l = X̂ k,l ,(15)Ũ i,k = C i,k,0 ,(16)

with respect to Ri,k,l,τ,q , U k,l and Ũ i,k . With respect to(q)
λq,i,l , (q) (q)
(q)
(q)
(q)
Fq,i,l Bq,i,l ,(17)
λq,i,l = Fq,i,l Bq,i,l /
q

(q)
(q)
Fq,i,l = p(Θ|zi,l = q)(q)

Bq,i,l =


q


q

(q)

(q)

(q)

Fq,i,l−1 ρi,q ,q ,
(q)

(q)

Bq ,i,l+1 p(Θ|zi,l+1 = q  )ρi,q,q ,

(18)
(19)

where
(q)

p(Θ|zi,l = q)1
∝ exp −
tr(y k,l+τ y k,l+τ H Ri,k,l+τ,τ,q
2
i,k,τ

C −1
i,k,τ Ri,k,l+τ,τ,q )/wi,k,q hi,l+ tr(U −1
k,l+τ C i,k,τ )wi,k,q hi,l


.

(20)

(j)

With respect to λj,i,l ,
(j)

(j)

(j)

λj,i,l = Fj,i,l Bj,i,l /
(j)

(j)


j

Fj,i,l = p(Θ|zi,l = j)

wi,k,q hi,l−τ+ tr(U −1
C
)w
h
|
+
log
|U
i,k,q i,l−τk,lk,l i,k (j)+λj,i,l (αj − 1) log hi,l − hi,l /βj(j)satisfy q λq,i,l = 1, j λj,i,l = 1 and o di,o = 1. We
denote the set of the auxiliary variables as Λ. CL+ is the sumof constant terms. tr(·) is the trace of a matrix. The equality
L(Θ) = L+ (Θ, Λ) is satisﬁed when(j)Bj,i,l =


j


j

(j)

(j)

(j)

Fj,i,l Bj,i,l ,
(j)

(21)
(j)

Fj,i,l−1 ρi,j  ,j ,
(j)

(j)

Bj  ,i,l+1 p(Θ|zi,l+1 = j  )ρi,j,j  ,

(22)
(23)

where
(j)

p(Θ|zi,l = j)
∝ exp (αj − 1) log hi,l − hi,l /βj− αj log βj − log Γ(αj ) .

j,i,l


− αj log βj − log Γ(αj )


1
−1
+
di,o − tr(C i,k,0 (J k,o + I) )
2
i,k,o

γ−M
−1
+
tr(Ũ i,k C i,k,0 )
2
 
(q)
λq,i,l log p(Z (q) )+ log |Ũ i,k | +(24)With respect to di,o ,(o)di,o = p(Θ|zi= o)1∝ exp(− tr(C i,k,0 (J ϑo ,k + I)−1 )),
2

q,i,l

2045

(25)

2015 23rd European Signal Processing Conference (EUSIPCO)

Table 1. The experimental condition of parameters
Q
15

α
1

O
30

β
β1 = 10−3 , β2 = 1010γ
3o di,o = 1.

andTable 2. SDR and SIR improvements obtained by the proposed method, the method in [7] and [8]



10−3

10−3

By setting the partial differential of L at zero, we canobtain the update rule of each parameter in Θ. Speciﬁcally,
the partial differential of L+ with respect to C i,k,0 is givenby∂L+∂C i,k,0 −1
C i,k,0 Ri,k,l,0,q y k,l y k,l H Ri,k,l,0,q C −1
1
i,k,0
=
λq,i,l2
wi,k,q hi,l
l,q

 

1
−1
(J
− U −1
w
h
d
+
I)
)
+
−
i,k,q
i,l
i,o
ϑ
,k
o
k,l
2
o

γ − M −1
(26)
Ũ i,k .
2
By setting this at zero, we obtain an algebraic Riccati equation:
C i,k,0 Ai,k,0 C i,k,0 = B i,k,0 ,
(27)
where

−1
λq,i,l wi,k,q hi,l X̂ k,lAi,k,0 =
+

l,q

+




di,o

o

B i,k,0 =C i,k,0 (method in [7]method in [8]proposed+− (J k,o + I)

−1



−1

) + (γ − M )Ũ i,k ,
(28)

−1
−1
λq,i,l wi,k,q hi,l X̂ k,l y k,l y k,l H X̂ k,l )C i,k,0 .

l,q

(29)
We can solve this equation with the same procedure proposedin [5] and obtain the update rule of C i,k,0 . For updating
C i,k,0 , we must calculate the inverse of Ũ i,k . Even if C i,k,0
equals J θi ,k + I, we can calculate the inverse of Ũ i,k which
is updated by Eq. (16). However, experimentally, Ũ i,k would
not always be a full-rank matrix in the iterations. Therefore,
we update Ũ i,k approximately as:(30)Ũ i,k ≈ C i,k,0 +  I,where is a small value. On account of the space, we omit the
update rules of the other parameters in Θ, however, they canbe also obtained analytically by setting the partial differentials
of L+ at zeros.
4. EXPERIMENTAL EVALUATION
We evaluated the performance of the present method in termsof the abilities of source separation. We used 5 mixed stereo
signals (therefore the number of the microphones M is 2) as
the experimental data, each of which we obtained by mixing three speech signals (speech of a male and two females)from the ATR database [11] and was convolved with the measured room impulse response from the RWCP database [12](in which the distance between the microphones was 11.48
cm and the reverberation time was 0 ms/380 ms). Thus, the2046

RT=0 ms
SDR
SIR-0.53 4.65
1.90
7.88
3.78 14.12

RT=380 ms
SDR SIR-1.55 3.73
-1.57 4.99
0.37 8.65

three signals were artiﬁcially located 30◦ , 90◦ and 130◦ fromthe microphones respectively. Fig. 1 shows an example ofa spectrogram of the observed signal (the reverberation timewas 380 ms). The sampling rate was 16 kHz. To compute the
STFT components of the observed signal, the STFT frame
length was set at 64 ms and a Hamming window was usedwith an overlap length of 48 ms. We set the hyperparameters as Table 1. We expected that j = 1 is an inactive state
and j = 2 is an active state. We set T as 0 for the anechoic
mixed signals, as 3 for echoic ones. d was√ initially set asshown in Fig. 2 (a).
C
was
set
to
1/
M × I, C i,k,1:T
i,k,0√
were set to 10−2 / M × I. W was initially randomized. H(j)(j)
was set as 1 initially. ρi was set as ρi,1 = (0.9, 0.1) and
(j)

(q)

ρi,2 = (0.1, 0.9), all components of ρi were set as 1/Q.
We set ρ(o) individually for each source so that directions of
3 sources did not overlap each other. The parameter estimation algorithm was run for 50 iterations. In order to avoid an
undesirable local optima, we iterated the proposed algorithm
25 times with setting T as 0, then gradually increased T up
to 3 according to the iteration. Moreover, in order to avoid
huge numerical errors, we ﬁxed W when T > 0. The numerical errors occurred probably because we ﬁtted the product of
W i,k and C i,k,1:T to observed signals, and so each of W i,kand C i,k,1:T was arbitrary scaled. We chose the methods proposed in [7] and [8] as comparisons. The method of [7] does
not use the information of DOAs, and that of [8] models aspatial correlation matrix by using a weighted sum of J . The
separated signal was obtained by Wiener ﬁltering. As evaluation measures, we used the signal-to-distortion ratio (SDR)
and the signal-to-interference ratio (SIR) [13]. The SDR andSIR are expressed in decibels (dB), and a higher SDR(/SIR)
indicates superior quality. The SDRs and SIRs were calculated with magnitude spectrograms. The averages of SDRs
and SIRs of the mixed signals were -1.98 and -1.04[dB] when
RT=0 ms, -4.77, -0.97[dB] when RT=380 ms respectively.
Table 2 shows the average SDR and SIR improvements
obtained by our previous and present methods. The average
SDRs and SIRs obtained with the present method were superior to these obtained with our previous methods. Fig. 2 (b)and (c) show examples of d obtained by the proposed method
when RT=0 ms and RT=380 ms respectively. The arrows
show the actual direction of the sources. We can see that the
DOA of the sources were estimated roughly by the proposed
method. Fig. 3 shows (a)a spectrogram of an anechoic source
signal, (b)that of a separated and dereverbed signal with theproposed method and (c)the result of λ(j) which corresponds
to the estimated activity of the source obtained with the mixed
signal showed in Fig. 1. We expected q = 1 was an inactive
state by setting the hyperparameters of the gamma distribution properly (see Table 1), and the result shows the voiceactivity was roughly detected.

2015 23rd European Signal Processing Conference (EUSIPCO)

Fig. 1. An example of the spectrogram of the mixed signal(RT=380 ms).
(a)
d

(b)

RT=0msd(c)

RT=380msdFig. 2. (a)First condition of d, (b)an example of the estimation
results of d obtained with the mixed signal (RT=0 ms) and(c)that obtained with the mixed signal (RT=380 ms). d ofeach source is colored differently. d corresponds to how likely
the source is at the direction.

5. CONCLUSION
In this paper, we incorporate the DOA mixture model into themultichannel factorial hidden Markov model. Parameter estimation allows us to solve the problems of BSS, VAD, dereverberation and DOA estimation simultaneously. The average of
SDRs obtained by the proposed method was 1.94 dB higher
than that obtained by our previous method.
REFERENCES[1] A. Hyvàrinen et al., Independent Component Analysis,
John Wiley & Sons, 2001.
[2] D. D. Lee, and H. S. Seung, “Learning the parts ofobjects with nonnegative matrix factorization,” Nature,
vol. 401, pp.788–791, 1999.
[3] P. Smaragdis, and J. C. Brown, “Non-negative matrix
factorization for polyphonic music transcription,” WASPAA 2003, pp. 177–180, Oct. 2003.
[4] A. Ozerov, and C. Févotte, “Multichannel nonnegative
matrix factorization in convolutive mixtures for audio
source separation,” IEEE Trans. ASLP, vol. 18, no. 3,
pp. 550–563, 2010.
[5] H. Sawada et al., “Efﬁcient algorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” ICASSP 2012, pp. 261–264, 2012.

2047

Fig. 3. (a)An example of a spectrogram of a source signal,(b)that of the separated and dereverbed signal obtained by the
proposed method and (c)the estimation result of λ(j) whichcorresponds to source’s activity obtained with the mixed signal showed in Fig. 1. Black indicates higher value and thestate is more likely to be assigned at that time in (c).

[6] T. Higuchi, et al., “A uniﬁed approach for underdetermined blind signal separation and source activity detection by multichannel factorial hidden Markov models,”Interspeech 2014, pp. 850–854, 2014.
[7] T. Higuchi and H. Kameoka, “Joint audio source separation and dereverberation based on multichannel factorial hidden Markov model,” MLSP 2014.
[8] T. Higuchi and H. Kameoka, “Uniﬁed approach forunderdetermined BSS, VAD, dereverberation and DOA
estimation with multichannel factorial HMM,” GlobalSIP 2014, pp. 725-729, 2014.
[9] J. Nikunen and T. Virtanen, “Multichannel audio separation by direction of arrival based spatial covariancemodel and non-negative matrix factorization,” ICASSP
2014, pp. 6727–6731, 2014.
[10] H. Kameoka et al., “Blind separation of inﬁnitely manysparse sources,” IWAENC 2012, pp. 1–4, 2012.
[11] A. Kurematsu et al., “ATR Japanese speech databaseas a tool of speech recognition and synthesis,” SpeechCommunication, pp. 357–363, 1990.
[12] S. Nakamura et al., “Acoustical sound database in realenvironments for sound scene understanding and handsfree speech recognition,” LREC 2000, pp. 965–968,
2000.
[13] E. Vincent et al., “Performance measurement in blindaudio source separation,” IEEE Trans. ASLP, pp. 1462–
1469, 2006.