2016 24th European Signal Processing Conference (EUSIPCO)

Reverberation-robust underdetermined source
separation with non-negative tensor double
deconvolution
Naoki Murata∗ , Hirokazu Kameoka∗† , Keisuke Kinoshita† , Shoko Araki† , Tomohiro Nakatani† ,Shoichi Koyama∗ and Hiroshi Saruwatari∗∗ 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-shi, Kanagawa 243-0198, Japan

Abstract—Source separation using an ad hoc microphonearray can be useful for enhancing speech in such applications as teleconference systems without the need to prepare
special devices. However, the positions of the sources (and themicrophones when using an ad hoc microphone array) canchange during recording, thus violating the commonly madeassumption in many source separation algorithms that the mixing
system is time-invariant. This paper proposes an extension of themultichannel nonnegative matrix factorization (NMF) approachto deal with the problem of underdetermined source separationin time-variant reverberant environments. The proposed methodmodels the mixing system as a non-negative convolutive mixture
based on the concept of a “semi-time-variant system” to handle
the reverberation in a room as well allowing for relatively smallchanges in the source/microphone positions. It also models thepower spectrogram of each sound source using the convolutive
NMF model to consider the local dynamics of speech.

I. I NTRODUCTION
An ad hoc microphone array is a microphone array thatuses built-in microphones in portable devices such as laptops
and smartphones as spatially distributed sensors in an ad hoc
fashion. It is particularly noteworthy in that it allows a flexible
and convenient way to acquire audio data compared withtraditional microphone arrays, which usually require special
devices. Since a traditional microphone array is typicallyarranged so that the built-in microphones are located close
to each other, the interchannel difference in arrival time (orphase difference) of each sound source provides an importantclue in such applications as sound source localization andsource separation. By contrast, since an ad hoc microphonearray allows each microphone to be spatially distributed, theinterchannel level difference of each sound source can also be
a useful clue. Motivated by this, this paper proposes a source
separation method specifically tailored for ad hoc microphonearrays that relies on the interchannel level difference of eachsound source.
The problem of separating the signals of individual soundsources from the signals observed with a microphone array canbe considered an inverse problem. In an underdetermined casewhere the sources outnumber the microphones, this inverse
problem generally has an infinite number of solutions. Thus,

978-0-9928-6265-7/16/$31.00 ©2016 IEEE

some reasonable assumption is usually needed to limit therange of possible solutions. In recent years, multichannel
extensions of non-negative matrix factorization (NMF) haveattracted particular attention after being proposed as a powerfulapproach for underdetermined source separation [1]–[4]. Notethat multichannel extensions of NMF for an overdeterminedcase have also been proposed with notable success [6], [7].
NMF was originally applied to monaural source separation
where the magnitude (or power) spectrogram of a mixture signal, interpreted as a non-negative matrix, is factorized into the
product of two non-negative matrices. This can be interpreted
as approximating the observed spectra at each time frame as alinear sum of basis spectra scaled by time-varying amplitudes,
and amounts to decomposing the observed spectrogram into
the sum of low rank spectrograms. Multichannel NMF is anextension of NMF to a multichannel input, which assumes thepower spectrogram of each underlying sound source to havea low rank structure.

While most of these approaches are formulated on theassumption that the mixing system is time-invariant, thisassumption does not necessarily hold when the relative positions of sources and microphones are likely to change during
recording. In such situations, methods based on the timeinvariant mixing system do not work very satisfactorily.

This paper proposes an extension of the multichannel NMF
approach to deal with the problem of underdetermined source
separation in time-variant reverberant environments. The proposed method models the mixing system as a non-negative
convolutive mixture based on the concept of the “semi-timevariant system” to take account of the reverberation in aroom as well as to allow for relatively small changes inthe source/microphone positions. The method also models thepower spectrogram of each sound source using the convolutive
NMF model [13] to take account of the local dynamics in thetime-frequency components of speech.

16482016 24th European Signal Processing Conference (EUSIPCO)

II. P ROPOSED MODEL
A. Convolutive mixture in time-frequency domain
We start by describing the mixing system with I sources andJ microphones as a convolutive mixture in the time-frequency
domain. When the transfer systems between the sources andthe microphones are linear time-invariant and their impulse
responses are sufficiently shorter than the window length of theshort-time Fourier transform (STFT), the signals obtained by
the microphones can be represented as an instantaneous mixture of source signals in the STFT domain, which is adopted
in many conventional BSS methods including multichannel
NMF [1]–[4]. However, when reverberation comes into play,
the length of the room impulse responses can be longer thanthe STFT window. In such cases, the signals observed at the
microphones can be well approximated by the convolutive
mixture model [8], [9]∑∑yk,l =
ai,k,n si,k,l−n ,
(1)
i

n

where i, k denote the source and frequency indices and l, ndenote the time frame indices, respectively. yk,l ∈ CJ is avector consisting of the complex time-frequency components
observed at the J microphones. We hereafter use j to denotethe microphone index. ai,k,n represents the steering vectorof each source, which corresponds to the contributions of
the direct component (when n = 0) and the reverberant
components arriving at the microphones with a time delayof n frames (when n ̸= 0). si,k,l represents the complex timefrequency component of source i.
B. Nonnegative tensor double deconvolution
In a general scenario, the positions of the sources (andthe microphones in an ad-hoc microphone array scenario)
can change during recording. To handle such situations, wemust consider a time-variant system. To achieve a source
separation algorithm that is robust against changes in the
acoustic environment, one possible approach involves treating
a time-variant factor as a latent variable to be marginalized out
so that the algorithm becomes less sensitive to that factor. [11]proposes a speech dereverberation method robust to speaker’s
movements that follows this idea. This idea is also adopted ina blind source separation method for an ad hoc microphonearray [12], which is designed to be robust against the sampling
rate mismatch of the microphones of the array. The proposed
method uses a similar approach to develop a robust source
separation algorithm that is designed to be less sensitive to the
movements of the microphones and sources. Specifically, when
the positions of sources and microphones are time-variant,i.e., the steering vector ai,k,n depends on time l, the mixing
process in (1) is described as a time-variant system∑∑
yk,l =
ai,k,n,l si,k,l−n .
(2)
i

n

We assume that the time-frequency component of each source
independently follows a complex Gaussian distribution, as in[15], i.e., si,k,l ∼ NC (0, Pi,k,l ) where Pi,k,l is the power

spectrogram of source i. Thus, the time-frequency components
observed at the J microphones follow( ∑∑)yk,l ∼ NC 0,
Pi,k,l−n ai,k,n,l aH(3)i,k,n,l .
inHere, we decompose the time-variant steering vector ai,k,n,l
into its magnitude and phase parts  jϕ|a1,i,k,n,l |0
e 1,i,k,n,l



..
..
ai,k,n,l =  . (4).
.
|aJ,i,k,n,l |0

ejϕJ,i,k,n,l

We assume that the magnitude part of the steering vector is lesssensitive to the changes in the acoustic environment made for
instance by a slight movement of the microphones/sources andassume a mixing system in which |aJ,i,k,n,l | is time-invariant
whereas ϕj,i,k,n,l is time-variant. We hereafter refer to thistype of mixture process as semi-time-variant system. Thus, byassuming |aj,i,k,n,l | = Aj,i,k,n , (4) can be rewritten as  jϕe 1,i,k,n,l
A1,i,k,n
0



..
..
(5)
ai,k,n,l = 
.

.
.
|

0

AJ,i,k,nejϕJ,i,k,n,l
{z
}|
{z
}
Ai,k,nψi,k,n,l

By substituting Ai,k,n ψi,k,n,l for ai,k,n,l in (3), we obtain thefollowing distribution( ∑)Hyk,l ∼ NC 0,
Pi,k,l−n Ai,k,n ψi,k,n,l ψi,k,n,lAHi,k,n . (6)i,nWhile the phase part of the steering vector is one of themost useful clues to source separation for traditional arrayswith spatially proximate microphones, the magnitude part ofthe steering vector, which can be seen as the interchannel level
differences between the direct and reverberant componentsof each source, can also be a useful clue when it comesto an ad hoc array with spatially distributed microphones(imagine a situation where each microphone is located very
close to a different speaker). Here, we instead consider thephase part of the steering vector to be a nuisance parameter
since it can vary sensitively according to the movements ofthe microphones/sources. Hence, we treat ϕj,i,k,n,l as a latentvariable to be marginalized out and make the following two
assumptions:′′• ϕj,i,k,n,l and ϕj ′ ,i,k,n,l (j ̸= j or l ̸= l ) are statisticallyindependent.
• ϕj,i,k,n,l follows a uniform distribution in [0, 2π).
H
By marginalizing out ϕj,i,k,n,l , E[ψi,k,n,l ψi,k,n,l] becomes anidentity matrix. Thus, we obtain the following distribution
)
( ∑∑
(7)
yj,k,l ∼ NC 0,
Pi,k,l−n A2j,i,k,n .
inUntil now, no assumptions have been made as regards Pi =(Pi,k,l )K×L . The conventional multichannel NMF [1]–[3]models the power spectrogram Pi,k,l as a product of two nonnegative matrices. This can be interpreted as representing the1649

2016 24th European Signal Processing Conference (EUSIPCO)

power spectrum of the source at each time frame as the nonnegative combination of spectrum templates. Although thisassumption may be relatively accurate for musical instrument
sounds, it may be inadequate when it comes to speech since
speech is characterized not only by the instantaneous spectra but also by the continuous transitions of these spectra
over time. Therefore, we consider it reasonable to treat the
concatenation of the spectra of multiple frames rather thana single-frame spectrum as a basic element constituting theentire spectrogram of speech. Motivated by this, we employ
the convolutive NMF model introduced in [13], [14] to model
the power spectrogram of each source. Namely, the power
spectrogram of each source is modeled as a convolution overtime of spectrogram templates and the corresponding temporal
activation functions
Pi,k,l =−1∑ T∑

Wi,k,m,τ Hi,m,l−τ ,(8)m τ =0where m and τ denote the indices of the spectrogram templates
and time frame, respectively. Wi,k,m,τ and Hi,m,l represent
the ith spectrogram template and the corresponding temporal
activation function. Note that this model reduces to the regular
NMF model when T = 1.
Our proposed model involves two convolutions, one describing the reverberation process, and the other describing
the power spectrogram model. Hence, the source separation
problem based on the proposed model is formulated as a
double deconvolution problem. We thus call our proposed
method Non-negative Tensor Double Deconvolution, (NTDD).
III. PARAMETER ESTIMATION
From (7), the negative log-likelihood function of the unknown parameters is given by()∑∑
CM L =
log NC yj,k,l |0,
Pi,k,l−n A2i,j,k,n
j,k,l
c

=

∑

j,k,l
c

(

i,n

dIS |yj,k,l | |2∑)
Pi,k,l−n A2i,j,k,n ,

(9)

i,n

where = denotes equality up to constant terms. dIS (y|x)
denotes the Itakura-Saito divergence between x and y [15].
Therefore, the parameter estimation using the maximum
likelihood estimation amounts to the minimization of theelement-wise Itakura-Saito divergence between the obtained
power spectrogram Y
j,k,l and the modeled power spec∑
trogramX=j,k,li,n Pi,k,l−n Ai,j,k,n where Pi,k,l =∑ ∑mτ Wi,k,m,τ Hi,m,l−τ .
Seen from a model-fitting perspective, it would be naturaland interesting to consider employing divergence measuresother than the Itakura-Saito divergence. Here, we generalize
the divergence measure using the β divergence, which involves
the Itakura-Saito divergence as a special case, and use themajorization-minimization (MM) approach to derive an iterative algorithm for finding the optimal parameters that minimize
the β divergence between Yj,k,l and Xj,k,l .

We follow the same idea proposed in [17], [18] to construct
an auxiliary function by bounding the convex term of theobjective function from above using Jensen’s inequality andbounding the concave term from above using a tangent line.
Owing to space limitations, here we only derive an update rule
for the magnitude part of the steering vector, namely Aj,i,k,n .
We can write the objective function as
)β
∑ β
1 ∑(∑
1
Yj,k,l +Aj,i,k,n Pi,k,l−nJ(Θ) =β(β − 1)
β
i,n
j,k,l
j,k,l
(
)β−1
∑
1 ∑
−
Yj,k,lAj,i,k,n Pi,k,l−n
,
(10)
β−1
i,n
j,k,l

where Θ denotes the set of the parameters. Here, both the
1 αsecond∑ and third terms of Eq. (10) involve the form α x1 with
x = i,n Aj,i,k,n Pi,k,l−n and α = β, β − 1. Note that α xα isconcave in x on the interval [0, ∞) when α ≥ 1 and convex
when α < 1. When α1 xα is convex, we can apply Jensen’s
inequality to build an upper bound function
(
)α
1 α
1∑
Aj,i,k,n Pi,k,l−n
,
(11)
x ≤
λj,k,l,i,nαα i,nλj,k,l,i,n∑where λj,k,l,i,n ≥ 0 and i,n λj,k,l,i,n∑ = 1. Equality holds
when λj,k,l,i,n = Aj,i,k,n Pi,k,l−n / i,n′ Aj,i,k,n′ Pi,k,l−n′ .
Next, when α1 xα is concave, we can use the fact that a concave
function lies below its tangent line to build an upper boundfunction
(∑
) Zα
1 α
j,k,l
α−1
x ≤ Zj,k,l
Aj,i,k,n Pi,k,l−n − Zj,k,l +, (12)ααi,nwhere x = Zj,k,l denotes
the point of tangency. Equality∑
holds when Zj,k,l = i,n Aj,i,k,n Pi,k,l−n . For simplicity ofnotation, we denote the right-hand sides of Eqs. (11) and(α)(α)(12) by Qj,k,l (Θ, ϑ) and Rj,k,l (Θ, ϑ), respectively, where ϑdenotes the set of auxiliary variables {λj,k,l,i,n } and {Zj,k,l }.
(α)(α)
By using Rj,k,l (Θ, ϑ) and Qj,k,l (Θ, ϑ), we arrive at thefollowing auxiliary function
∑ (β)
∑ β
1
J + (Θ, ϑ) =Yj,k,l +Sj,k,l (Θ, ϑ), (13)β(β − 1)
j,k,l

j,k,l

where

 (β)
(β−1)

Rj,k,l (Θ, ϑ) − Yj,k,l Qj,k,l (Θ, ϑ) β < 1
(β)
(β)
(β−1)
Sj,k,l = Qj,k,l(Θ, ϑ) − Yj,k,l Qj,k,l (Θ, ϑ) 1 ≤ β ≤ 2 .

 (β)
(β−1)
Qj,k,l (Θ, ϑ) − Yj,k,l Rj,k,l (Θ, ϑ) β > 2(14)According to the principle of the auxiliary function method,
we can show that iteratively minimizing this function withrespect to Θ and ϑ does not increase Eq. (10). The update
equation for Aj,i,k,n can thus be obtained as
)φ(β)
(∑
β−2
l Yj,k,l Xj,k,l Pi,k,l−n, (15)
Aj,i,k,n ← Aj,i,k,n∑ β−1
l Xj,k,l Pi,k,l−n1650

2016 24th European Signal Processing Conference (EUSIPCO)

6

SDR improvement [dB]

[m] 5

4

S2
3

S3

M2

M1

S1

2
1

Height: 4 m
0

1

2

4
3
2
0.75

0
3

4

5

6

7

8

9

10 [m](16)

Wi,k,m,τ ← Wi,k,m,τ
(∑
)φ(β)
β−2
j,l,n Yj,k,l Xj,k,l Aj,i,k,l−n Hi,m,n−τ,∑
β−1j,l,n Xj,k,l Aj,i,k,l−n Hi,m,n−τβ−2j,k,l,n Yj,k,l Xj,k,l Aj,i,k,l−n Wi,k,m,n−τ∑β−1j,k,l,n Xj,k,l Aj,i,k,l−n Wi,k,m,n−τ

0.85

0.9

0.95

Fig. 2. SDR improvements of proposed method and MNMF with different
settings of the reflection coefficient of walls

Similarly, the update rules for other parameters are given as

Hi,m,τ ← Hi,m,τ
(∑

0.8

Reflection coefficient of walls

Fig. 1. Room configuration for numerical simulation

where φ(β) is defined as
1/(2 − β) (β < 1)
φ(β) = 1(1 ≤ β ≤ 2) .


1/(β − 1) (2 < β)Proposed (T=4)
MNMF

5

)φ(β)
.

IV. E XPERIMENT
We conducted experiments on semi-blind source separationin reverberant environments to confirm the effectiveness of our
proposed method. We used 45 utterances extracted from the
ATR Japanese speech database [19] as the source signals. The
number of speakers was 3, and the number of utterances per
speaker was 15. First, to verify the robustness of the proposed
method against reverberation, the proposed method was tested
on synthetic signals simulating signals recorded in a room
where the walls had different reflection coefficients. Second,
to verify the robustness against the change in the positions of
the microphones, the proposed method was tested on synthetic
signals simulating signals recorded by moving microphones.
There were three sources and two microphones. The impulse responses were generated by using the image method.
Figure 1 shows the two-dimensional configuration of theroom. “•” (M1, M2) represent the microphone locations,and “×” (S1∼S3) represent the source locations. The signalsof sources 1 and 2 were female utterances, and those ofsource 3 were male utterances. We compared the proposed
method (Proposed) with the multichannel nonnegative matrix
factorization (MNMF) proposed in [3]. Since MNMF is basedon an anechoic (instantaneous) mixture model, its performanceis expected to degrade as the reverberation time of the test data
increases.
Out of the 15 utterances, one utterance was used for the
separation and the remaining 14 utterances were used for pretraining the spectrogram templates of each source. We iterated
the separation task 15 times for each selected utterance. For
pre-training, convolutive NMF was used for the proposed
method whereas NMF was used for MNMF. The frame number of each spectrogram template for both convolutive NMF
and the proposed method was set experimentally at T = 4.
The numbers of the templates of NMF and convolutive NMF
were 40 and 20, respectively. The generalized KL divergence
was adopted as a divergence measure for pre-training with
NMF and convolutive NMF (where the parameter of the βdivergence was β = 1). The STFT window length was 32 ms,
and the shift length was 16 ms. We evaluated the source
separation performance using the improvement of the signalto-distortion ratio (SDR).
Figure 2 shows the SDR improvement results with thedifferent reflection coefficients of the walls obtained withthe proposed method and MNMF. The SDR improvements
were averaged over the SDR improvements of each sourceon each microphone. For instance, the reverberation times(RT60 ) were 620 ms, 720 ms and 1980 ms when the reflection
coefficient of the walls were 0.75, 0.8 and 0.95, respectively.
When the reflection coefficient of the walls was low, MNMF
outperformed the proposed method. On the other hand, theproposed method outperformed MNMF as the reverberation
time increased.
Figure 3 shows the SDR improvement results with different
T settings. Note that when T = 1, the source spectrogram
model reduces to the regular NMF model. As shown by
the result, T = 4 outperformed T = 1, indicating that the
convolutive NMF model is more suitable for expressing speech
spectrograms than the regular NMF model.
Second, the proposed method was tested on synthetic signalsthat simulated signals recorded by moving microphones. The
test signals were created by concatenating the signals recordedby the microphones located as shown in Figure 1 and by those
moved a distance of ∆x m along the x-axis. The reflection
coefficient of the walls was set at 0.8. Figure 4 shows the1651

2016 24th European Signal Processing Conference (EUSIPCO)R EFERENCES

SDR improvement [dB]6
Proposed (T=4)
Proposed (T=1)

5
4
3
2
0.75

0.8

0.85

0.9

0.95

Reflection coefficient of walls

Fig. 3. SDR improvements of proposed method with different T settings

SDR improvement results with different ∆x settings. The
degradation of the MNMF performance, which assumes thetime-invariance of the mixing system, was more significant
than that of the proposed method as ∆x increased. This
may imply that the phase part of the steering vector is moresensitive to the microphone movements than the magnitude
part.
5
Proposed
MNMF

SDR improvement [dB]

4.9
4.8
4.7
4.6
4.5
4.4
4.3
0

0.04

0.08

0.12

0.16

0.2

∆ x [m]Fig. 4. SDR improvements of proposed method and MNMF with different∆x settings

V. C ONCLUSION
This paper proposed an extension of the multichannel nonnegative matrix factorization (NMF) approach to deal with the
problem of underdetermined source separation in time-variant
reverberant environments tailored for an ad hoc microphonearray. The proposed method uses a non-negative convolutive
mixture model to take account of the reverberation in a room
as well as to allow for relatively small changes of the microphone/source positions. The method also uses the convolutive
NMF model to express the power spectrogram of speech.
We verified through experiments that the proposed method isrobust against reverberation and microphone movements.

[1] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix
factorization in convolutive mixtures for audio source separation,” IEEE TASLP, 18(3), 550–563, 2010.
[2] A. Ozerov et al., “Multichannel nonnegative tensor factorization
with structured constraints for user-guided audio source separation,” in Proc. ICASSP, 257–260, 2011.
[3] H. Sawada et al., “Multichannel extensions of non-negative
matrix factorization with complex-valued data,” IEEE TASLP,
21(5), 971–982, 2013.
[4] T. Higuchi and H. Kameoka, “Joint audio source separation anddereverberation based on multichannel factorial hidden Markov
model,” in Proc. MLSP, 1–6, 2014.
[5] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” in Proc. WASPAA,
177–180, 2003.
[6] H. Kameoka et al., “Statistical model of speech signals basedon composite autoregressive system with application to blindsource separation,” in Proc. LVA/ICA, 245–253, 2010.
[7] D. Kitamura et al., “Efficient multichannel nonnegative matrix
factorization exploiting rank-1 spatial model,” in Proc. ICASSP,
276–280, 2015.
[8] T. Nakatani et al., “Blind speech dereverberation with multichannel linear prediction based on short time Fourier transform
representation,” in Proc. ICASSP, 85–88, 2008.
[9] T. Yoshioka et al., “Blind separation and dereverberation of
speech mixtures by joint optimization,” IEEE TASLP, 19(1), 69–
84, 2011.
[10] S. Miyabe et al., “Blind compensation of inter-channel sampling
frequency mismatch with maximum likelihood estimation in
STFT domain,” in Proc. ICASSP, 674–678, 2013.
[11] H. Kameoka et al., “Robust speech dereverberation based on
non-negativity and sparse nature of speech spectrograms,” inProc. ICASSP, 45–48, 2009.
[12] H. Chiba et al., “Amplitude-based speech enhancement with
nonnegative matrix factorization for asynchronous distributed
recording,” in Proc. IWAENC, 203–207, 2014.
[13] P. Smaragdis, “Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs,” inProc. ICA, 494–499, 2004.
[14] P. D. O’Grady and B. A. Perlmutter, “Convolutive non-negative
matrix factorisation with a sparseness constraint,” in Proc.
MLSP, 427–432, 2006.
[15] C. Févotte, N. Bertin, and J. L. Durrieu, “Non-negative matrix
factorization with the Itakura-Saito divergence: With application
to music analysis,” Neural Computation, 21(3), 793–830, 2009.
[16] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix
factorization,” in Adv. NIPS, 556–562, 2001.
[17] H. Kameoka et al., “Selective amplifier of periodic and nonperiodic components in concurrent audio signals with spectral
control envelopes,” in IPSJ SIG Technical Reports, 2006-MUS66-13,77–84, in Japanese, Aug. 2006.
[18] M. Nakano et al., “Convergence-guaranteed multiplicativealgorithms for non-negative matrix factorization with betadivergence,” in Proc. MLSP, 283–288, 2010.
[19] A. Kurematsu, et al., “ATR Japanese speech database as a toolof speech recognition and synthesis,” Speech Communication,
9(4), 357–363, 1990.

ACKNOWLEDGMENT
This work was supported by JSPS KAKENHI Grant Numbers 26730100 and 26280060.

1652