JOINT SEPARATION AND DEREVERBERATION OF REVERBERANT MIXTURES WITH
DETERMINED MULTICHANNEL NON-NEGATIVE MATRIX FACTORIZATION
Hideaki Kagami† , Hirokazu Kameoka‡ , Masahiro Yukawa†† Dept. Electronics and Electrical Engineering, Keio University, Japan.
‡ NTT Communication Science Laboratories, NTT Corporation, Japan.
sum of a limited number of basis spectra scaled by time-varying
amplitudes. Multichannel NMF (MNMF) is an extension of thisapproach to a multichannel case in order to allow for the use of
spatial information as an additional clue for separation. It can also
be viewed as an extension of frequency-domain BSS that allows theuse of spectral templates as a clue for both frequency-wise source
separation and permutation alignment.
The original MNMF [3–5] was formulated under an underdetermined problem setting where sources can outnumber microphones.
A determined version of MNMF (DMNMF) [6–8] specializes in theoverdetermined case and can be implemented without using any matrix inversion, shortening the computation time more than 30 times
than the underdetermined MNMF. In [7, 8], a theoretical relation of
DMNMF to IVA is discussed, which has naturally allowed for the incorporation of a fast update rule of the separation matrix developed
for IVA, called the “iterative projection (IP)” [11], into the parameter optimization process of DMNMF. It has been shown that this has
contributed not only to further accelerating the entire optimization
process but also improving the separation performance.
One drawback as regards methods using instantaneous mixture
models is that they are weak against long reverberation. It has beenshown in [6, 12] that source separation of highly reverberant mixtures can be solved effectively by using frequency-domain convolutive mixture models and that the parameters of frequency-domain
convolutive mixture models can be efficiently estimated by iteratively updating the separation matrix, the dereverberation filters andthe spectral parameters of each source. In this paper, we take thisapproach to develop an extension of DMNMF using a frequencydomain convolutive mixture model, where the optimization process
consists of iteratively updating (i) the NMF parameters using the
majorization-minimization (MM) algorithm, (ii) the separation matrix using iterative projection (IP), and (iii) the dereverberation filters
using multichannel linear prediction.

ABSTRACTThis paper proposes an extension of multichannel non-negative matrix factorization (MNMF) that simultaneously solves source separation and dereverberation. While MNMF was originally formulated
under an underdetermined problem setting where sources can outnumber microphones, a determined counterpart of MNMF, which wecall the determined MNMF (DMNMF), has recently been proposedwith notable success. This approach is particularly notable in thatthe optimization process can be more than 30 times faster than theunderdetermined version owing to the fact that it involves no matrix
inversion computations. One drawback as regards all methods basedon instantaneous mixture models, including MNMF, is that they are
weak against long reverberation. To overcome this drawback, thispaper proposes an extension of DMNMF using a frequency-domain
convolutive mixture model. The optimization process of the proposed method consists of iteratively updating (i) the spectral parameters of each source using the majorization-minimization algorithm,(ii) the separation matrix using iterative projection, and (iii) the dereverberation filters using multichannel linear prediction. Experimental results showed that the proposed method yielded higher separation performance and dereverberation performance than the baseline
method under highly reverberant environments.
Index Terms— Blind source separation, blind dereverberation,
non-negative matrix factorization, independent component analysis
1. INTRODUCTION
Blind source separation (BSS) is a technique for separating out individual source signals from microphone array inputs when the transfer characteristics between sources and microphones are unknown.
The frequency-domain BSS approach provides the flexibility of allowing us to utilize various models for the time-frequency representations of source signals and/or array responses. For example, independent vector analysis (IVA) [1, 2] allows us to efficiently solve
frequency-wise source separation and permutation alignment in ajoint manner by assuming that the magnitudes of the frequency components originating from the same source tend to vary coherently
over time.
With a different approach, multichannel extensions of nonnegative matrix factorization (NMF) have attracted a lot of attention
in recent years [3–8]. NMF was originally applied to monauralsource separation tasks [9,10]. The idea is to approximate the power(or magnitude) spectrogram of a mixture signal, interpreted as a
non-negative matrix, as the product of two non-negative matrices.
This amounts to assuming that the power spectrum of a mixture signal observed at each time frame can be approximated by the linear

2. PROBLEM FORMULATION
We start by introducing a frequency-domain convolutive mixture
model employed in [6, 12]. Let us consider a situation in which
signals emanate from M sources and are captured by M microphones. Let xi (f, n) be the short-time Fourier transform (STFT)of the signal observed at the i-th microphone, where f and n are
the frequency and time indices, respectively. Also let x(f, n) =[x1 (f, n), ..., xM (f, n)]T ∈ CM be a vector of observed data. We
use sj (f, n) to denote the j-th source conponents and s(f, n) =[s1 (f, n), ..., sM (f, n)]T ∈ CM to denote a vector of M sourcecomponents. Under a determined problem setting, many conventional BSS systems employ frequency-domain instantaneous mixture model:This work was supported by JSPS Grants-in-Aid (15K06081,
15K13986, 15H02757, 17H01763).

978-1-5386-4658-8/18/$31.00 ©2018 IEEEs(f, n) = W H (f )x(f, n)

31

(1)

ICASSP 2018

where W H (f ) is called the separation matrix. However, in a highly
reverberant condition where the length of the room impulse responses can be longer than the STFT frame length, the ability ofthis model to separate sources will be limited. In this paper, we use
a separation system that has a multichannel finite-impulse-response
form in the time-frequency domain [6, 12] such that:three minimization steps:s(f, n) = W (f, 0)x(f, n) +

H





W (f, n )x(f, n − n ),where W (f, n ), 0 ≤ n ≤ N are the coefficient matrices of size
M ×M , and (·)H stands for Hermitian transpose. When W H (f, 0) is
invertible, (2) can be written equivalently by the following process:GH (f, n )x(f, n − n ),The update rules for Θ̂V can be derived by the MM algorithm. Dropping the constant terms with respect to Θ̂V from (8), we obtain


K


|sj (f, n)|2
hj,k (f )uj,k (n) + K
C1 (ΘV ):=log,
k=1 hj,k (f )uj,k (n)
f,n,j
k=1

(3)

n =1

s(f, n) = W H (f, 0)y(f, n),
H



H

−1

(4)
H

(12)



where G (f, n ) := −(W (f, 0)) W (f, n ). Equation (3) canbe seen as a dereverberation process of the observed mixture signal x(f, n), whereas (4) can be seen as an instantaneous demixing
process of the dereverberated mixture signal y(f, n).

where vj (f, n) in (8) is replaced by (7). To minimize this function,
we can design a majorization function of C1 with a convenient form
[13–15]K
1hj,k (f )uj,k (n)−ξ(f, n)
C1+(ΘV , Θ̄)=
ξ(f, n)Let us assume that sj (f, n) is a random variable:sj (f, n) ∼ NC (sj (f, n) | 0, vj (f, n))

f,n,j

k=1

λ2k (f, n), (13)
hj,k (f )uj,k (n)where Θ̄ := {ξ, λ} with ξ = {ξ(f, n)} and λ = {λk (f, n)}. It canbe verified that the majorization function C1+ satisfies the following
properties:
1. C1 (ΘV ) ≤ C1+ (ΘV , Θ̄)2.

hj,k (f )uj,k (n),K+ log ξ(f, n) + |sj (f, n)|2

denotes the complex Gaussian distribution. We further assume thatthe power spectral density vj (f, n) can be modeled as
K


k=1

(5)

and that sj (f, n) and sj  (f  , n ) are statistically independent when(f, n, j) = (f  , n , j  ). Here,|sj (f, n)|21
exp −NC (sj (f, n) | 0, vj (f, n)) ∝(6)
vj (f, n)vj (f, n)vj (f, n) =(11)

3.1. Update procedure for Θ̂VNΘ̂G ← argmin Q(Θ̂V , Θ̂W , ΘG ).

In the next subsections, we focus on each step and derive the update
rules.

y(f, n) = x(f, n) −

(10)

ΘG



(2)


Θ̂W ← argmin Q(Θ̂V , ΘW , Θ̂G )
ΘW

n =1
H

(9)

ΘV



N


H

Θ̂V ← argmin Q(ΘV , Θ̂W , Θ̂G )(7)

C1 (ΘV ) = min C1+ (ΘV , Θ̄).
Θ̄

(14)
(15)

Here, (14) holds with equality and, at the same time, the minimumin (15) is achieved when

k=1

where hj,k (f ) ≥ 0 is the (j, k) element of the basis matrix anduj,k (n) ≥ 0 is the (j, k) element of the activation matrix for the jth source. Multichannel source separation methods using the power
spectrogram model (7) or its variants are called MNMF [3,5–8]. The
negative log-likelihood of the signal yi (f, n) is given by





log det W H (f, 0)
Q(ΘV , ΘW , ΘG ) := − 2N

ξ(f, n) =Khj,k (f )uj,k (n),

(16)

k =1

λk (f, n) = hj,k (f )uj,k (n)/Khj,k (f )uj,k (n).

(17)

k =1

The function C1 is expected to be minimized indirectly by repeating
the following two steps.

f|sj (f, n)|2+, (8)log vj (f, n) +vj (f, n)

1. Minimize C1+ with respect to Θ̄ by (16) and (17), whichmakes C1 (ΘV ) = C1+ (ΘV , Θ̄).

f,n,j

where ΘV := {H, U } with H = {hj,k (f )} and U = {uj,k (n)},ΘW := {W H (f, 0)}, and ΘG := {GH (f, n )}.

2. Minimize C1+ with respect to ΘV .
For the second step, making the partial derivatives of C1+ with respect to each element of H and U to be zero yields the following
update rules for H and U :

2
−2


j  |sj (f, n)| uj ,k (n)rj  (f, n)
hj,k (f ) ← hj,k (f ), (18)
−1j  uj ,k (n)rj  (f, n)

3. OPTIMIZATION PROCESSuj,k (n) ← uj,k (f )The cost function (8) can be iteratively decreased by using a coordinate descent method in which each iteration comprises the following

where rj  (f, n) :=32

j  |sj (f, n)|




K

2

hj  ,k (n)rj−2
 (f, n)−1j  hj ,k (n)rj  (f, n)k =1 hj ,k (f )uj ,k (n).








, (19)

3.2. Update procedure for Θ̂W
Dropping the constant terms with respect to Θ̂W from (8), we obtain





C2 (ΘW ) := −2N
log det W H (f, 0)
+N



SIR
DRR

f

wH
j (f )Σy/vj (f ) w j (f ),(20)
Table 2. Average computational time (sec.)
Proposed (IP)
371.1974
Proposed (FICA) 396.7370
Baseline [7, 8]
303.5906
Sequential
321.8167f,jwhere wj (f ) is the j-th column vector of W (f, 0) and Σy/vj (f ) = y(f,n)y H (f,n)
1, which is assumed positive definite (i.e.,nN
vj (f,n){y(f, n)}n is a linearly independent set). As noted above, when
the parameters Θ̂V and Θ̂G are fixed, the cost function (8), whichis equal to (20), can be seen as an instantaneous demixing processof the dereverberated mixture signal y(f, n). To update Θ̂W , several determined BSS techniques are available, such as the natural
gradient method [16], FastICA (FICA) [17] and IP [2].
IP [2] is a block coordinate descent algorithm, consisting of sequentially optimizing one column vector of W H (f, 0) at a time
while keeping the other column vectors fixed. Differentiating
C2 (ΘW ) with respect to the conjugate w∗j (f ) of wj (f ), and setting
the result to zero, we have∂(21)
log det W H (f, 0) = 0.
Σy/vj (f ) wj (f ) − 2∂w ∗j (f )where g m (f, n ) is the m-th column of G(f, n ). By using g(f ),
 
H


we can rewrite the term Nn =1 G (f, n )x(f, n − n ) in (26) asN

wHj (f )Σy/vj (f ) wj (f ) = 1,

GH (f, n )x(f, n − n ) = X(f, n)g ∗ (f ),(28)

X(f, n) = [I ⊗ xT (f, n − 1) I ⊗ xT (f, n − 2)

(29)

n =1

where
T

for j = j,f,nwj (f )
wj (f ) ← ,Hwj (f )Σy/vj (f ) wj (f )

(25)

(30)

Our objective in this part is now to minimize this function with respect to g ∗ (f ).
Equation (30) is readily solved because the cost function isquadratic with respect to g ∗ (f ). We calculate the partial derivativesof C3 (ΘG ) to be zero, the update rule for g ∗ (f ) is obtained as

−1
 H
∗
X (f, n)Σw/v(f,n) X(f, n)g (f ) ←A solution to (22) and (23) can be found by performing the following updates(24).

× Σw/v(f,n) (x(f, n) − X(f, n)g ∗ (f )) .

(23)

wj (f ) ← (W H (f, 0)Σy/vj (f ) )−1 ej ,M ×M 2 N 

Here, ⊗ stands for Kronecker product. Substituting (28) into (26),the cost function is rewritten asH
C3 (ΘG ) =(x(f, n) − X(f, n)g ∗ (f ))

(22)




· · · I ⊗ x (f, n − N )] ∈ C∂
det(A) = A−T det(A),By using the matrix derivatives formula ∂A(21) can be rearranged in the following simultaneous vector equations:wHj  (f )Σy/vj (f ) wj (f ) = 0

Table 1. Input average SIRs and DRRs (dB)Male + Male Male + Female Female + Female
0.21900.11600.07382.91722.20314.4501n

for all f and j. Here, ej denotes the j-th column of the M × M
identity matrix I.

×





H

X (f, n)Σw/v(f,n) x(f, n) ,

(31)

n

for all f .

3.3. Update procedure for Θ̂GDropping the constant terms with respect to Θ̂G from (8), we obtain
2


N




H



x(f, n) −C3 (ΘG ):=
G
(f,
n
)x(f,
n
−
n
)
,



f,n 
n =1

3.4. Summary of optimization process
Overall, the proposed algorithm is summarized as follows.
1. Initialize H, U , W H (f, 0) and GH (f, n ).
2. Update H, U using (18), (19) for all f, n, j.
3. Update W H (f, 0) using (24), (25) for all f, j.
4. Update GH (f, n ) using (31) for all f .
Repeat 2–4 until convergence.

Σw/v(f,n)(26)

where x Σw/v(f,n) :=w/v(f,n) x with Σw/v(f,n) = w(f )wH (f )which is assumed positive definite, and O denotesjvj (f,n)the zero matrix. It is clear that, for each f , all elements of thematrices {G(f, n )} minimizing (26) are mutually dependent. To
update G(f, n ) for each f independently, we vectorize {G(f, n )}as follows:xH Σ

4. EXPERIMENTAL RESULTS
We evaluate the effect of the convolutive model under highly reverberant environments and the difference between computational timeby separating matrix update method. We quantitatively compared
the source separation performance of the proposed approach using IP
and FICA, conventional DMNMF and sequential method (Sequential). All of the examples use the same two-input four-output impulse response, which was measured in an anechoic chamber whereg(f ) := vec({G(f, n )}):= [g T1 (f, 1), ..., g TM (f, 1), g T1 (f, 2), ..., g TM (f, 2), ...,
2

g TM (f, N  − 1), g T1 (f, N  ), ..., g TM (f, N  )]T ∈ CM N ,
(27)

33

Proposed (FICA)

20

20
15
10

Input

5

15

Baseline [7, 8]
Sequential

0
-5
0

10

1000

14

4000

5000

6000

7000

Proposed (FICA)
Proposed (IP)15

DRR [dB]

Male + Male Male + Female Female + Female
Fig. 1. Average SIR for each gender pair.

16

3000

Fig. 3. Average SIR for each frequency.
20

0

2000

Frequency [Hz]5

Proposed (IP)
Proposed (FICA)
Baseline [7, 8]
Sequential

10
5
0
-5

-10
0

Sequential
Input
Baseline [7, 8]
1000

2000

3000

4000

5000

6000

7000

Frequency [Hz]12

Output DRR [dB]

Proposed (IP)25

SIR [dB]

Output SIR [dB]2530

Proposed (IP)
Proposed (FICA)
Baseline [7, 8]
Sequential

Fig. 4. Average DRR for each frequency.

10
8that performs dereverberation and DMNMF sequentially. Figs. 1
and 2 show average output SIRs and DRRs and the 95% confidence
interval for each gender pair. One factor of the performance difference is that the filter length N  was fixed as described above in allcases. Figs. 3 and 4 show average output SIRs and DRRs for eachfrequency. Fig. 3 shows that the proposed approach attains higher
SIR for all frequencies. Referring to Fig. 4, the proposed approachattains significantly higher DRRs than the baseline method over thewide range below 3 kHz, although its DRRs are lower (significantly
only for frequency between 6.2 – 6.7 kHz) above 3 kHz. Overall,these results show that the proposed approaches perform better thanthe conventional DMNMF and sequential methods.

6
4
2
0

Male + Male Male + Female Female + Female
Fig. 2. Average DRR for each gender pair.

the reverberation time was 0.6 sec. With this impulse response, we
mixed two speech signals into four mixtures. The two speech signals
of male or female speakers, taken from the ATR speech database,were sampled at 16 kHz and band limited to the 50 Hz to 7 kHz frequency range. We generate 10 speech combinations for each gender
pair, male+male, male+female and female+female. The average input Signal-to-Interference ratios (SIRs) and Direct-to-Reverberation
ratios (DRRs) are shown in Table 1. Time-frequency representations were obtained using the polyphase filterbank analysis with aframe length of 32 ms and a hop size of 8 ms. For the DMNMF, weperformed time-frequency analysis with a frame length of 256 ms
and a hop size of 64 ms. The filter length N  was set as follows:N  = 25 for Pf < 0.8; N  = 20 for 0.8 ≤ Pf < 1.5; N  = 15 for
1.5 ≤ Pf < 3; N  = 10 for Pf ≥ 3; ,where Pf is the frequency in
kHz of the f -th frequency bin. The NMF basis number K was set at
20. The iterative algorithm was run for 20 iterations. For each stepof update the separation matrix W H (f, 0), the FICA algorithm was
run for 100 iterations.
Table 2 shows IP is faster then FICA algorithm, where Baseline
refer to conventional DMNMF [7, 8]. The “Sequential” is a method

5. CONCLUSION
One drawback as regards all methods based on instantaneous mixture models, including MNMF, is that they are weak against longreverberation. To overcome this drawback, this paper proposed an
extension of DMNMF using a frequency-domain convolutive mixture model, which allows us to solve source separation and dereverberation sumultaneously. The optimization process of the proposed
method consists of iteratively updating (i) the spectral parameters ofeach source using the majorization-minimization algorithm, (ii) the
separation matrix using IP, and (iii) the dereverberation filters using
multichannel linear prediction. Experimental results showed that the
proposed method yielded higher separation performance and dereverberation performance than the baseline method under highly reverberant environments.

34

[1] T. Kim, T. Eltoft, and T.-W. Lee, “Independent vector analysis: An extension of ica to multivariate components,” in International Conference on Independent Component Analysis andSignal Separation. Springer, 2006, pp. 165–172.

[14] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama,“Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with βdivergence,” in Machine Learning for Signal Processing(MLSP), 2010 IEEE International Workshop on. IEEE, 2010,
pp. 283–288.

[2] N. Ono and S. Miyabe, “Auxiliary-function-based independent
component analysis for super-gaussian sources,” in LVA/ICA.
Springer, 2010, pp. 165–172.

[15] C. Févotte and J. Idier, “Algorithms for nonnegative matrix
factorization with the β-divergence,” Neural computation, vol.
23, no. 9, pp. 2421–2456, 2011.

[3] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “New formulations and efficient algorithms for multichannel nmf,” in Applications of Signal Processing to Audio and Acoustics (WASPAA), 2011 IEEE Workshop on. IEEE, 2011, pp. 153–156.

[16] S. Amari, A. Cichocki, and H. H. Yang, “A new learning algorithm for blind signal separation,” in Advances in neural
information processing systems, 1996, pp. 757–763.

6. REFERENCES[17] A. Hyvärinen, J. Karhunen, and E. Oja, Independent component analysis, vol. 46, John Wiley & Sons, 2004.

[4] N.Q.K. Duong, E. Vincent, and R. Gribonval, “Underdetermined reverberant audio source separation using a fullrank spatial covariance model,” IEEE Transactions on Audio,
Speech, and Language Processing, vol. 18, no. 7, pp. 1830–
1840, 2010.
[5] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix
factorization in convolutive mixtures for audio source separation,” IEEE Transactions on Audio, Speech, and Language
Processing, vol. 18, no. 3, pp. 550–563, 2010.
[6] H. Kameoka, T. Yoshioka, M. Hamamura, J. Le Roux, and
K. Kashino, “Statistical model of speech signals based on composite autoregressive system with application to blind source
separation.,” in LVA/ICA. Springer, 2010, pp. 245–253.
[7] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Efficient multichannel nonnegative matrix
factorization exploiting rank-1 spatial model,” in Acoustics,
Speech and Signal Processing (ICASSP), 2015 IEEE International Conference on. IEEE, 2015, pp. 276–280.
[8] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Determined blind source separation unifying
independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 24, no. 9, pp. 1626–1641, 2016.
[9] P. Smaragdis, “Non-negative matrix factorization for polyphonic music transcription,” in Proc. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, 2003,
pp. 177–180.
[10] 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.
[11] N. Ono, “Stable and fast update rules for independent vector
analysis based on auxiliary function technique,” in Applications of Signal Processing to Audio and Acoustics (WASPAA),
2011 IEEE Workshop on. IEEE, 2011, pp. 189–192.
[12] T. Yoshioka, T. Nakatani, M. Miyoshi, and H. G. Okuno,“Blind separation and dereverberation of speech mixtures by
joint optimization,” IEEE Transactions on Audio, Speech, andLanguage Processing, vol. 19, no. 1, pp. 69–84, 2011.
[13] H. Kameoka, M. Goto, and S. Sagayama, “Selective amplifierof periodic and non-periodic components in concurrent audio
signals with spectral control envelopes,” Proc. SIG Tech. Reports on Music and Computer of IPSJ, pp. 77–84, 2006.

35