1

Multichannel Linear Prediction for Blind
Reverberant Audio Source Separation

arXiv:1702.07713v1 [cs.SD] 24 Feb 2017

İlker Bayram and Savaşkan Bulek
Abstract—A class of methods based on multichannel linear
prediction (MCLP) can achieve effective blind dereverberation of
a source, when the source is observed with a microphone array.
We propose an inventive use of MCLP as a pre-processing step
for blind source separation with a microphone array. We show
theoretically that, under certain assumptions, such pre-processing
reduces the original blind reverberant source separation problem
to a non-reverberant one, which in turn can be effectively tackled
using existing methods. We demonstrate our claims using real
recordings obtained with an eight-microphone circular array in
reverberant environments.
Index Terms—Microphone array, source separation, blind
dereverberation, multichannel linear prediction, beamforming,
post-filtering.

I. I NTRODUCTION
In anechoic environments, audio source separation with a
microphone array can be achieved by beamforming and postfiltering [7]. On the other hand, in reverberant environments,
reverberation causes significant degradation in the quality and
intelligibility of the beamformer output, especially if the room
impulse responses (RIR) are not taken into account or are
not known precisely. This is a serious limitation because
in practice, RIRs are usually not known with precision and
they can vary wildly with respect to the positions of the
microphone/source pair.
A related problem addressing reverberation without explicit
RIR knowledge is called blind dereverberation [17]. A successful class of methods for blind dereverberation fall under
the framework of multi-channel linear prediction (MCLP)
[16], [11], [10]. These methods require a microphone array
and employ linear prediction to extract information about
the properties of reverberation. In this paper, we propose
to employ MCLP, which is primarily used for single source
dereverberation, as a pre-processing step in source separation.
We show theoretically and demonstrate numerically that such
a scheme not only dereverberates the sources but also helps
in the subsequent separation.
In order to make this statement more precise, let us briefly
discuss an MCLP formulation. Consider an array of microphones, recording a source in a reverberant room. Recent
MCLP formulations [16], [11], [10] select one of the microphones as the reference and use delayed versions of the whole
set of observations to linearly predict the reference signal.
The residual of prediction, which may also be interpreted
as the ‘innovations’ component, forms the dereverberated
İ. Bayram is with the Dept. of Electronics and Communications Eng.,
Istanbul Technical University, Istanbul, Turkey. E-mail : ibayram@itu.edu.tr.
S. Bulek is with Qualcomm Atheros, Inc., Auburn Hills, MI, USA. E-mail :
sbulek@gmail.com.

estimate of the source. In our scenario, we have multiple
sources recorded with a microphone array and we would like
to separate and dereverberate the sources. For this, suppose
we apply MCLP for each possible choice of the reference
microphone. If we have N microphones in the array, this
results in N hopefully dereverberant but mixed observations.
We say ‘hopefully dereverberant’ because the use of MCLP
has previously been considered and justified solely for the case
of a single source [11]. Further, even if the MCLP outputs thus
obtained are dereverberated, in order to achieve separation,
we have to know the relation between the recordings. In
this paper, we show that, under a mild assumption on the
RIRs, MCLP actually suppresses reflections other than the
direct source and preserves the phase of the direct source.
Therefore, the MCLP preprocessing step essentially reduces
the reverberant separation problem to an anechoic one. This
greatly simplifies the separation step because we can now use
the known array geometry to estimate the direction of arrival
and perform separation using this estimate. We argue that this
effect is not dependent on the following DOA estimation and
separation algorithms. In that sense, the proposed MCLP stage
complements currently available source separation algorithms,
improving their performance. Consequently, we do not have
any reservations about which method to use following the
MCLP stage. Nevertheless, for the sake of completeness, we
include a brief description of the geometric source separation
(GSS) algorithm [18] followed by a post-filter.
Related Work
Literature on dereverberation/source separation with a microphone array is vast. Here, we focus on blind 1 methods for
dereverberation and source separation for arrays with a known
geometry.
For source separation, provided the directions of arrival
(DOA) of the sources are known, and reverberation is limited, beamformers like the minimum variance distortionless
response (MVDR) beamformer, along with a post-filter [2],
[7], [23] can be effective. In MVDR, one aims to minimize
the energy of the reconstructed signals, while preserving
the signals from the DOAs. However, reverberation causes
a source to contribute to different directions, reducing the
effectiveness of such approaches. Geometric source separation
(GSS) [18] addresses this problem by also enforcing statistical
independence of the reconstructed sources. This is achieved
by selecting the beamformer weights so as to minimize a cost
1 We use the term ‘blind’ to refer to methods that do not require to know
the RIRs. However, the blind methods we consider may use the known array
geometry, which is independent of the characteristics of the environment.

2

function that consists of the sum of a term imposing geometric
constraints based solely on the array geometry and a term that
promotes the statistical independence of the sources (motivated
by blind source separation techniques [9]). In practice, GSS
outputs are still reverberant for moderate to highly reverberant
rooms. GSS has been complemented with a post-filter in [28],
[27] to devise a high performance source-separation method
that can operate in real time. An alternative approach, based
on the assumption that the activity of the different sources
do not overlap in the time-frequency domain [31], has been
proposed in [24]. However, the model in [24] assumes that the
room impulse responses are relatively short, compared to the
analysis window size used in the STFT.
Despite its simplicity, blind dereverberation is a challenging
problem due mainly to the characteristics of room impulse
responses. For a general overview and discussion of relatively
recent methods, we refer to the edited book [17]. Here, we
focus on a successful class of methods that can be collected
under the heading ‘multichannel linear prediction’ (MCLP), as
they are more directly related to the content of our work. These
methods usually fit in a schema as follows. Given multiple
observations yi , for i = 1, 2, . . . , M , one of the observations
is chosen as the reference. Then, this reference observation is
linearly predicted by the delayed version of all of the observations. Once the linear prediction coefficients are determined,
this information has been used differently in prior work. For
instance, in [12], the authors use the coefficients to estimate the
power of late reverberations, which are subsequently reduced
by spectral subtraction. In more recent work [16], [11], [10],
the residual of linear prediction is taken as the dereverbed
audio signal. Prior work also differentiates based on the criteria
used for determining the weights used in linear estimation. In
[12], the authors propose to select the linear estimation weights
so as to minimize the energy of the residual. In [16], the source
signal is modeled as a time-varying Gaussian process and the
weights are selected so as to maximize the likelihood of such
a model. More recently, [11], [10] employ a sparsity- based
model in the time-frequency domain and select the weights so
that the residual minimizes a sparsity promoting cost function.
An interesting approach that employs MCLP to achieve
suppression of late reverberation and source separation is
presented in [26]. The authors use MCLP to suppress late
reverberation, using a statistical framework similar to that in
[16]. In contrast to [16], a time-varying scenario with multiple
sources is considered and the MCLP weights are determined
along with the statistics of the late reverberant part of the
sources in an iterative manner.
Contribution
Our proposal in this paper is of a complementary nature and
is likely to improve the performance of existing separation
algorithms. We present an inventive way to employ MCLP
as a pre-processing step for blind source separation using a
microphone array. We show theoretically that, under certain assumptions, such a pre-processing step converts the reverberant
source separation problem into a non-reverberant one, which in
turn can be solved by existing methods such as [24], [18], [27]

outlined above. Such an approach is inherently different than
previous uses of MCLP which aim to dereverb a single source
as in [12], [16], [11], [10]. Our proposal also differs from that
of [26] because [26] requires the source statistics to cancel the
late reverberation and thus resorts to a joint estimation of the
linear prediction weights and source statistics. This approach
in [26] requires an iterative procedure and is handled with an
EM algorithm. In contrast, we propose to employ the MCLP
step once, in order to reduce the reverberant problem to an
approximately non-reverberant one.
In addition to theoretical justification, we demonstrate the
validity of the claims using real recordings. For the sake of
completeness, we include a description of GSS along with
a simple post-filter, to be applied following the MCLP step,
in order to achieve separation. However, after MCLP, other
methods that assume short RIRs, or operate under a free
propagation model may also be used, as noted above.
Notation
Sequences and continuous functions of time are denoted
with a small letter as in h(n) or h(t), respectively. The DFT
and STFT of sequences are denoted with a capital letter as
H(ω) and H(s, ω) respectively – which one is referred to
will be clear from the context. The z-transform of sequences
are denoted with calligraphic letters as in H(z). We refer to
H(z) as ‘causal’ if it is a polynomial of z −1 . For a causal
H(z), we denote the constant h(0) also as H(∞).
Throughout the paper, K denotes the number of sources
and N denotes the number of microphones in the array.
II. D ESCRIPTION OF THE M ODEL
Throughout the paper, we are interested in far-field processing. Thus, the distance of the source to the microphones
is assumed to be significantly greater than the array size.
Our proposed formulation/algorithm is in the STFT domain.
However, for simplicity of exposition, we first introduce a
model in the continuous-time domain and then present an
approximation of this model in the STFT domain.
A. A Time-Domain Observation Model
We introduce the time-domain model in three stages, where
complexity increases with each stage. We start with the case
of a single source, single microphone. In the second stage, we
consider the single source, multiple microphone case. Finally,
we discuss the multiple source, multiple microphone case,
which is the main scenario of interest in this paper.
1) Single Source, Single Microphone: For a given
microphone-source pair, let x(t) denote the signal produced
by the source. The observation y(t) can be approximated as,
y(t) ≈ x(t) ∗ h̃(t),

(1)

where h̃(t) is the room impulse response (RIR) for the given
source-microphone positions. Now suppose that the distance
from the source to the microphone is D meters and the speed
of sound in the environment is c m/sec. Then, the acoustic
signal emitted by the source reaches the microphone after a

3
3

delay of D/c seconds. We declare the distance compensated
room impulse
response We
(DCIR)
to be
delay
of D/c seconds.
declare
the distance compensated
(DCIR) to be
room h(t)
impulse
response
= h̃ t + D/c .
(2)

h(t)
=for
h̃ tpractical
+ D/c .purposes, we can work with the DCIR
(2)
Notice
that
insteadthat
of the
for apurposes,
single microphone-source
For
Notice
for RIR
practical
we can work with setup.
the DCIR
multiple
microphones
or
sources,
we
need
to
take
into
account
instead of the RIR for a single microphone-source setup. For
relative distances,
as discussed
multiple
microphones
or sources,next.
we need to take into account
2)
Single
Source,
Multiple
Consider a setup
relative distances, as discussedMicrophones:
next.
consisting
of
a
microphone
array
and
a
source
as shown
in
2) Single Source, Multiple Microphones: Consider
a setup
Fig.
1.
Suppose
we
set
the
center
of
the
array
as
the
origin
for
consisting of a microphone array and a source as shown in
a coordinate
system.
this
origin,
letarray
the vector
pi denote
Fig.
1. Suppose
we setFor
the
center
of the
as the origin
for
th
the
position
of
the
the
i
microphone.
Also,
let
u
denote
the
x
a coordinate system. For this origin, let the vector pi denote
th
unit
vector
in
the
direction
of
the
source.
Finally,
let
d
=
i
the position of the the i microphone. Also, let ux denote the
hp
,
u
i.
Note
that
d
denotes
the
coordinate
of
the
projection
i
x
i
unit vector
in the direction of the source. Finally, let di =
th
ofi ,the
microphone’s
position
the subspace
by
hp
ux i.i Note
that di denotes
the to
coordinate
of thespanned
projection
thinduced time-delay is τi = −di /c. Notice that τi may
u
.
The
x
of the i microphone’s position to the subspace spanned by
or negative,
i.e., may
advance
ubex .positive
The induced
time-delay
is ⌧i represent
= di /c.anNotice
thator⌧idelay.
may
th
microphone
Finally,
let
the
DCIR
from
the
source
to
the
i
be positive or negative, i.e., may represent an advance or delay.
be denoted
as hDCIR
we source
model the
observation
at the
i (t). Then,
Finally,
let the
from the
to the
ith microphone
ibeth denoted
microphone,
namely
y
,
as,
i
as hi (t). Then, we model the observation at the
ith microphone,
yi ,τas,
yi (t) ≈ x(t)namely
∗ hi (t −
i ) = x(t − τi ) ∗ hi (t).

(3)

⇡ x(t)
⇤ hi (t ⌧amounts
hi (t).
i (t)field
i ) = x(t to ⌧the
i ) ⇤ following
The yfar
assumption
: If (3)
the
distance
of theassumption
source to the
center of
the
array
is D, then
the
The
far field
amounts
to
the
following
:
If
the
distance of
of the
the source
source to
to the
thecenter
ith microphone
is approximately
distance
of the array
is D, then the
D − di . of
Under
the far-field
the isdescribed
model
distance
the source
to the assumption,
ith microphone
approximately
preserves
the ability
of the array
to differentiate
differentmodel
direcD
di . Under
the far-field
assumption,
the described
tions. We the
nextability
add another
layertotodifferentiate
this model different
by considering
preserves
of the array
direcmultiple
tions.
Wesources.
next add another layer to this model by considering
3) Multiple
Sources, Multiple Microphones: Suppose now
multiple
sources.
that
are Sources,
K sources
and NMicrophones:
microphones.Suppose
Let thenow
k th
3) there
Multiple
Multiple
source
signal
(t).microphones.
Also, let the Let
unitthe
vector
that
there
are be
K denoted
sourcesth by
andxkN
k th
in the direction
the k by
source
denoted
by uunit
the
k . For
source
signal be ofdenoted
xk (t).beAlso,
let the
vector
ith the
microphone,
positioned
at pi , be
wedenoted
set τi,kby=uhp
,
u
i/c.
i
k
in
direction of
the k th source
.
For
the
k th
th
let hi,k (t)
denote the
DCIR
between
the
k i , usource
iFinally,
microphone,
positioned
at
p
,
we
set
⌧
=
hp
i
i,k
k i/c.
th
th
and
the
i
microphone.
Then,
we
model
the
observation
at
Finally,
let
h
(t)
denote
the
DCIR
between
the
k
source
i,k
th
th
the
i
microphone
as,
and the i microphone. Then, we model the observation at
K
K
the ith microphone
as,
X
X
yi (t) ≈ K xk (t)∗hi,k (t−τi,k ) = K xk (t−τi,k )∗hi,k (t) (4)
X
X
yi (t) ⇡ k=1 xk (t)⇤hi,k (t ⌧i,k ) = k=1 xk (t ⌧i,k )⇤hi,k (t) (4)
k=1 the absolute delay from
k=1 the k th source is not
Notice that,
th
included
in
the
model.
However,
this
is not
omission,
Notice that, the absolute delay from
thea kcareless
source
is not
because
this
delay
is
in
fact
not
known
and
is
not
estimated
included in the model. However, this is not a careless omission,
in practice.
for
because
this Therefore,
delay is in the
fact described
not knownmodel
and isisnotsuitable
estimated
practical
far-field
processing
scenarios.
in practice. Therefore, the described model is suitable for

practical far-field processing scenarios.
B. An Approximate Model in the STFT Domain

B. We
An Approximate
the STFT
Domain to be briefly
work in the Model
STFT in
domain
for reasons
clarified
in
the
sequel.
We
are
specifically
interested
in the
We work in the STFT domain for reasons
to be briefly
representation
of
an
LTI
system
in
the
STFT
domain.
Consider
clarified in the sequel. We are specifically interested in the
a convolution of
equation
representation
an LTI as,
system in the STFT domain. Consider
a convolution
equation
as,
y(n) = x(n)
∗ h(n).

(5)

= x(n)
⇤ h(n).
Let y(n)
Y (s, ω),
X(s,
ω) denote the STFTs of y and x (5)
respectively,
obtained
using
a
window
w(n).
We
differentiate
Let Y (s, !), X(s, !) denote the STFTs of y and x respectively, obtained using a window w(n). We differentiate

m2
x
m3 |d3 |

|d1 | m1
m4

Fig. 1. Microphones and the source are denoted by mi , and x respectively.
Fig. center
1. Microphones
andmarked
the source
denoted
by mto
x respectively.
i , and
The
of the array,
with are
a circle,
is taken
be the
origin. The
The center
of the
with is
a circle,
is taken
to be
the spanned
origin. The
position
vector
for array,
the iththmarked
microphone
projected
onto the
space
by
position
vector for
the i inmicrophone
is projected
onto to
theobtain
space the
spanned
by
the
unit vector
pointing
the direction
of the source
relative
the unit vector
pointing
of are
the positive
source towhereas
obtain dthe
relative
coordinate
di . Note
that din1 , the
d2 direction
(not shown)
3 , d4 (not
coordinate
i . Note that d1 , d2 (not shown) are positive whereas d3 , d4 (not
shown)
are dnegative.
shown) are negative.

between
between two
two cases,
cases, following
following the
the discussion
discussion in
in [13]
[13] (see
(see
Lemma
1
in
[13]).
Lemma 1 in [13]).
(I)
(I) If
If the
the support
support of
of hh is
is contained
contained in
in the
the support
support of
of w,
w, then
then

Y (s, !) ⇡ X(s, !) H(!),
(6)
Y (s, ω) ≈ X(s, ω) H(ω),
(6)
where H(!) denotes the DFT of h.
where
H(ω) denotes
the DFT oflarger
h. than the support
(II) If
the support
of h is significantly
(II) of
If the
support
of
h
is
significantly
larger
than the support
w, then
of w, then
s

s
YY (s,
(s, !)
ω) ⇡
≈ X(s,
X(s, !)
ω) ⇤∗ H̃(s,
H̃(s, !),
ω),

(7)
(7)

ss

for
for some
some H̃(s,
H̃(s, !),
ω), where
where ⇤∗ denotes
denotes convolution
convolution with
with
respect
to
the
time
variable
s,
defined
respect to the time variable s, defined as
as
X
X
s
X(s, !)
ω) ⇤s∗ H̃(s,
H̃(s, !)
ω) =
=
X(s − l,l, !)
ω) H̃(l,
H̃(l, !).
ω). (8)
(8)
X(s,
X(s
ll

We remark that H̃(s, ω) is not necessarily the STFT of
We remark that H̃(s, !) is not necessarily the STFT of
h, but is chosen so as to make the approximation as
h, but is chosen so as to make the approximation as
close as possible – see [20] for a discussion.
close as possible – see [20] for a discussion.
We remark that, as a special case of (I), if τ is small
We remark
that,
as aofspecial
if ⌧ is small
enough,
then the
STFT
x(n − τcase
) canofbe(I),approximated
as
enough,
then
the
STFT
of
x(n
⌧
)
can
be
approximated
as
X(s, ω) exp(−jω τ ). Therefore, the model introduced in (4)
X(s,
!)approximated
exp( j! ⌧ ). in
Therefore,
modelasintroduced in (4)
can be
the STFTthe
domain
can be approximated in the STFT domain as
K


X
s
K Xk (s, ω) ∗ ⇣ H̃i,k (s, ω) exp(−jω τi,k )⌘ . (9)
Yi (s, ω) ≈ X
s
Yi (s, !) ⇡ k=1 Xk (s, !) ⇤ H̃i,k (s, !) exp( j! ⌧i,k ) . (9)
k=1

Such a model is preferred over the time domain model
because
or is
penalty
functions
that time
enforce
(structured)
Such aprior
model
preferred
over the
domain
model
sparsity
of
the
audio
signals
are
typically
defined
in
the STFT
because prior or penalty functions that enforce (structured)
(or someofother
time-frequency)
domain. defined
If we use
theSTFT
time
sparsity
the audio
signals are typically
in the
domain
convolutional
model
(4)
in
an
iterative
algorithm,
(or some other time-frequency) domain. If we use the time
we haveconvolutional
to work in themodel
time (4)
domain
to iterative
undo the algorithm,
effects of
domain
in an
reverberation,
and
the
STFT
domain
to
take
advantage
of the
we have to work in the time domain to undo the effects
of
signal
prior.
Switching
from
the
time
domain
representation
reverberation, and the STFT domain to take advantage of the
to the prior.
STFTSwitching
representation
at each
increases the
signal
from the
time iteration
domain representation
computational
cost
significantly.
Working
with
a
convolutional
to the STFT representation at each iteration increases the
model directlycost
in the
STFT domain,
such
allows to
computational
significantly.
Working
withasa (9),
convolutional
carry
out
the
iterations
in
the
STFT
domain,
eliminating
model directly in the STFT domain, such as (9), allows the
to
need to
the forward
STFT operations
at each
carry
outapply
the iterations
in and
the inverse
STFT domain,
eliminating
the
iteration.
need
to apply the forward and inverse STFT operations at each
iteration.

4

III. M ULTI -C HANNEL L INEAR P REDICTION FOR
D EREVERBERATION
In this section, we discuss the implications of using multichannel linear prediction (MCLP) as a preliminary step for the
source separation problem. We derive conditions under which
MCLP transforms the original blind separation/dereverberation
problem to a non-blind one.
We start the discussion with a single source and then
consider extensions to multiple sources.
A. MCLP for a Single Source
as

Consider a signal x that is observed via different channels
yi = x ∗ hi , for i = 1, 2, . . . , N,

(10)
th

where hi denotes the impulse response of the i channel. In
the z-domain, this can be written as,

 

Y1 (z)
H1 (z)
 ..   .. 
(11)
 .  =  .  X (z).
YN (z)

HN (z)

We would like to find causal filters G1 (z), . . . , GN (z) such
that


H1 (z)



G1 (z) . . . GN (z)  ...  = 1.
(12)
HN (z)

MINT [15] ensures that such filters exist provided that Hi ’s
are co-prime, i.e., they do not have a common non-constant
factor. Even though it is relieving to know that such filters
exist, finding the filters is still a challenging problem when
hi ’s are not known. MCLP addresses this challenge.
In MCLP, the idea is to estimate x by subtracting from yr
(for a chosen reference index r) delayed weighted combinations of yi ’s. This is equivalent to using in (12) a set of filters
of the form
Gr (z) = 1 − z −d Ur (z) for a reference index r,
Gi (z) = −z −d Ui (z) for i = 1, . . . , N, i 6= r,

(13)

for a positive integer d (determining the amount of delay) and
causal Uk ’s. In [11], the selection of the filters is guided by
the assumption that the original x is sparse. More specifically,
suppose yi (s) denotes the 1D time-series derived from a
frequency band of the observations, that is yi (s) = Yi (s, ω)
for some ω. The filters ui are found by solving the following
minimization problem.
!
N X
L
X
min P yr (s) −
yi (s − d − l) ui (l) , (14)
u1 ,...,uN

i=1 l=0

where P (·) is a sparsity promoting function. We will have
more to say about P , but first we would like to stress a
point about the choice (13). MINT ensures the existence of
dereverberating filters but does not guarantee that they can
be of the form (13). Below, we propose modified conditions
which ensure the existence of such Gi ’s. Another issue is the
role of the reference index r. It turns out that repeating the

linear prediction for different r values leads to an interesting
outcome, summarized in Prop. 2 below. In order to derive
Prop. 2, we need an auxiliary result (Prop. 1) and an assumption.
Proposition 1. Suppose that the N − 1 polynomials Hi , i ∈
{1, . . . , N } \ {r}, are co-prime. Then, for d = 1, there exist
a set of filters of the form (13) such that


H1 (z)

 . 
G1 (z) . . . GN (z)  ..  = c,
(15)
HN (z)

for a constant c. Further, for a set of causal filters of the form
(13), if (15) holds, then c = Hr (∞).
Proof. See Appendix A.

Consider now the application of this result to the model (9).
Assume that there is a single source, i.e., K = 1. For a fixed
frequency ω, let
(16)

hi (s) = H̃i,1 (s, ω) exp(−jω τi,1 ).

Under the hypotheses of Prop. 1, if, for d = 1, MCLP can
recover a set of filters satisfying (15), then Prop. 1 implies
that,
(17)

ŷr (s) = hr (0) exp(−jω τr,1 ) x(s).
Now, if further,
hi (0) = hi0 (0) for any distinct pair (i, i0 ),

(18)

then,
ŷr (s) = c0 exp(−jω τr,1 ) x(s),

(19)

for some constant c , independent of r (actually, c = hi (0),
for any i). We remark that the difference of (17) and (19) is
that the rhs of (19) depends only on the direction of arrival,
whereas this is not the case in (17). In other words, MCLP
eliminates reflections from different directions.
The assumption in (18) states that, around the origin,
distance compensated impulse responses are similar. This
assumption is key to the development in this paper and we state
it explicitly for later reference. We will verify this assumption
experimentally in Section V.
0

0

Assumption 1. For a given source - microphone array pair,
let H̃i (s, ω) denote the STFT representation of the distance
compensated impulse response for the ith microphone. Then,
H̃i (0, ω) = H̃i0 (0, ω),

(20)

for i 6= i and all ω.
0

We state the foregoing observations as a proposition.

Proposition 2. For a fixed frequency ω, let hi (s) =
H̃i,1 (s, ω) exp(−jω τi,1 ). Assume also that the collection of
polynomials Hi , i ∈ {1, . . . , N } \ {r} are co-prime. Then, for
d = 1, there exist a set of filters of the form (13) such that
X
gi (n) ∗ hi (n) = cr δ(n)
(21)
i

5

for a constant cr .
Further, for a set of filters of the form (13), if (21) holds,
and Assumption 1 is in effect, then cr = c0 exp(−jωτr,1 ) for
a constant c0 , independent of r.
Next, we extend the discussion to the case of multiple
sources.

Consider an observation model with K sources x1 ,. . . xK
and N > K observations defined as,
K
X

k=1

xk ∗ hi,k , for i = 1, 2, . . . , N.

In the z-domain, this can be

 
H1,1 (z)
Y1 (z)
 ..  
 . =
YN (z)

(22)

written as,


. . . H1,K (z)  X1 (z) 

..
. . .  (23)

.
HN,1 (z) . . . HN,K (z) XK (z)
{z
}
|
H

For this model, MINT also ensures that we can find filters
Gk,i , for k = 1, . . . , K, i = 1, . . . , N such that



H1,1 (z) . . . H1,K (z)
G1,1 (z) . . . G1,N (z)



..
..
 = I,


.
.
GK,1 (z) . . .

GK,N (z)

HN,1 (z) . . .

HN,K (z)

(24)

provided that the (Smith) canonical form of H (see Appendix B or [8], Sec.VI.2, Defn.3) has a non-zero constant
(i.e., no polynomial higher than zero order) diagonal. However,
unlike the single source case, we cannot use filters of the form
(13) to achieve (24). We have the following extension in this
case.
Proposition 3. For a fixed reference microphone index r ∈
{1, . . . , N }, let H̄r be the matrix obtained by removing the rth
row of H. If the greatest common divisor of all the K-minors
of H̄r is a non-zero constant, then, for d = 1, there exist a
set of filters of the form (13) such that


H1,1 (z) . . . H1,K (z)



(r)
..
G1 (z) . . . GN (z) 
=c ,
.
HN,1 (z) . . .

HN,K (z)

(25)

for a constant vector c(r) (i.e., containing no polynomial of
order greater than zero).
Further, for a set of filters of the form (13), if (25) holds,
then the components of c satisfy,
(r)

ck = Hr,k (∞), for k = 1, 2, . . . , K.

When we apply this result to our observation model with
K sources, we obtain the following proposition.
Proposition 4. For a fixed frequency ω, let
hi,k (s) = H̃i,k (s, ω) exp(−jω τi,k ),

B. MCLP for Multiple Sources

yi =

is shown in Prop. 8 in Appendix B. Therefore, we can replace
the hypothesis in Prop. 3 with an alternative, stated in terms
of the Smith canonical form of H̄r .

(26)

Proof. See Appendix C.

Remark 1. For an N × K matrix H with K ≤ N , it can be
shown that, the greatest common divisor of all the K-minors
of H is a non-zero constant if and only if the Smith canonical
form of H has a non-zero constant diagonal. This equivalence

(27)

for i ∈ {1, . . . , N }, k ∈ {1, . . . , K}. Also, let H be defined as
in (23), and H̄r be defined as in Prop. 3. Finally, let yi (s) =
Yi (s, ω).
If the greatest common divisor of all the K-minors of H̄r
is a non-zero constant, then, for d = 1, there exist a set of
filters of the form (13) such that
N
X
i=1

(r)

gi (s) ∗ hi,k (s) = ck δ(s)

(28)

(r)

for some constants ck .
Further, for a set of filters of the form (13), if (28) holds,
and Assumption 1 is in effect, then we have
N
X
i=1

gi (s) ∗ yi (s) = c0

K
X

xk (s) exp(−j ω τr,k ),

(29)

k=1

where c0 is a constant independent of r.
Proof. See Appendix D.
This proposition suggests that, in a scenario with multiple sources, MCLP, if successful, suppresses reflections and
preserves the direct signal even when there exist multiple
sources from different directions. Therefore, after the MCLP
step, the recordings appear as if they are recorded with the
same microphone array, but in a non-reverberant environment.
This in turn allows to read off the directions of the sources
with a simple procedure like MUSIC [22]. Once the source
directions are determined, we can then achieve separation by
using known manifold vectors of the microphone array.
IV. I MPLEMENTATION D ETAILS
The discussion after Prop. 4 sets the stage for the proposed
framework. The main building blocks of the framework are
shown in Fig. 2. We first input the N observations y1 , . . . , yN ,
recorded in a reverberant environment with unknown reverberation, to an MCLP block. This block applies MCLP for every
possible choice of the reference microphone. The output of this
block is N dereverberated signals ŷ1 , . . . , ŷN . These signals
are input to a DOA estimation algorithm (such as MUSIC
[22]) and the directions of the K sources are determined. Then
these estimated directions are used to separate the sources in
ŷi using a non-blind source separation method. The estimated
sources are denoted as x̂i in Fig. 2. In the following, we briefly
present the algorithm we used for a practical implementation
of the components. However, we remark that following MCLP,
the problem turns into a non-reverberant source separation
problem (with known array geometry) and one can use any
method of choice for this task.

6
6

y1
yN

(i)

ŷ1
..
.

MCLP

..
.

ŷN

Source ..
Sep. .

···

(i)
where
uuk(i) =
ui (k). Also,
where uu(i) denotes
denotes the
the column
column vector
vector with
with(r)
k = ui (k). Also,
(r)
let
y
(r) denote the column vector with yk(r) = yr (k). Using
let y
denote the column vector with yk = yr (k). Using
these,
these, we
we can
can rewrite
rewrite the
the problem
problem (31)
(31) as,
as,

x̂1
x̂K

min
kzkw , subject to z = y rr − X
min
X uu
z,u kzkw , subject to z = y
z,u

Est.
DOA

Fig. 2.
2. The
The building
building blocks
blocks of
of the
the proposed
proposed framework.
framework. The
The input
input to
to
Fig.
the system
system is
is yy11,, ......yyN
which are
are N
N reverberant
reverberant observations
observations of
of K
K
N,, which
the
sources. Source
Source directions
directions are
are estimated
estimated under
under the
the assumption
assumption of
of known
known array
array
sources.
geometry. The
The estimated
estimated source
source directions
directions are
are used
used in
in aa non-blind
non-blind sourcesourcegeometry.
separation algorithm.
algorithm.
separation

A. Multi-Channel
Multi-Channel Linear
Linear Prediction
Prediction
A.
We
use
the
formulation
given in
in Sec.
Sec. III-B.
III-B. We
We work
work in
in the
the
We use the formulation given
STFT
domain
and
process
the
frequency
channels
separately.
STFT domain and process the frequency channels separately.
For aa fixed
fixed value
value of
of ω,
!, we
we set
set yyii(s)
(s) =
= YYii(s,
(s, ω).
!). Then
Then for
for each
each
For
value
of
r
2
{1,
.
.
.
,
N
},
the
dereverbed
estimate
ŷ
is
defined
value of r ∈ {1, . . . , N }, the dereverbed estimate ŷrr is defined
as,
as,

(36)
(36)

We
We split
split variables
variables and
and introduce
introduce zz in
in order
order to
to make
make the
the
formulation
formulation feasible
feasible for
for employing
employing DR.
DR.
For
For this
this formulation,
formulation, the
the Douglas-Rachford
Douglas-Rachford algorithm
algorithm is
is of
of
the
the form
form outlined
outlined in
in Algorithm
Algorithm 1.
1. We
We remark
remark that
that in
in this
this
algorithm,
algorithm, at
at each
each M
M iterations,
iterations, the
the weights
weights w
w used
used to
to compute
compute
kk ·· kkww are
are updated,
updated, so
so as
as to
to improve
improve the
the sparsity
sparsity inducing
inducing
property
property of
of the
the penalty
penalty function,
function, following
following [11].
[11]. We
We refer
refer to
to
[11]
for
more
intricate
update
strategies.
[11] for more intricate update strategies.

Algorithm
Algorithm 11 Multi-Channel
Multi-Channel Linear
Linear Prediction
Prediction
(r)
(r)
Require:
Require: Input
Input Data
Data yy ,, X,
X, aa constant
constant for
for weight
weight updates
updates
✏ >
> 00
1:
1: Initialize
Initialize w
wii ← 11 for
for all
all ii
2:
2: repeat
repeat
N L−1
L
N
X
X1
X
X
3:
for
.. Update
3:
for m
m == 11 :: M
M do
do
Update w
w every
every M
M iterations
iterations
(r)
(s) =
= yyrr(s)
(s) −
(s − dd − l)
l) ûû(r)
(l).
(30)
ŷŷrr(s)
yyii(s
(30)
ii (l).
4:
z̃z̃ii ← 22 soft(z,
4:
soft(z, α
↵ w)
w) − zz

i=1 l=0
l=0
i=1
5:
ũũ ← (I
5:
(I +
+X
XTT X)
X)−11 uu +
+X
XTT (y
(yrr −z̃)
z̃)
(r)
6:
z̃z̃ ← (1
6:
(1 − λ)) zz +
+ λ 22 (y
(y − X
X ũ)
ũ) − z̃z̃
where ûû(r)
denote filters
filters determined
determined by
by solving
solving
where
ii denote
7:
u
←
(1
−
2λ)
u
+
2λ
ũ
7:
u
(1
2
)
u
+
2
ũ
⇣
⌘
(r)
(r)
(r)
8:
z̃z̃ ← soft(z,
ûû1(r)
ûN
8:
soft(z, α
↵ w)
w)
1 ,, .. .. .. ,, û
N
!
!
9:
w
←
/(
+
.. update
9:
wii
✏/(✏ + |z
|zii|)
|) for
for all
all ii
update w
w
L
N
N
L
X
XX
X
10:
until
convergence
yyii(s
=
= arg
arg min
min P
P yyrr(s)
(s) −
(s − dd − l)
l) uuii(l)
(l) .. 10: until convergence
u
u11,...,u
,...,uN
N
11:
11: return
return ŷŷrr =
= soft(z,
soft(z, α
↵ w)
w)
i=1
i=1 l=0
l=0
(31)
(31)
We
We use
use aa weighted
weighted ``11 norm
norm as
as the
the sparsity
sparsity promoting
promoting function
function
P
P(·).
(·). More
More precisely,
precisely, for
for uu ∈
2C
Cnn,, we
we set
set
P
P(u)
(u) =
= kuk
kukww =
=

n
n
X
X
ii

w
wii |u
|uii|,|,

(32)
(32)

n
where
where w
w∈
2R
R+n+ is
is aa given
given weight
weight vector.
vector. We
We update
update the
the weight
weight
vector
vector w
w as
as the
the iterations
iterations progress,
progress, to
to improve
improve the
the sparsity
sparsity
promoting
promoting property
property of
of P
P..
We
We remark
remark that
that for
for MCLP,
MCLP, aa similar
similar formulation
formulation has
has been
been
discussed
discussed in
in [11],
[11], along
along with
with an
an algorithm
algorithm based
based on
on the
the
alternating
alternating direction
direction method
method of
of multipliers
multipliers (ADMM)
(ADMM) algorithm
algorithm
[3].
employ the
the Douglas-Rachford
Douglas-Rachford(DR)
(DR)algorithm
algorithm[14,
[14],
[3]. We employ
5]
[5]
solving
problem.ToTothat
thatend,
end, we
we first
first define the
for for
solving
thisthis
problem.
the
(i)
matrices
matrices X
X (i) for
for ii ∈
2 {1,
{1, .. .. .. ,, N
N}} such
such that
that
(i)

(i) = y (s − d − l), for l = 0, 1, . . . , L − 1,
X
i
Xs,l
d l), for l = 0, 1, . . . , L 1,
s,l = yi (s

(33)
(33)

and
and as
as ss ranges
ranges over
over the
the domain
domain of
of the
the observation
observation yyrr.. Notice
Notice
(i)
(i) is Toeplitz and has L columns and as many rows
that
X
that X is Toeplitz and has L columns and as many rows
(i)
as
as the
the number
number of
of samples
samples in
in yyrr.. By
By concatenating
concatenating X
X (i)’s,
’s, we
we
construct
the
data
matrix
X
as,
construct the data matrix X as,
⇥ (1)

(n) ⇤
X
(34)
X=
= X
(34)
X (1) ·· ·· ·· X
X (n) ..
Similarly,
Similarly, we
we construct
construct the
the variable
variable vector
vector uu as
as


2 u(1)
3
u(1)

... 
6
7,
uu =
=
4 ... 
5,
)
uu(N
(N )

(35)
(35)

B. DOA
DOA Estimation
Estimation
B.
Because audio
audio is
is aa wideband
wideband signal,
signal, narrowband
narrowband DOA
DOA
Because
estimation techniques
techniques like
like MUSIC
MUSIC [22],
[22], or
or ESPRIT
ESPRIT [21]
[21] are
are
estimation
not directly
directly applicable.
applicable. However,
However, simple
simple extensions
extensions of
of MUSIC
MUSIC
not
to wideband
wideband scenarios
scenarios exist
exist and
and we
we employ
employ one
one of
of them.
them. In
In aa
to
nutshell, the
the idea
idea is
is to
to estimate
estimate the
the DOA
DOA for
for each
each frequency
frequency
nutshell,
separately and
and average
average the
the estimates.
estimates. This
This approach
approach is
is referred
referred
separately
to as
as the
the ‘incoherent
‘incoherent signal
signal subspace
subspace method’
method’ (ISSM)
(ISSM) in
in the
the
to
recent literature
literature [32]
[32] –– see
forfor
early
work.
Despite
its
recent
see [30,
[30],25]
[25]
early
work.
Despite
simplicity,
ISSM
is noted
to perform
well well
in high
its
simplicity,
ISSM
is noted
to perform
in SNR
high [32].
SNR
Subspace methods for DOA estimation make use of two
[32].
essential
components,
namely
theestimation
covariancemake
matrixuse
of the
Subspace
methods for
DOA
of data
two
and the manifold
vector,
determined
by thematrix
array ofgeometry.
essential
components,
namely
the covariance
the data
For the
a fixed
frequency
!, determined
and time instance
s (in geometry.
the STFT
and
manifold
vector,
by the array
domain),
wefrequency
estimate the
matrix RsY(in
(s, the
!) as
For
a fixed
ω, correlation
and time instance
STFT
domain),
we!)estimate
correlation
ω) as (37)
RY (s,
= (1 the
⌘) M
(s, !) + ⌘matrix
RY (s RY1,(s,
!),

RY (s, ω) = (1 − η) M (s, ω) + η RY (s − 1, ω),
(37)
where
2
3
Ŷ1 (s, !)
where
⇤
6
7⇥ ⇤

M (s, !) = 4 Ŷ1 (s,... ω) 
5 Ŷ1 (s, !) · · · ŶN⇤ (s, !) . (38)


.. !)  ∗
M (s, ω) =  ŶN (s,
 Ŷ1 (s, ω) · · · ŶN∗ (s, ω) . (38)
.
In this definition,
⌘ 2
ŶN (s,
ω) [0, 1] is a smoothing (or forgetting)
factor justified by the assumption that the sources move slowly.
In Let
thisusdefinition,
∈ manifold
[0, 1] is avector.
smoothing
now defineη the
Recall (or
fromforgetting)
Sec. II-A
factor
justified
by
the
assumption
that
th the sources move slowly.
that pi denotes the position of the i microphone with respect

7

Let us now define the manifold vector. Recall from Sec. II-A
that pi denotes the position of the ith microphone with respect
to the center of the array. Also, let uθ denote the unit vector
in the direction θ. Finally, let c denote the speed of sound, fs
denote the sampling frequency and
τi,θ = −

hpi , uθ i
fs
c

(39)

denote the relative delay (in samples) of the ith microphone
for a source in the direction θ. The manifold vector is defined
as


exp −jτ1,θ ω


..
(40)
aθ (ω) = 
.
.

exp −jτN,θ ω

Assuming there are K sources, let CK (s, ω) denote the (N −
K) × N unitary matrix whose columns are the eigenvectors of
RY (s, ω) corresponding to the smallest N − K eigenvalues.
The DOA function is defined as
!−1
X
∗
D(s, θ) =
aθ (ω) CK (s, ω) aθ (ω)
.
(41)
ω

For each time instance s, the K dominant peaks of D(s, θ)
with respect to θ are taken to be the directions of the sources.
C. Source Separation
For source separation, we use the geometric source separation (GSS) method [18], [27] complemented with a simple
post-filter. We assume that the sources are stationary, so that
the estimated directions are constant with respect to time.
Extension to moving sources is straightforward, at least in
principle, and will not be further discussed.
1) Geometric Source Separation: Let θi , for i = 1, . . . , K
denote the directions of the sources determined by the method
described in Sec. IV-B. For each ω, we define a matrix A(ω)
as,


A(ω) = aθ1 (ω) · · · aθK (ω)
(42)
We model the MCLP outputs as




Ŷ1 (s, ω)
X1 (s, ω)




..
..
.

 = A(ω) 
.
.
|

ŶN (s, ω)
{z
}
Ŷ (s,ω)

|

(43)

XK (s, ω)
{z
}

for a given weight matrix W . Also, let RS (ω) denote the
empirical autocorrelation matrix of S(s, ω), obtained by averaging S(s, ω) S ∗ (s, ω) over the time variable s. Given this
notation, Ŵ (ω) is defined as the solution of
n
 2
 γ
RS (ω) − diag RS (ω) F
min C W (ω) =
2
W (ω)

1
2
+ W (ω) A(ω) − I F , (45)
2

where k · kF denotes the Frobenius norm.
Thanks to the differentiability of both terms in the cost
function, gradient descent can be applied for obtaining a
stationary point. For a fixed ω, we remark that the gradient
[4] of the complex valued C(W ) is given as [18],

∂C(W )
= γ RS −diag(RS ) W RŶ +(W A−I) A∗ , (46)
∗
∂W

where RŶ denotes the empirical autocorrelation matrix of Ŷ ,
obtained by averaging Ŷ ∗ (s, ω) Ŷ (s, ω) over the time variable
s.
2) Post-Filter: The GSS algorithm aims to achieve source
separation by considering the best linear combination of the
inputs. In practice, separation performance can be improved
by post-filtering. Specifically, we view the separated signals as
contaminated with noise and eliminate this noise using a softthreshold, where the time-frequency varying threshold value
depends on an estimate of the remaining noise variance for
each time-frequency sample.
More specifically, if Ŝ1 , . . . , ŜK denotes the K source estimates obtained by GSS, for each value of the triplet (i, s, ω),
we regard Ŝi (s, ω) as an unbiased estimate of Xi (s, ω)with
variance σi2 (s, ω). The post-filtered estimate is formed as,


X̂i (s, ω) = soft Ŝi (s, ω), α σi (s, ω) ,
(47)
where α is a parameter.
In order to realize such a post-filter, we need to estimate
σi2 (s, ω). We estimate σi2 (s, ω) as,

N
X
1
|Wi,n (s, ω)|
N − 1 n=1

2
× Ŷn (s, ω) − aθi ,n (ω) Ŝi (s, ω)
, (48)

where aθi ,n (ω) denotes the nth entry of aθi (ω). While this
expression may make intuitive sense, further justification is
provided in Appendix E. We also remark that this is a
generalized version of the post-filter employed in [1].

X(s,ω)

This model is not exact because despite the significant dereverberation by MCLP, Ŷ still contains reflections other than
the direct signal.
For this model, GSS proposes to estimate X(s, ω) as
Ŵ (ω) Ŷ (s, ω), where the K × N weight matrix Ŵ (ω) is
determined by solving a minimization problem. To describe
the problem, let us denote
S(s, ω) = W (ω) Y (s, ω),

(44)

V. E XPERIMENTS
Our primary claim in this paper is that an MCLP step suppresses reflections other than the direct source, under a certain
assumption (Assumption 1). The experiments in this section
are geared towards demonstrating the claims and verifying
Assumption 1. There are three experiments. In the first, we
verify Assumption 1 using real impulse responses obtained in
a classroom. In the second experiment, we demonstrate that
an MCLP preprocessing step eliminates reflections other than

8
88

m77
m

Normalized Variance Fields v̂(s, !) for Different DOAs

m66
m

m88
m

(a) DOA = 0 radians
−20
−20

✓✓

m55
m

m44
m

m11
m

m22
m

Frequency
(kHz)
Frequency (kHz)

88

−22
−22

66
5

−24
−24

44
3

−26
−26

22

−28
−28

1

m33
m

00

Fig.
Fig. 3.
3. We
We use
use an
an array
array of
of eight
eight microphones
microphones distributed
distributed uniformly
uniformly on
on aa circle.
circle.
✓θ denotes
denotes the
the angle
angle of
of arrival.
arrival. In
In particular,
particular, ✓θ =
= 00 radians
radians if
if the
the source
source is
is in
in
the
direction
of
m
,
✓
=
⇡/2
radians
if
the
source
is
in
the
direction
of
the direction of m111 , θ = π/2 radians if the source is in the direction of m
m777 ,,
etc.
etc.

0.05

0.1
0.1

0.15

0.2
0.25
0.3
0.2
0.3
Time (seconds)
Time (seconds)

0.35

0.4
0.4

0.45

(b) DOA = ⇡/4 radians
(b) DOA = ⇡/4 radians

8

Frequency
Frequency(kHz)
(kHz)

−30
−30

−20

8

−20

6

−22

7

−22

6

5

−24

−24

4

4

−26

3

−26

2

−28

2

−28

1
0

0.05

0

0.1

0.15

0.1

0.2
0.25
Time
0.2 (seconds)

0.3

0.35

0.3

0.4

0.45

0.4

−30

−30

Time (seconds)

(c) DOA = ⇡/2 radians
(c) DOA = ⇡/2 radians

8

−20

−20

87
Frequency
(kHz)
Frequency
(kHz)

claim
about
the effect
of employing
MCLP
preprocessing
the direct
source,
by showing
how theanDOA
function
changes
step
in
a
source
separation
problem.
after the MCLP step. Finally, the third experiment validates our
All of
the the
experiments
performed
real preprocessing
data. In order
claim
about
effect of are
employing
an on
MCLP
to
obtain
the
data,
we
used
a
circular
microphone
array with
step in a source separation problem.
eight
omni-directional
microphones.
The
configuration
and
the
All of the experiments are performed on real data. In
order
definition
of
DOA,
namely
✓,
are
as
shown
in
Fig.
3.
The
to obtain the data, we used a circular microphone array with
microphones
are distributed
uniformly
a circle ofand
radius
eight omni-directional
microphones.
Theon
configuration
the
5definition
cm.
of DOA, namely θ, are as shown in Fig. 3. The
In the experiments,
we use uniformly
an STFT with
microphones
are distributed
on aa smooth
circle ofwindow
radius
of
length
60
ms,
and
hop-size
is
15
ms.
The
window
is selected
5 cm.
such
thatexperiments,
the STFT is we
self-inverting
(i.e.,
a Parseval
In the
use an STFT
with
a smoothframe).
window
Some
of
the
audio
signals
used
in
the
experiments,
as well
of length 64 ms, and hop-size is 16 ms. The window is selected
as
those
from
experiments
not
contained
in
the
paper
can be
such that the STFT is self-inverting (i.e., a Parseval frame).
listened
to
at
“http://web.itu.edu.tr/ibayram/MCLP/”.
The sampling frequency for the audio signals is 16 KHz. For
the MCLP block,
the parameters
(see Algorithm
chosen
Experiment
1 (Verification
of Assumption
1).1)Inareorder
to
as d =Assumption
2, L = 30, M
= 50,
α=
0.05, λ =
0.5. We
run the
verify
1, we
used
estimated
impulse
responses
algorithm
for 2000
iterations.
from
different
DOAs.
The impulse responses are estimated
Some
the audio
signalsinused
the experiments,
well
using
theofmethod
described
[6], inusing
loudspeakers as
placed
as those
experiments
contained
in the
can we
be
at
knownfrom
DOAs
in the farnotfield.
Knowing
thepaper
DOAs,
listened to the
at “http://web.itu.edu.tr/ibayram/MCLP/”.
computed
distance compensated impulse responses and
obtained filters H̃ii (s, !) that approximate their action in the
Experiment 1 (Verification of Assumption 1). In order to
STFT domain. For each s and ! value, we computed the
verify Assumption 1, we used estimated impulse responses
variance of H̃ii (s, !)’s with respect to i (recall i takes eight
from different DOAs. The impulse responses are estimated
different values), which we denote as v(s, !). We normalize
using the method described in [6], using loudspeakers placed
v(s, !) by P
the total power at the location s, !, defined by
at known DOAs in the2 far field. Knowing the DOAs, we
p(s, !) = ii |Hii (s, !)|2 . Thus for each s and ! value, we
computed the distance compensated impulse responses and
obtained a number v̂(s, !) = v(s, !)/p(s, !). We remark that
obtained filters H̃i (s, ω) that approximate their action in the
Assumption 1 is valid if and only if v̂(0, !) = 0 for all !. In
STFT domain. For each s and ω value, we computed the
Fig. 4, we show (in decibels) the fields v̂(s, !) for different
variance of H̃i (s, ω)’s with respect to i (recall i takes eight
DOAs, namely 0, ⇡/4 and ⇡/2 radians. Observe that for all
different values), which we denote as v(s, ω). We normalize
DOAs and almost all ! shown, v̂(0, !) (the first column in
v(s, ω) by the total power at the location s, ω, defined by
each image)Pis significantly
lower than v̂(s, !) for s 6= 0 (the
p(s, ω) = i |Hi (s, ω)|2 . Thus for each s and ω value, we
remaining columns).
Therefore, Assumption 1 is seen to be
obtained a number v̂(s, ω) = v(s, ω)/p(s, ω). We remark that
valid for practical purposes.
Assumption 1 is valid if and only if v̂(0, ω) = 0 for all ω. In
Experiment
2 (Direction
of Arrival).
In v̂(s,
this ω)
experiment,
we
Fig. 4, we show
(in decibels)
the fields
for different
demonstrate
the
effect
of
an
MCLP
front-end
on
estimating
DOAs, namely 0, π/4 and π/2 radians. Observe that for all
the
DOA.
DOAs
and almost all ω shown, v̂(0, ω) (the first column in
Theimage)
sourceissignal
is produced
loudspeaker
each
significantly
lower with
than av̂(s,
ω) for s 6=located
0 (the
at
✓ = ⇡/2 columns).
radians in Therefore,
a reverberant
environment.
spectroremaining
Assumption
1 isThe
seen
to be
gram
of the
signal purposes.
observed at the first microphone is shown
valid for
practical
in Fig. 5a. Notice that the harmonics are seriously smeared.
Experiment 2 (Direction of Arrival). In this experiment, we
For this signal, MCLP significantly removes the effects of
demonstrate the effect of an MCLP front-end on estimating
reverberation – see Fig. 5b.

7

−22

6

−22

65

−24

−24

4

−26

43

−26

2

−28

21

−28
0

0

0.05

0.1

0.1

0.15

0.2
0.25
Time (seconds)

0.3

0.2
0.3
Time (seconds)

0.35

0.4

0.4

0.45

−30

−30

Fig. 4. Experiment 1 : ‘Correlation figures’ (see the text for a description),
in dB,
the circular
and different
directions
arrival.
DOA = 0
Fig.
4. for
Experiment
1 :array
‘Correlation
figures’
(see theoftext
for a(a)
description),
rad,dB,
(b)for
DOA
⇡/4 rad,
(c)and
DOA
= ⇡/2 directions
rad.
in
the =circular
array
different
of arrival. (a) DOA = 0
π/4 rad, (c) DOA = ⇡/2
π/2 rad.
rad, (b) DOA = ⇡/4

In order to visualize the effects of MCLP on DOA estimation,
we show
the DOA
✓) (defined
(41))
order
to visualize
thefunction
effects ofD(k,
MCLP
on DOA inestimatheInDOA.
averaged
over
the
time
variable
k
in
Fig.
6,
for
both
the
tion,
we
show
DOA
function
D(k,
✓)
(defined
in
(41))
The source signal is produced with a loudspeaker located
reverberant
observations
and
MCLP
outputs.
We
remark
that
averaged
over
the
time
variable
k
in
Fig.
6,
for
both
the
at θ = π/2 radians in a reverberant environment. The spectroaveraging
isobservations
meaningful
here,at
since
themicrophone
source
stationary.
reverberant
and
MCLP
outputs.
Weisremark
that
gram of the
signal
observed
the first
is shown
For
comparison,
we
also
normalize
the
averaged
DOA
funcaveraging
is
meaningful
here,
since
the
source
is
stationary.
in Fig. 5a. Notice that the harmonics are seriously smeared.
tions
so that
theirMCLP
unity. Thethe
correct
position
offuncthe
For
wepeak
alsois
normalize
averaged
For comparison,
this
signal,
significantly
removes
theDOA
effects
of
source
is
also
indicated
with
a
dotted
line.
We
see
that
tions
so
that
their
peak
is
unity.
The
correct
position
of
the
reverberation – see Fig. 5b.
DOA function
clearly improved
by the line.
application
MCLP.
source
is also isindicated
with a dotted
We seeofthat
the
In order to visualize the effects of MCLP on DOA estimaThe variance
reduced and
the application
function isofcloser
to
DOA
functionisis clearly improved
by the
MCLP.
tion, we show the DOA function D(k, θ) (defined in (41))
an impulse
at thereduced
correctand
position.
The
variancelocated
is clearly
the function is closer to
averaged over the time variable k in Fig. 6, for both the
an impulse located at the correct position.
reverberant observations
and MCLP outputs.
remark that
Experiment
3 (Source Separation).
In this We
experiment,
we
averaging is the
meaningful
here,gained
since the
source
is stationary.
demonstrate
improvement
by
MCLPwe
as
Experiment
3 (Source
Separation).
In employing
this
experiment,
comparison,
we also
the
aFor
pre-processing
step
in anormalize
blind
separation
setup. DOA
We
usefuncthe
demonstrate
the improvement
gained
by averaged
employing
MCLP
as
tions
so that theirstep
peak
unity.
The
correct
position
offrom
microphone
array
in
Exp.
2. We
use
the
source
asame
pre-processing
in as
ais blind
separation
setup.
We use
the
source2microphone
isasalso
a dotted
We
seeproduced
thatfrom
the
Exp.
the indicated
firstarray
source.
second
source
is
same
aswith
in The
Exp.
2. Weline.
use
the source
DOA
clearly
improved
the⇡ application
MCLP.
by
a loudspeaker
placed
at ✓The
= by
1.75
radians
(soof
that
the
Exp.
2function
as the isfirst
source.
second
source is
produced
Thea variance
is clearly
and
the
is closer
to
sources
are separated
byreduced
0.75
The spectrogram
of
by
loudspeaker
placed
at ✓ ⇡=radians).
1.75
⇡ function
radians
(so
that the
an impulse
at by
the
correct
position.
the
observation
at the
first
microphone
isThe
shown
in Fig. 7a.
sources
are located
separated
0.75
⇡ radians).
spectrogram
of
The observation
output of MCLP
the first microphone
as
the
at the when
first microphone
is shown is
in taken
Fig. 7a.
Experiment
3 (Source
In this experiment, we
the
reference
shown when
inSeparation).
Fig.the7b.first microphone
The
output ofisMCLP
is taken as
demonstrate the improvement gained by employing MCLP as
averageis (over
DOA
theThe
reference
showntime)
in Fig.
7b.functions obtained from the
a pre-processing step in a blind separation setup. We use the
original
recordings
and
MCLP
outputs
are shown
in Fig.
8. We
The average (over time) DOA functions
obtained
from
the
same microphone array as in Exp. 2. We use the source from
observe
that
MCLP
has
considerably
sharpened
the
functions.
original recordings and MCLP outputs are shown in Fig. 8. We

99
9

(a) Reverberant Observation
(a) Reverberant Observation

0
0
-10
-10

2
2

-20
-20
-20

1.5
1.5

-30
-30

1
1

-40
-40

0.5
0.5
0
00
0

1
1

2
2

3
3

5
5

Time (seconds)
Time (seconds)
(b) MCLP Output
(b) MCLP Output

3
3

Frequency
Frequency
(kHz)
Frequency(kHz)
(kHz)

4
4

6
6

7
7

2
2

-30
-30

1
1

-40
-40

0.5
0.5
1
1

2
2

3
3

4
4

5
5

6
6

Time (seconds)
Time (seconds)

7
7

Fig. 5. Experiment 2 : Spectrograms of the signals used in the experiment.
Fig. 5. Experiment 2 : Spectrograms of the signals used in the experiment.
(a) The observation at the first microphone. (b) MCLP output when the first
(a)
The observation
microphone. (b) MCLP output when the first
microphone
is used at
as the
the first
reference.
microphone is used as the reference.

-40
-40

0.5
0.5
0
00
0

1
1

2
2

3
3

0.75
0.75

6
6

7
7

-50
-50
-50

8
8

0
0
0
-10
-10

2
2

-20
-20

1.5
1.5

-30
-30
-30

1
1

-40
-40

0.5
0.5
0
00
0

1
1

2
2

3
3

4
4

5
5

6
6

Time (seconds)
Time (seconds)

7
7

-50
-50

8
8

Fig. 7. Experiment 3 : (a) Spectrogram of the two mixture signal observed
Fig. 7. Experiment
Experiment 33 :: (a)
(a) Spectrogram
Spectrogram of
of the
the two mixture signal observed
at
at the
the first
first microphone.
microphone. (b)
(b) MCLP
MCLP applied
applied to
to the
the observations,
observations, when
when the
the first
first
microphone
microphone is
is taken
taken as
as the
the reference.
reference.

DOA Functions
DOA Functions
Original
Original
MCLP output
MCLP output

0.75
0.75

0.5
0.5

5
5

2.5
2.5

1
1

Original
Original
MCLP output
MCLP output

4
4

Time (seconds)
Time (seconds)
(b) Dereverbed Observations
(b) Dereverbed Observations

DOA Functions
DOA Functions

1
1

0.5
0.5

0.25
0.25
0
00
0

-30
-30

1
1

-50
-50

8
8

-20
-20
-20

1.5
1.5

3
3

-20
-20

1.5
1.5

-10
-10

2
2

0
0
-10
-10

0
0

2.5
2.5

-50
-50
-50

8
8

2.5
2.5

0
00
0

Frequency
Frequency
(kHz)
Frequency(kHz)
(kHz)

2.5
2.5

(a) Reverberant Observation
(a) Reverberant Observation

3
3

Frequency
Frequency
(kHz)
Frequency(kHz)
(kHz)

Frequency
Frequency
(kHz)
Frequency(kHz)
(kHz)

3
3

0.25
0.25
0.2
0.2

0.4
0.4

0.6
0.6

0.8
0.8

1
1

θ/π
θ/π

1.2
1.2

1.4
1.4

1.6
1.6

1.8
1.8

2
2

0
00
0

0.2
0.2

0.4
0.4

0.6
0.6

0.8
0.8

1
1

θ/π
θ/π

1.2
1.2

1.4
1.4

1.6
1.6

1.8
1.8

2
2

Fig. 6. Experiment 2 : The DOA functions D(k, ✓) (see (41)) averaged over
Fig.
Experiment
: The DOA functions
D(k,
✓) (see (41))
averaged
time 6.
s) for the2
observations
and
outputs.
We over
also
Fig.
6.(i.e.,
Experiment
2 :reverberant
The DOA functions
D(k,
θ) MCLP
(see (41))
averaged
over
time (i.e., s) for the reverberant observations and MCLP outputs. We also
apply
a normalization
to each DOA
function and
so that
the outputs.
peak is equal
to
time (i.e.,
s) for the reverberant
observations
MCLP
We also
apply a normalization to each DOA function so that the peak is equal to
unity.
apply a normalization to each DOA function so that the peak is equal to
unity.
unity.

Fig. 8. Experiment 3 : The DOA functions D(k, ✓)
θ) (see (41)) averaged over
Fig. 8. Experiment 3 : The DOA functions D(k, ✓) (see (41)) averaged over
time (i.e., k) for the reverberant observations and MCLP outputs. We also
time
k) for the reverberant
observations
MCLP
apply(i.e.,
a normalization
to each DOA
function and
so that
the outputs.
peak is We
equalalso
to
apply a normalization to each DOA function so that the peak is equal to
unity.
unity.

observe that MCLP has considerably sharpened the functions.
observe
thatthe
MCLP
considerably
sharpened
Exp.
as
first has
source.
The second
sourcethe
is functions.
produced
Also, 2without
MCLP,
the peaks
of the DOA
functions
are
Also,
without
MCLP,
the
peaks
of
the
DOA
functions
are
by
a loudspeaker
placed
at the
θ =
1.75directions,
π radians which
(so that
the
offset
by some amount
from
actual
would
offset byare
some
amountby
from
theπactual
directions,
which would
sources
separated
0.75
radians).
The spectrogram
of
lead
to
aa biased
estimate
of
the
true
DOA.
lead
to
biased
estimate
of
the
true
DOA.
theThe
observation
at the first
microphone
Fig. 7a.
source separation
results
are shownis inshown
Fig. 9.inPerceptuThe source
separation
results
are
shown in Fig. 9.isPerceptuThe
MCLPespecially
when
thefor
first
as
ally, output
withoutofMCLP,
themicrophone
second source,taken
we can
ally,
without
MCLP,
especially
for
the
second source, we can
the
reference
is shown
in Fig. despite
7b.
clearly
hear the
first source,
the post-filter. Further,
clearly
hear the(over
first time)
source,
despite
the post-filter.
Further,
The average
DOA
functions
obtained
from
the
reverberation
is not eliminated.
These
artifacts
are eliminated
reverberation
is
not
eliminated.
These
artifacts
are
eliminated
original
recordings and MCLP outputs are shown in Fig. 8. We
using MCLP.
using MCLP.
observe
thattoMCLP
has the
considerablywesharpened
thePESQ
functions.
In
also
report
In order
order
to quantify
quantify
the results,
results,
wethe
also
report
PESQ and
and
Also,
without
MCLP,
the
peaks
of
DOA
functions
are
SIR
scores
[29].
PESQ
gives
an
objective
perceptual
distance
SIR
scores
[29].
PESQ
gives
an
objective
perceptual
distance
offset
by
some
amount
from
the
actual
directions,
which
would
between
an input
speech
and
aa reference.
between
input
speechofsignal
signal
andDOA.
reference. In
In order
order
lead
to a an
biased
estimate
the true
to
obtain
meaningful
distances,
we
use
the
multiple
channel
to The
obtain
meaningful
distances,
we
use
the
multiple
channel
source
separation
results
are shown
inmixing,
Fig. 9. Percepturecordings
of source-1
and
source-2,
prior to
to extract
recordings of source-1 and source-2, prior to mixing, to extract

ally, without MCLP, especially for the second source, we can
th
the
signals.
take
as
clearly
hear the
first We
source,
despite
the post-filter.
Further,
the reference
reference
signals.
We
take the
the iith microphone
microphone
as reference
reference
and
apply
MCLP
to
the
multiple
microphone
recordings
reverberation
is
not
eliminated.
These
artifacts
are
eliminated
and apply MCLP to the multiple microphone recordings to
to
obtain
the
using
obtainMCLP.
the references
references for
for each
each source.
source. We
We then
then compute
compute the
the
PESQ
scores
for
described
this
paper.
We
In order
to quantify
the results,
we alsoin
PESQ
scores
for the
the methods
methods
described
inreport
this PESQ
paper. and
We
finally
average
over
i
to
obtain
a
reliable
PESQ
score.
The
SIR
scores
[29].over
PESQ
an objective
distance
finally
average
i togives
obtain
a reliableperceptual
PESQ score.
The
resulting
scores
are
tabulated
in
remark
the
between
input
signal
and I.
a We
reference.
resulting an
scores
arespeech
tabulated
in Table
Table
I.
We
remarkInthat
thatorder
the
inclusion
of
an
MCLP
step
(denoted
as
LP
in
the
table)
to
obtain meaningful
distances,
we use as
the LP
multiple
inclusion
of an MCLP
step (denoted
in thechannel
table)
improves
PESQ
noticeably.
recordings
of source-1
and source-2,
prior to mixing, to extract
improves the
the
PESQ scores
scores
noticeably.
th
theThe
reference
signals.
We
take
the
i
microphone as reference
SIR
for
separation
with
different
of
The
SIR scores
scores
for
separation
withmicrophone
different combinations
combinations
of
and
apply
MCLP
to
the
multiple
recordings
to
the
discussed
methods
are
also
provided
in
Table
I.
We
observe
the
discussed
methods
are
also
provided
in
Table
I.
We
observe
obtain
the
references
for
each
source.
We
then
compute
the
that
including
an
MCLP
preprocessing
step
also
improves
that
including
an the
MCLP
preprocessing
step
alsopaper.
improves
PESQ
scoresthe
for
methods
described
in by
this
We
significantly
separation
quality
measured
SIR.
significantly the separation quality measured by SIR.

10
10

TABLE
TABLE II
PESQ
PESQ AND
AND SIR
SIR SCORES
SCORES FOR
FOR E
EXP
XP.. 3.
3. (LP
(LP :: M
MULTICHANNEL
ULTICHANNEL L
LINEAR
INEAR
PPREDICTION
REDICTION,, G
G :: G
GEOMETRIC
EOMETRIC S
SOURCE
OURCE S
SEPARATION
EPARATION,, P
P :: PPOST
OST FILTER
FILTER))

(a) Source-1 : GSS + PF

3

0

Frequency (kHz)

2.5

-10

2
-20

-30

1
-40

0

-50

0

1

2

3

4

5

6

7

8

Time (seconds)
(b) Source-2 : GSS + PF

3

0

Frequency (kHz)

2.5

-10

2
-20

1.5
-30

1
-40

0.5
0

-50

0

1

2

3

4

5

6

7

8

Time (seconds)
(c) Source-1 : MCLP + GSS + PF

3

0

2.5

Frequency (kHz)

Obs.
Obs.
1.07
1.07
1.04
1.04

LP
LP
1.04
1.04
1.08
1.08

1.5

0.5

-10

2
-20

1.5
-30

1
-40

0.5
0

-50

0

1

2

3

4

5

6

7

8

Time (seconds)
(d) Source-2 : MCLP + GSS + PF

3

0

2.5

Frequency (kHz)

Reference
Reference
Source-1
Source-1
Source-2
Source-2

-10

2
-20

1.5
-30

1
-40

0.5
0

-50

0

1

2

3

4

5

6

7

8

Time (seconds)

Fig. 9. Experiment 3 : Spectrograms of the sources separated with and without
Fig. 9. Experiment 3 : Spectrograms of the sources separated with and without
MCLP. (a,b) Source-1,2 obtained with GSS and a post-filter, without MCLP.
MCLP. (a,b) Source-1,2 obtained with GSS and a post-filter, without MCLP.
(c,d) Source-1,2 obtained when MCLP is applied first.
(c,d) Source-1,2 obtained when MCLP is applied first.

G
G
5.0
5.0

G+P
G+P
13.0
13.0

PESQ
PESQ
G
G+P
G
G+P
1.05
1.05 1.08
1.08
1.05
1.05 1.14
1.14
SIR
SIR
LP+G
LP+G
7.8
7.8

LP+G
LP+G
1.10
1.10
1.22
1.22

LP+G+P
LP+G+P
1.20
1.20
1.49
1.49

LP+G+P
LP+G+P
19.5
19.5

finally average over iVI.
to obtain
a reliable PESQ score. The
C ONCLUSION
resulting
scores
are
tabulated
in
Table
WeMCLP
remarkasthat
the
In this paper, we advocated the useI. of
a preinclusion
of
an
MCLP
step
(denoted
as
LP
in
the
table)
processing step for blind reverberant source separation. We
improves
PESQ scores
noticeably.
provided the
theoretical
justification
for such a scheme and
The
SIR
scores
for
separation
combinations
of
demonstrated numerically that itwith
can different
significantly
improve the
the
discussed
methods
are
also
provided
in
Table
I.
We
observe
performance of readily available separation algorithms.
thatThe
including
an MCLPforpreprocessing
step also improves
current algorithm
MCLP is computationally
intensignificantly
the
separation
quality
measured
SIR. moving
sive. However, for a practical scenario withby
slowly

sources, if we have approximate MCLP weights at a certain
time, we think it would
beONCLUSION
possible to update the weights
VI. C
in a computationally efficient manner by warm-starting the
In this paper, we advocated the use of MCLP as a preoptimization procedure. Therefore, the proposed scheme has
processing step for blind reverberant source separation. We
potential for real-time applications as well.
provided theoretical justification for such a scheme and
Another issue involves the correlation between the weights
demonstrated numerically that it can significantly improve the
obtained for different reference microphones. In the current
performance of readily available separation algorithms.
formulation, MCLP is applied independently for each choice
The current algorithm for MCLP is computationally intenof the reference microphone. However, we would expect
sive. However, for a practical scenario with slowly moving
some correlation between these weights, especially if the
sources, if we have approximate MCLP weights at a certain
microphones are close with respect to the room dimensions,
time, we think it would be possible to update the weights
which is usually the case for a microphone array. We plan to
in a computationally efficient manner by warm-starting the
explore the implications of such correlation in future work.
optimization procedure. Therefore, the proposed scheme has
potential for real-time applications as well.
PPENDIX A
Another issue involvesAthe
correlation between the weights
P
ROOF OF P ROP. 1
obtained for different reference microphones. In the current
Without loss
of generality,
r = 1 infortheeach
following.
formulation,
MCLP
is appliedsuppose
independently
choice
H2 , . . . , microphone.
HN are co-prime,
find expect
causal
of Since
the reference
However,wewecanwould
polynomials
Ũ2 , . .between
. , ŨN such
thatweights, especially if the
some
correlation
these
microphones
are close with respect to the room dimensions,
N
X
which is Ũ
usually
case for a microphone array. We plan(49)
to
(z)
Hthe
k
k (z) = 1.
explore
k=2the implications of such correlation in future work.
Now set, U1 (z) = 1 and
A PPENDIX A
H1 (z) H1 (1)
P ROP.. 1
Uk (z) = Ũk (z) P ROOF OF
(50)
z 1
Without loss of generality, suppose r = 1 in the following.
Then for G as in (13), and d = 1, we have
Since H2 , . . . , HN are co-prime, we can find causal
N
X
polynomials
Ũ2 , . . . , ŨN such that
Gk (z) Hk (z) = H1 (z)
H1 (z) H1 (1)
(51)
N
X
k=1
Ũk (z) Hk (z) = 1.
(49)
N
X
k=2
⇥
Ũk (z) Hk (z)
(52)
k=2
Now set, U1 (z) = 1 and
= H1 (z)
H1 (z)
H1 (z) − H1 (∞)
Uk (z) = Ũk (z) = H1 (1).
.
z −1

H1 (1)

(53)

(54)
(50)

11

Then for G as in (13), and d = 1, we have
N
X

k=1


Gk (z) Hk (z) = H1 (z) − H1 (z) − H1 (∞)
×

N
X

k=2

Ũk (z) Hk (z)

= H1 (z) − H1 (z) − H1 (∞)
= H1 (∞).

(51)
(52)



(53)
(54)

For the converse, observe that for causal Uk and Hk , the
polynomial
C(z) =

N
X

k=1

z −d Uk (z) Hk (z)

(55)

does not have any constant terms, i.e., C(∞) = 0. Observe
also that G1 (∞) = 1. Thus if (15) holds, then
N
X

k=1

Gk (z) Hk (z)

z=∞

= G1 (∞) H1 (∞) − C(∞)

(56)

= H1 (∞).

(57)

A PPENDIX B
C ANONICAL F ORM OF A P OLYNOMIAL M ATRIX
In this appendix, we briefly present, for the sake of completeness, the (Smith) canonical form of a polynomial matrix, along with some of its properties. For a more detailed
exposition, we refer to the fine treatment in Sec. VI.1–3 of
Gantmacher’s book [8].
We refer to a quantity or matrix as non-zero constant if it
is a polynomial of order zero.
Let us start by introducing three elementary row operations
on a polynomial matrix A(z) :
1) Multiplication of a row with a non-zero constant.
2) Addition, to the j th row, the product of the ith row and a
polynomial p(z).
3) Interchange of two rows.
These operations can be performed by left multiplication with
a square matrix whose determinant is a non-zero constant. We
refer to matrices that can be used to perform elementary operations as elementary matrices in this appendix. We also remark
that similar operations can be defined on the columns instead
of rows. Elementary column operations can be performed by
right-multiplication with elementary matrices.
We say that two polynomial matrices A and A0 are equivalent if A can be transformed to A0 by elementary row/column
operations.
Let us now introduce the (Smith) canonical form.
Definition 1. An N × K (with K ≤ N ) diagonal polynomial
matrix A(z) of rank r with diagonal ai (z) for i ∈ {1, . . . , K}
is said to be in (Smith) canonical form if
1) ai (z) 6= 0 if i ≤ r, ai (z) = 0 if i > r,
2) ai (0) = 1 for i ≤ r,
3) ai (z) divides ai+1 (z) for i < r.

A fundamental result is that we can transform any polynomial matrix to one in canonical form via elementary
row/column operations.
Proposition 5. Each N × K polynomial matrix is equivalent
to a unique matrix in canonical form.
In general, the canonical form is reached by performing
both row and column operations. However, we are interested
in row operations only, for the development in this paper. We
have the following result.
Proposition 6. Suppose an N × K matrix A has rank K and
the diagonal of its canonical form is constant. Then, A can be
reduced to its canonical form via row operations.
This result is a corollary of the discussion in [8] but is
not stated explicitly. We provide a proof for the sake of
completeness. We need a few intermediate results for the
proof. All of the non-obvious claims stated before the proof
of Prop. 6 can be found in [8].
Given an N × K polynomial matrix A(z), let Dr (z) denote
the greatest common divisor of all minors of A of order r,
normalized such that the leading coefficient is unity. Notice
that Dr is divisible by Dr−1 . Therefore, for r > 1, the ratios
ir = Dr /Dr−1 are polynomials. We also set i1 = D1 . It can
be shown that Dr ’s are invariant under elementary row/column
operations. This property transfers to ir as a consequence of
its definition.
Suppose now that the canonical form of A(z) consists of
a diagonal matrix with elements ar (z) as in Defn. 1. It turns
out that ar = ir , giving us an alternative definition of ar .
We need a final result before proving Prop. 6, which can be
shown by a simple elimination argument.
Proposition 7. Using elementary row operations, any polynomial matrix B(z) can be reduced to an upper triangular matrix
U (z) such that, at each column of U (z), if the polynomial
on the diagonal is non-zero, then the degree of any offdiagonal polynomial is strictly smaller than the degree of the
polynomial on the diagonal.
We are now ready for the proof of Prop. 6.
Proof of Prop. 6. Suppose the canonical form of A is B.
Observe that the only non-zero minor of B of order K is given
by the product of its diagonals, which is unity. Therefore, the
invariant polynomial DK of A is a constant.
Suppose we reduce A to an upper triangular form U as in
Prop. 7 using elementary row operations only. Observe that
the only non-zero minor of U of order K is given by the
product of its diagonals. Let us denote this non-zero minor by
UK (z). But since this is the only non-zero minor of order
K, and elementary matrices have constant determinants, it
follows that we must have UK (z) = c DK for some constant c.
Thus, UK (z) is a constant. Therefore, U is a constant diagonal
matrix with a non-zero diagonal. Each constant on the diagonal
can be made unity by an elementary row operation, yielding
the canonical form B.
In a final proposition, we show the equivalence stated in
Remark 1.

12

Proposition 8. Suppose A is an N × K polynomial matrix
with K ≤ N . The greatest common divisor of all the Kminors of A is a non-zero constant if and only if the Smith
canonical form of A has a non-zero constant diagonal.
Proof. As in the discussion above, let Dr (z) denote the
greatest common divisor of all minors of A of order r, where
the leading coefficient is unity, ir = Dr /Dr−1 for r > 1,
with i1 = D1 . Recall that ir ’s form the diagonal in the Smith
canonical form of A(z).
Suppose now that DK (z) = 1. This implies that ir = 1 for
all r ≤ K. Therefore the diagonal of the canonical form of A
is the identity matrix.
For the converse, suppose the canonical form of A has a
constant diagonal. This means that ir (z)’s are constant for
r ≤ K. But this implies that Dr (z) = 1 for r ≤ K.
A PPENDIX C
P ROOF OF P ROP. 3

A PPENDIX D
P ROOF OF P ROP. 4
We remark that (28) follows by expressing (25) in the timedomain.
Let us consider (29). If (25) holds, and Gi ’s satisfy (13),
then by (26), we have
N
X


gi ∗ hi,k (s) = hr,k (0) δ(s).

N
X


gi ∗ hi,k (s) = c0 exp(−j ω τr,k ) δ(s),

i=1

By Assumption 1 and (27), it follows that

i=1

i=1

=
=

K
X

k=1
K
X

k=1

xk (s) ∗

N
X
i=1

gi (s) ∗ hi,k (s)



xk (s) ∗ c0 exp(−j ω τr,k ) δ(s)

k=1
K
X
0

=c

(58)

V H̄1 = I.

xk (s) exp(−j ω τr,k ).

k=1

For k = 1, . . . , K, let us define the polynomials
Pk =
and set

U2

H1,k (z) − H1,k (∞)
z −1
 
UN = P1

···

···

(59)


PK V.

(60)

Also, set U1 (z) = 0. Then for G defined as in (13), we have

G1

···




GN H = H1,1 (z) · · · H1,K (z)


−z −1 P1 (z) · · · PK (z) (61)


= H1,1 (∞) · · · H1,K (∞)

For the converse, observe that, for causal Uk and Hi,k , the
polynomial
Ck (z) =

N
X
i=1

z −d Ui (z) Hi,k (z)

(62)

does not have any constant terms, regardless of k, i.e.,
Ck (∞) = 0 for all k. Observe also that G1 (∞) = 1. Thus
if (25) holds, then,
N
X
i=1

Gi (z) Hi,k (z)

z=∞

= G1 (∞) H1,k (∞) − Ck (∞)
(63)
= H1,k (∞).

(64)

(66)

where c0 is independent of r and k. Using (22), and the
properties of convolution, we finally obtain
!
N
N
K
X
X
X
gi (s) ∗ yi (s) =
gi (s) ∗
xk (s) ∗ hi,k (s)
i=1

Without loss of generality, suppose r = 1 in the following.
By Remark 1, which is proved in Prop. 8, we first note that
the gcd of the K-minors of H̄1 is a non-zero constant if and
only if the canonical form of H̄1 has a non-zero constant
diagonal. Thus, by Prop. 6, we can find a K × (N − 1)
polynomial matrix V such that

(65)

(67)

A PPENDIX E
M OTIVATION FOR THE P OST-F ILTER IN S EC . IV-C2
In order to justify the variance expression in (48), we will
make use of Prop. 9 below. For simplicity, this proposition
uses a seemingly independent notation. We will identify the
elements in the proposition with our separation scenario right
after the proposition. Recall that a complex valued Gaussian
random variable Z is said to be circular [19] and have variance
2 σ 2 if its real and imaginary parts are independent and
Gaussian random variables with variance σ 2 .
Proposition 9. Suppose γ is an unknown constant and we
observe Yn = cn γ +un Zn , for n = 1, . . . , N , where |cn | = 1,
un ∈ C and Zn ’s are circularly symmetric, independent, unit
variance, zero-mean Gaussian random variables.
The uniformly minimum variance unbiased (UMVU) estimator for γ is of the form
γ̂ =

N
X

(68)

wn Yn

n=1

where

wn = c∗n β |un |−2 ,

(69)

with
β=

N
X

n=1

−2

|un |

!−1

.

(70)

13

We also have that,
N
X

1
|wn | Yn − cn γ̂
N − 1 n=1

2

(71)

is an unbiased estimator for the variance of γ̂.
Proof (Sketch). It can be shown, by taking into account the
Gaussian pdf of Yn , that a complete sufficient statistic for γ
is
U=

N
X

n=1

c∗n |un |−2 Yn .

(72)

Observing that E(U ) = β −1 γ, it follows by the RaoBlackwell theorem that the UMVUE is as given in (68) for
wn as in (69) and β as in (70). It can also be checked that the
variance of γ̂ is,
var(γ̂) =

N
X

n=1

|wn |2 |un |2 = β 2

Let us now set ν =
E(ν) =

N
X

n=1

=

N
X

n=1

=

N
X

n=1

=
=

N
X

n=1
N
X
n=1

PN

n=1

n=1

|un |−2 = β.

Yn − cn γ̂



2
|wn | E Yn − cn γ̂


|wn |  1 − |wn |


|wn |  1 − |wn |
|wn |

|wn |

N
X



2
2

|un |2 +

X

k6=n

|un |2 + β



2

(73)

. Then,



|wk |2 |uk |2 

X

k6=n



|wk |


2 β
+ β(1 − |wn |)
1 − |wn |
|wn |

β(1 − |wn |)

= β(N − 1).

(74)

In order to relate this proposition with the discussion
in Section IV-C2, we make the following identifications :
Yn → Ŷn (s, ω); cn → aθ,n ; wn → Wi,n ; γ̂ → Ŝi (s, ω).
For the ith source, if we regard all of the remaining sources
as noise terms with a Gaussian distribution, then the optimal
estimate of the ith source is obtained by a linear combination
of the observations. We assume that GSS approximately finds
these optimal weights. In that case, the expression in (48) gives
an estimate of the variance of these unwanted components (or,
‘noise’) that remain after GSS.
R EFERENCES
[1] İ. Bayram. A multichannel audio denoising formulation based on
spectral sparsity. IEEE/ACM Transactions on Audio, Speech, and
Language Processing, 23(12):2272–2285, Dec 2015.
[2] J. Bitzer and K. U. Simmer. Superdirective microphone arrays. In
M. Brandstein and D. Ward, editors, Microphone Arrays. Springer, 2001.

[3] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed
optimization and statistical learning via the alternating direction method
of multipliers. Foundations and Trends in Machine Learning, 3(1):1–
122, 2011.
[4] D. H. Brandwood. A complex gradient operator and its application
in adaptive array theory. Microwaves, Optics and Antennas, IEE
Proceedings H, 130(1):11–16, February 1983.
[5] P. L. Combettes and J.-C. Pesquet. Proximal splitting methods in signal
processing. In H. H. Bauschke, R. S. Burachik, P. L. Combettes,
V. Elser, D. R. Luke, and H. Wolkowicz, editors, Fixed-Point Algorithms
for Inverse Problems in Science and Engineering. Springer, New York,
2011.
[6] A. Farina. Simultaneous measurement of impulse response and distortion
with a swept-sine technique. In Proc. 108th Convention of the AES,
2000.
[7] S. Gannot and I. Cohen. Adaptive beamforming and postfiltering. In
Handbook of Speech Processing. Springer, 2008.
[8] F. R. Gantmacher. The Theory of Matrices - I. Chelsea Publishing
Company, 1959.
[9] A. Hyvärinen, J. Karhunen, and E. Oja. Independent Component
Analysis. Wiley, 2001.
[10] A. Jukić, T. van Waterschoot, T. Gerkmann, and S. Doclo. Multichannel linear prediction-based speech dereverberation with sparse
priors. IEEE/ACM Transactions on Audio, Speech, and Language
Processing, 23(9):1509–1520, Sept 2015.
[11] A. Jukić, T. van Waterschoot, T. Gerkmann, and S. Doclo. A general
framework for incorporating time - frequency domain sparsity in multichannel speech dereverberation. J. Audio Eng. Soc, 65(1/2):17–30,
2017.
[12] K. Kinoshita, M. Delcroix, T. Nakatani, and M. Miyoshi. Suppression of
late reverberation effect on speech signal using long-term multiple-step
linear prediction. IEEE Transactions on Audio, Speech, and Language
Processing, 17(4):534–545, May 2009.
[13] M. Kowalski, E. Vincent, and R. Gribonval. Beyond the narrowband approximation: Wideband convex methods for under-determined
reverberant audio source separation. IEEE Trans. Audio, Speech and
Language Proc., 18(7):1818–1829, Sept 2010.
[14] P. L. Lions and B. Mercier. Splitting algorithms for the sum of two
nonlinear operators. SIAM Journal on Numerical Analysis, 16(6):964–
979, 1979.
[15] M. Miyoshi and Y. Kaneda. Inverse filtering of room acoustics. IEEE
Transactions on Acoustics, Speech, and Signal Processing, 36(2):145–
152, Feb 1988.
[16] T. Nakatani, T. Yoshioka, K. Kinoshita, M. Miyoshi, and B. H. Juang.
Speech dereverberation based on variance-normalized delayed linear prediction. IEEE Transactions on Audio, Speech, and Language Processing,
18(7):1717–1731, Sept 2010.
[17] P. A. Naylor, E. A. P. Habets, N. D. Gaubitch, J. Y.-C. Wen, M. R. P.
Thomas, S. Gannot, A. W. H. Khong, J. Hopgood, C. Evers, S. Fortune,
M. Miyoshi, J. Buchner, and W. Kellermann. Speech Dereverberation.
Springer, 2010.
[18] L. C. Parra and C. V. Alvino. Geometric source separation: Merging
convolutive source separation with geometric beamforming. IEEE
Transactions on Speech and Audio Processing, 10(6):352–362, Sep
2002.
[19] B. Picinbono. On circularity. IEEE Transactions on Signal Processing,
42(12):3473–3482, Dec 1994.
[20] J. P. Reilly, M. Wilbur, M. Seibert, and N. Ahmadvand. The complex subband decomposition and its application to the decimation of
large adaptive filtering problems. IEEE Trans. Signal Processing,
50(11):2730–2743, Nov 2002.
[21] R. Roy and T. Kailath. ESPRIT-Estimation of signal parameters via
rotational invariance techniques. IEEE Transactions on Acoustics,
Speech, and Signal Processing, 37(7):984–995, Jul 1989.
[22] R. Schmidt. Multiple emitter location and signal parameter estimation.
IEEE Transactions on Antennas and Propagation, 34(3):276–280, Mar
1986.
[23] K. U. Simmer, J. Bitzer, and C. Marro. Post-filtering techniques. In
M. Brandstein and D. Ward, editors, Microphone Arrays. Springer, 2001.
[24] M. Souden, S. Araki, K. Kinoshita, T. Nakatani, and H. Sawada. A
multichanneMMSE-based framework for speech source separation and
noise reduction. IEEE Transactions on Audio, Speech, and Language
Processing, 21(9):1913–1928, Sept 2013.
[25] G. Su and M. Morf. The signal subspace approach for multiple wideband emitter location. IEEE Transactions on Acoustics, Speech, and
Signal Processing, 31(6):1502–1522, Dec 1983.

14

[26] M. Togami, Y. Kawaguchi, R. Takeda, Y. Obuchi, and N. Nukaga.
Optimized speech dereverberation from probabilistic perspective for time
varying acoustic transfer function. IEEE Transactions on Audio, Speech,
and Language Processing, 21(7):1369–1380, July 2013.
[27] J. M. Valin, J. Rouat, and F. Michaud. Enhanced robot audition based on
microphone array source separation with post-filter. In 2004 IEEE/RSJ
International Conference on Intelligent Robots and Systems (IROS)
(IEEE Cat. No.04CH37566), volume 3, pages 2123–2128 vol.3, Sept
2004.
[28] J.-M. Valin, J. Rouat, and F. Michaud. Microphone array post-filter for
separation of simultaneous non-stationary sources. In Proc. IEEE Int.
Conf. on Acoustics, Speech and Signal Proc. (ICASSP), 2004.
[29] E. Vincent, R. Gribonval, and C. Fevotte. Performance measurement
in blind audio source separation. IEEE Transactions on Audio, Speech,
and Language Processing, 14(4):1462–1469, July 2006.
[30] M. Wax, Tie-Jun Shan, and T. Kailath. Spatio-temporal spectral analysis
by eigenstructure methods. IEEE Transactions on Acoustics, Speech, and
Signal Processing, 32(4):817–827, Aug 1984.
[31] O. Yilmaz and S. Rickard. Blind separation of speech mixtures via
time-frequency masking. IEEE Transactions on Signal Processing,
52(7):1830–1847, July 2004.
[32] Y.-S. Yoon, L. M. Kaplan, and J. H. McClellan. TOPS: New DOA
estimator for wideband signals. IEEE Transactions on Signal Processing,
54(6):1977–1989, June 2006.

