1

Performance Analysis of Source Image
Estimators in Blind Source Separation
Zbyněk Koldovský and Francesco Nesta
1 Faculty

of Mechatronics, Informatics, and Interdisciplinary Studies, Technical University of
Liberec,

Studentská 2, 461 17 Liberec, Czech Republic. E-mail: zbynek.koldovsky@tul.cz,

arXiv:1603.04179v3 [cs.SD] 22 May 2017

fax:+420-485-353112, tel:+420-485-353534
2 Conexant

System, 1901 Main Street, Irvine, CA (USA). E-mail: francesco.nesta@conexant.com

Abstract
Blind methods often separate or identify signals or signal subspaces up to an unknown scaling factor.
Sometimes it is necessary to cope with the scaling ambiguity, which can be done through reconstructing signals as
they are received by sensors, because scales of the sensor responses (images) have known physical interpretations.
In this paper, we analyze two approaches that are widely used for computing the sensor responses, especially,
in Frequency-Domain Independent Component Analysis. One approach is the least-squares projection, while
the other one assumes a regular mixing matrix and computes its inverse. Both estimators are invariant to the
unknown scaling. Although frequently used, their differences were not studied yet. A goal of this work is to
fill this gap. The estimators are compared through a theoretical study, perturbation analysis and simulations. We
point to the fact that the estimators are equivalent when the separated signal subspaces are orthogonal, and vice
versa. Two applications are shown, one of which demonstrates a case where the estimators yield substantially
different results.
Index Terms
Beamforming, Blind Source Separation, Independent Component Analysis, Principal Component Analysis, Independent
Vector Analysis

I. I NTRODUCTION
The linear instantaneous complex-valued mixture model
x = Hs

or

X = HS

(1)

describes many situations where multichannel signals are observed, especially those considered in the field of
array processing [1] and Blind Source Separation (BSS) [2], [3]. The former equation is a vector-symbolic
description of the model while the latter equation describes a batch of data.
The vector x = [x1 , . . . , xd ]T represents d observed signals on sensors, s = [s1 , . . . , sr ]T represents original
signals, and H is a d × r complex-valued mixing matrix representing the linear mixing system. Upper-case

2

bold letters such as X and S will denote matrices whose columns contain concrete samples of the respective
signals; let the number of columns (samples) be N where N  d. We will focus on the regular case when the
number of the observed signals d is the same as that of the original signals r, but later in the article we will
also address an undetermined case where r > d. From this point forward, let H be a d × d full rank matrix.
Consider a situation where only a subset of the original signals is of primary interest (e.g., only one
particular source or the subspace spanned by some sources, so-called multidimensional source). Without a loss
of generality, let s be divided into two components [s1 ; s2 ] where the sub-vectors s1 and s2 have, respectively,
length m and d − m, 1 ≤ m < d. The former component will be referred to as target component, and the
latter as interference. Correspondingly, let H be divided as [H1 H2 ] where the sub-matrices H1 and H2 have
dimensions d × m and d × (d − m), respectively. Then (1) can be written as
x = H1 s1 + H2 s2 .

(2)

The terms H1 s1 and H2 s2 correspond to the contributions of s1 and s2 , respectively, for the mixture x, and
4

will be denoted as si = Hi si , i ∈ {1, 2}. If, for example, s2 is not active, then x = s1 , which is equal to the
observations of s1 on the sensors, that is, the sensor response or source image of s1 .
In audio applications, s1 is often a scalar signal (m = 1) representing a point source located in the room,
and the model (1) describes linear mixing in the frequency domain for a particular frequency bin [4], [5], [6],
[7]. For m > 1, s1 can correspond to a subgroup of speakers [8]. In biomedical applications, s1 or s2 can
consist of components related to a target activity such as muscular artifacts in electroencephalogram (EEG) [9],
maternal or fetal electrocardiogram (ECG) [10], and so forth.
The problem of retrieving si from x is often solved with the aid of methods for Blind Source Separation. The
objective of BSS is to separate the original signals based purely on their general properties (e.g., independence,
sparsity or nonnegativity). In a general sense, BSS involves Principal Component Analysis (PCA), Independent
Component Analysis (ICA) [2], [11], Independent Vector Analysis (IVA) [12], [13], Nonnegative Matrix
Factorization [14], etc. Some methods separate all of the one-dimensional components of s [15], [16], extract
selected components only [17], or separate multidimensional components; see, e.g., [18], [19], [20], [21], [22].
The separation can also proceed in two steps where a steering vector/matrix (e.g., H1 ) is identified first, while the
signals are separated in the second step using an array processor such as the minimum variance distortion-less
(MVDR) beamformer [23].
The blind separation or identification is often not unique. For example, the order and scaling factors of
the separated components are random and cannot be determined without additional assumptions. Throughout
this paper, we will always assume that the problem of the random order (known as the permutation problem)
has already been resolved, so the subspaces in (2) are correctly identified; for practical methods solving the
permutation problem, see, e.g., [5], [24], [25].
The reconstruction of the signal images is a way to cope with the scaling ambiguity [7]. The advantage
is that si can be retrieved without prior knowledge of the scale of si in (1) or in (2). The scale of si has
clear physical interpretation (e.g., voltage), so the retrieval is highly practical. For example, in the Frequency
Domain ICA for audio source separation, the scaling ambiguity must be resolved within each frequency bin,

3

which is important for reconstructing the spectra of the separated signals in the time domain [5]. Some recent
BSS methods aim to consider the observed signals directly as the sum of images of original sources, by which
the scaling ambiguity is implicitly avoided and the number of unknown free parameters in the BSS model is
decreased [4], [26], [27]. This motivates us for this study, because the way to reconstruct the signal images
(either s1 or s2 ) is an important topic.
In this paper, we study two widely used methods to estimate the source images: One approach performs the
least-squares projection of a separated source on the subspace spanned by X. The other approach assumes that a
blind estimate of a demixing transform is available and exploits its inverse to compute the sources’ images. Both
estimators are invariant to the unknown scaling of the separated sources. The goal of this study is to compare the
estimators, which was not performed yet, and to provide a guidance which estimator is advantageous compared
to the other from different aspects. We also show conditions under which the estimators are equivalent.
The following section introduces the estimators and points to their important properties and relations. Section
III contains a perturbation analysis that studies cases where the estimated demixing transform contains “small”
errors. Section IV studies properties of the least-squares estimator in underdetermined situations, that is, when
there are more original signals than the observed ones. Section V presents results of simulations, and, finally,
Section VI demonstrates two applications.
II. S OURCE I MAGE E STIMATORS
Consider an exact demixing transform as a regular d × d matrix W defined as such that
WH = bdiag(Λ1 , Λ2 ),

(3)

where Λ1 and Λ2 are arbitrary nonsingular matrices representing the random scaling factors of dimensions
m × m and (d − m) × (d − m), respectively; bdiag(·) denotes a block-diagonal matrix with the arguments on
its block-diagonal. By applying W to x, the outputs are


 
y
 =  1 .
y = Wx = WHs = 
Λ2 s2
y2
Λ1 s1



(4)

The components y1 = Λ1 s1 and y2 = Λ2 s2 are separated in the sense that each is a mixture only of s1 and
s2 , respectively.
Let W1 and W2 be sub-matrices of W such that W = [W1 ; W2 ], and W1 contains the first m rows of
W, i.e., y1 = W1 x. From (3) it follows that W is demixing if and only if1 W1 H2 = 0 and W2 H1 = 0.
Throughout the paper we will occasionally use the following assumptions. Consider an estimated demixing
c
matrix W.
c i Hj = 0 where j ∈ {1, 2}, j 6= i. Assuming an exact demixing transform thus
A1(i) The assumption that W
corresponds to A1(1) simultaneously with A1(2).
A2 The assumption of uncorrelatedness of s1 and s2 means that E[s1 sH
2 ] = 0.
1A

more general definition is that W is demixing if and only if W1 H2 s2 = 0 and W2 H1 s1 = 0. However, we will assume that the

mixing model (1) is determined, so cases where W1 H2 s2 = 0 while W1 H2 6= 0 and similar do not exist.

4

A3 The assumption of orthogonality of Y1 and Y2 , that is,
Y1 Y2H /N = 0,

(5)

in BSS also known as the orthogonal constraint [28], means that the sample-based estimate of E[y1 y2H ] is
exactly equal to zero.
In the determined case r = d and under A1(1) and A1(2), the condition (5) corresponds with
S1 SH
2 /N = 0,

(6)

but not generally so when r > d. The latter condition could be seen as a stronger alternative to A2.
For example, the orthogonal constraint A3 is used by some ICA methods such as is Symmetric or Deflation
FastICA [15]. There are several reasons for this. First, A2 is the necessary condition of independence of the
original signals, so A3 is a practical way to decrease the number of unknown parameters in ICA. Second, A3
helps to achieve the global convergence (to find all independent components) and prevents algorithms from
c i is already determined
finding the same component twice. Finally, in the model (2) with the A3 constraint, W
c j is given, j 6= i, and vice versa.
up to a scaling matrix when W
A. Reconstruction Using Inverse of Demixing Matrix
The estimator to retrieve si described here will be abbreviated as INV.
c = [W
c 1; W
c 2 ] denote an estimated demixing matrix by a BSS method, and let
Definition 1 (INV): Let W
4
b be its inverse matrix, i.e., A
b =
c −1 . Let A
b = [A
b1 A
b 2 ] be divided in the same way as the system matrix
A
W

H. Then, the INV estimator is defined as
b iW
c ix
b
siINV = A

or

bi
b c
S
INV = Ai Wi X,

(7)

for i ∈ {1, 2}.
In particular, INV is popular in the frequency-domain ICA for audio source separation; see, e.g., [7], [6],
[33]. The following two propositions point to its important properties.
Proposition 1 (consistency of INV): Consider an exact demixing transform W satisfying A1(1) and A1(2),
that is, satisfying (3); let A = [A1 A2 ] be its inverse matrix. For i ∈ {1, 2}, it holds that
b
siINV = Ai Wi x = si .

(8)

Proof: By (3) it holds that Ai = Hi Λi−1 . Then,
Ai Wi x = Ai Wi (H1 s1 + H2 s2 )

(9)

= Ai Wi Hi si

(10)

i
= Hi Λ−1
i Λi si = Hi si = s .

(11)

c be an estimated demixing transform and A
b =W
c −1 .
Proposition 2 (scaling invariance of INV): Let W
c ← bdiag(Λ1 , Λ2 )W,
c where Λ1 and Λ2 are arbitrary
The INV estimator is invariant to substitution W
nonsingular matrices of dimensions m × m and (d − m) × (d − m), respectively.

5

c 1 ; Λ2 W
c 2 ]−1 = [A
b 1 Λ−1 A
b 2 Λ−1 ].
Proof: The proof follows from the fact that [Λ1 W
1
2
One advantage is that the transform Ai Wi is purely a function of W and does not explicitly depend on the
signals or on their statistics, e.g., on the sample-based covariance matrix. This makes the approach suitable for
real-time processing [29].
On the other hand, Ai Wi is a function of the whole W through the matrix inverse; it does not depend
solely on Wi , as one would expect when only si should be estimated. Formula (7) can thus be used only if
the whole demixing W is available. BSS methods extracting only selected components (e.g., one-unit FastICA
[15]) cannot be applied together with (7). Next, it follows that potential errors in the estimate of W2 can have
an adverse effect on the estimation of s1 . This is analyzed in Section III.
B. Least-squares reconstruction
Another approach to estimate si is to find an optimum projection of the separated components back to the
observed signals x in order to find their contribution to them. A straightforward way is to use the quadratic
criterion, that is, least squares, which gives two estimators that will be abbreviated by LS.
c i denote an estimated part of a demixing matrix, i ∈ {1, 2}, yi = W
c i x and
Definition 2 (LS): Let W
c i X. The theoretical LS estimator of si is defined as
Yi = W



b
siLS = arg min E kx − Vyi k2 x
V

c H (W
c i CW
c H )−1 W
c i x,
= CW
i
i

(12)
(13)

where C = E[xxH ]. The practical LS estimator of Si is defined as


b i = arg min kX − VYi k2 X
S
LS
F

(14)

bW
c H (W
c iC
bW
c H )−1 W
c i X,
=C
i
i

(15)

V

b = XXH /N , and k · kF denotes the Frobenius norm.
where C
Proposition 3 (scaling invariance of LS): The estimators (12) and (14) are invariant to a scaling transform
c i ← Λi W
c i where Λi is a nonsingular square matrix.
W
The proof of Proposition 3 is straightforward. It is worth pointing out that the LS estimators are purely
c i as compared to INV. Also, they involve the covariance matrix or its sample-based estimate.
functions of W
However, their consistency is not guaranteed under A1(i) even if the assumption holds for both i = 1, 2 as
assumed in Proposition 1. In fact, additional assumptions are needed as is shown by the following proposition.
Proposition 4: Let Wi be a part of an exact demixing matrix, so A1(i) holds. Let Wi Hi = Λi be
nonsingular. Then, under A2 it holds that
b
siLS = CWiH (Wi CWiH )−1 Wi x = si ,

(16)

b i = CW
b H (Wi CW
b H )−1 Wi X = Si .
S
LS
i
i

(17)

and under A3 it holds that

Proof: The proof will be given for (16) while the one for (17) is analogous.

6

According to (1) it holds that C = E[xxH ] = H Cs HH where Cs = E[ssH ]. Under A2 it follows that Cs
4

4

H
has the block-diagonal structure Cs = bdiag(Cs1 , Cs2 ) where Cs1 = E[s1 sH
1 ] and Cs2 = E[s2 s2 ] are regular

(because Cs is assumed to be regular). Without a loss of generality, let i = 1. Since W1 H = (Λ1 0),
b
siLS = CW1H (W1 CW1H )−1 W1 x

(18)

= HCs HH W1H (W1 HCs HH W1H )−1 W1 Hs

(19)

H −1
= H1 Cs1 ΛH
Λ1 s1
1 (Λ1 Cs1 Λ1 )

(20)

= H1 s1 = s1 .

(21)

b H . This
It is worth pointing out that LS involves a matrix inverse, namely, of Wi CWiH or of Wi CW
i
matrix (actually, the (sample) covariance of (Yi ) yi ) has a lower dimension than W and is more likely well
conditioned so that the computation of its inverse is numerically stable.
C. On the equivalence between INV and LS under the orthogonal constraint
c 1 is given. Under the assumption A3, W
c 2 is already determined
Without a loss of generality, assume that W
c is actually available, and the INV estimator (7) can be
up to a scaling matrix through (5), so the whole W
applied. The goal here is to verify that, in this case, INV coincides with LS.
c = [W
c 1 ; B]. Then,
Let B denote the unknown lower part of the entire demixing W

  
c X
W
Y
c =  1  =  1 .
WX
BX
Y2

(22)

The condition (5) requires that
bW
c H = 0,
BC
1

(23)

bW
c H . It can be verified that any B of the
which means that the rows of B are orthogonal to the columns of C
1
form
bW
c H (W
c 1C
bC
bW
c H )−1 W
c 1 C),
b
B = Q(I − C
1
1

(24)

where Q can be an arbitrary (d − m) × m full-row-rank matrix such that B has full row-rank, meets the
condition (23).
b 1 must be computed, which consists of first m columns of A
b =W
c −1 , so it satisfies
Now, to apply (7), A
c 1A
b 1 = I,
W

(25)

bW
c H (W
c 1C
bC
bW
c H )−1 W
c 1 C)
b A
b 1 = 0.
BA1 = Q(I − C
1
1

(26)

b1 = C
bW
c H R where R is an m × m matrix. To satisfy also (25),
The latter equation is satisfied whenever A
1
c 1C
bW
c H )−1 . Finally,
R = (W
1
b 1W
c 1X = C
bW
c H (W
c 1C
bW
c H )−1 W
c 1 X,
S1INV = A
1
1

(27)

which allows us to conclude this section by the following proposition.
c i be a part of an estimated demixing matrix, i ∈ {1, 2}, and let A3 be assumed. Then,
Proposition 5: Let W
bi
bi
S
INV = SLS .

7

III. P ERTURBATION A NALYSIS
Throughout this section, let W denote the exact demixing transform, that is W = H−1 . We present an
analysis of the sensor response estimators (7) and (14) when W is known up to a small deviation2 . Let
V = W + Ξ be the available observation of W where Ξ is a “small” matrix;V1 will denote the sub-matrix of
V containing the first m rows; similarly Ξ = [Ξ1 ; Ξ2 ]; A = V−1 and A1 contains the first m columns of A;
b − C be “small” and of the same asymptotic order as Ξ (typically, ∆C = Op (N −1/2 ) where
let also ∆C = C
Op (·) is the stochastic order symbol).
4
4
b H (V1 CV
b H )−1 V1 estimating S1 from
Now, consider the transform matrices T1 = A1 V1 and T2 = CV
1
1

X, respectively. The analysis resides in the computation of their squared distances (the Frobenius norm) from
the ideal transform, that is, from H1 W1 . Using first-order expansions and neglecting higher-order terms, it is
derived in Appendix that the following approximations hold.
2

2

kH1 W1 − T1 kF ≈ kH Ξ H1 W1 − H1 Ξ1 kF ,

(28)

2

−1
kH1 W1 − T2 kF ≈ H1 (Ξ1 CW1H + W1 C ΞH
1 )Cs1 W1
H −1
(H1 W1 − I)∆CW1H C−1
s1 W1 − H1 Ξ1 − C Ξ1 Cs1 W1

2

.
F

(29)

To provide a deeper insight, we will analyze a particular case where H = W = I, Cs1 = σ12 I, and
Cs2 = σ22 I. Let the elements of Ξ all be independent random variables with zero mean such that the variance
of each element of Ξi is equal to λ2i . For further simplification, let the elements of ∆Cij , which denotes the
ijth block of ∆C, i, j = 1, 2, be also independent random variables whose variance is equal to σi σj C. Then,
the expectation values of the right-hand sides of (28) and (29), respectively, are equal to
i
h
2
E kH1 W1 − T1 kF ≈ (λ21 + λ22 )m(d − m),


h
i 
σ2
σ2
2
E kH1 W1 − T2 kF ≈
1 + 22 λ21 + C m(d − m).
σ1
σ1

(30)
(31)

Comparing (30) and (31) shows the pros and cons of the estimators. The latter depends on σ22 /σ12 , which
reflects the ratio between the power of s1 and that of s2 . The expression (30) does not depend on this ratio
explicitly3 . For simplicity, let us assume that σ22 /σ12 = 1.
Next, (31) depends on C while (30) is independent of it. Since C captures the covariance estimation error
in ∆C, it typically decreases with the length of data N . Usually, C has asymptotic order O(N −1/2 ); see, e.g.,
Appendix A.B in [30]. For the sake of the analysis, we will assume that C = 0. Then, (31) changes to
i
h
2
E kH1 W1 − T2 kF ≈ 2λ21 m(d − m).

(32)

The expressions (30) and (32) point to the main difference of the estimators when H = W = I: While the
performance of INV depends on λ21 and λ22 , that of LS depends purely on λ11 . In the special case λ11 = λ12 , the
performances coincide.
2 The

analysis of (12) follows from that of (14) when ∆C = 0.

3 Typically,

σ22 .

there is an implicit dependency of the estimation error Ξ on σ22 /σ12 . Therefore, λ21 as well as λ22 are influenced by σ12 and

8

20
10

MSE [dB]

0

INV - theory
INV - experimental
LS - theory
LS - experimental
LS - experimental, N=106

-10
-20
-30
-40
-50
10-3

10-2

10-1

1
Fig. 1.

Comparison of theoretical and experimental mean square error of INV and LS when d = 20, m = 5, H = W = I, σ12 = σ22 = 1,

C = 0, N = 105 , λ22 = 10−4 .

To verify the theoretical expectations, we conducted a simple simulation where d = 20, m = 5, H = W = I,
σ12

= σ22 = 1, C = 0, N = 105 , λ22 = 10−4 . Ξ and S were drawn from complex Gaussian distribution with

zero mean and unit variance. The average squared errors of INV and LS averaged over 100 trials for each value
of λ1 are compared, respectively, with the expressions (30) and (32) in Fig 1.
The theoretical error of INV is in a good agreement with the experimental one for all values of λ1 . The
same holds for that of LS until λ1 ≥ 0.003. For smaller values of λ1 , the experimental error of LS is limited
while the theoretical one is decreasing. This is caused by the influence of ∆C, which is fully neglected in (32)
and modeled through C in (31). The experimental error of LS evaluated for N = 106 confirms that (31) is a
more accurate theoretical error of LS than (32).
In this example, INV outperforms LS when λ1 > λ2 , and vice versa. However, these results are valid only
in the special case when H = W = I. Simulations in Section V consider general mixing matrices, thereby
compare the estimators in more realistic situations.

9

IV. N OISE E XTRACTION FROM U NDERDETERMINED M IXTURES
A. Mixture Model
Now we focus on a more realistic scenario that appears in most array processing problems. Let the mixture
be described as
x = H1 s1 + s2 ,

(33)

where H1 is a d × m matrix having full column rank, s1 is an m × 1 vector of target components, and s2 is a
d × 1 vector of noise signals. Note that, in this model, s2 is simultaneously equal to s2 .
The mixture model corresponds with (1), but H is equal to [H1 Id×d ] and has dimensions d × (m + d),
which makes the problem underdetermined (r = m + d).
In general, a linear transform that separates s1 or s2 from x does not exist, unless (33) is implicitly regular
(e.g., when Cs2 has rank d − m)4 . From now on, we focus on the difficult case where, generally speaking,
neither s1 nor s2 can be separated.
B. Target Signal Cancelation and Noise Extraction
Since the separation of s1 is not possible, multichannel noise reduction systems follow an inverse approach:
the target components s1 are first linearly canceled from the mixture in order to estimate a reference of the
noise components s2 . Second, a linear transform or adaptive filtering is used to subtract the noise from the
mixture as much as possible; see, e.g., [37], [38], [39], [40], [42], [43].
Specifically, the cancelation of the target component is achieved through a matrix W such that
WH1 = 0.

(34)

Since H1 has rank m, the maximum possible rank of W is d − m, which points to the fundamental limitation:
The maximum dimension of the subspace spanned by the extracted noise signals Wx = Ws2 is d − m.
Assume for now that any (d − m) × d matrix W having full row-rank has been identified (e.g., using BSS).
To estimate s2 , LS can be used (INV cannot be applied in the underdetermined case), so
b
s2LS = CWH (WCWH )−1 Wx,

(35)

b 2 = CW
b H (WCW
b H )−1 WX.
S
LS

(36)

or

Proposition 6: Let W be a (d − m) × d transform matrix having rank d − m and satisfying (34), and let
Q denote a d × (d − m) matrix. Under A2, s2LS is a minimizer of


min E ks2 − b
sk2 .

b
s=QWx

(37)

Assuming (6), S2LS is a minimizer of
min
b
S=QWX

4 For

b 2.
kS2 − Sk
F

(38)

example, model (33) is often studied under the assumption that at most d signals out of s1 and s2 are active at a given time instant;

see, e.g. [35], [36].

10

Proof: Under A2 it holds that




min E ks2 − b
sk2 = min E kx − b
sk2 .

b
s=QWx

b
s=QWx

When (6) holds, than
b 2 =
min kS2 − Sk
F

b
S=QWx

b 2.
min kX − Sk
F

b
S=QWx

The statements of the proposition follow, respectively, by the definitions (12) and (14).
The latter proposition points to important limitations of LS that should be taken into account in the underdetermined scenario. First, LS estimates a d-dimensional signal only from a (d − m)-dimensional signal subspace.
Second, (35) is optimal under A2 in the least-squares sense of (37). Third, (36) is optimal in the sense of (38)
only when (6) is valid, which is much stronger assumption than A2.
V. S IMULATIONS
This section is devoted to extensive Monte Carlo simulations where the signals and system parameters are
randomly generated. Real and complex parts of random numbers are always generated independently according
to the Gaussian law with zero mean and unit variance. Each trial of a simulation consists of the following steps.
1) The dimension parameters d and m are chosen.
2) N samples of the original components s1 and s2 are randomly generated according to the Gaussian law.
b = XXH /N .
3) The mixing matrix H is generated, W = H−1 , X = HS, and C
c 1 = W1 + Ξ1
4) The estimation of W is simulated by adding random perturbations to its blocks, that is, W
c 2 = W2 + Ξ2 , where the elements of Ξ1 and Ξ2 have, respectively, variances λ2 and λ2 ; W
c =
and W
1
2
c1 W
c 2 ]5 . Then, W
c 1 and W
c 2 are multiplied by random regular scaling matrices of corresponding
[W
dimensions.
c is evaluated through the normalized mean5) The accuracy of the INV and LS estimates of S1 using W
squared error defined as
NMSEj =

kH1 W1 − Tj k2F
,
kH1 W1 k2F

(39)

j = 1, 2.
The following subsection reports results of simulations assuming the determined model. The next subsection
considers the underdetermined model (33).
A. Determined model
c The experiment is done with d = 5, m = 2, and N = 104 ; λ2
1) Influence of the Estimation Errors in W:
2
is equal to one of four constants (10−1 , 10−2 , 10−3 , and 10−4 ) while λ21 is varied. Each simulation is repeated
in 105 trials. The average NMSE achieved by INV and LS are shown in Fig. 2.
c 2 . For example, for
The results of INV are highly influenced by λ22 that controls the perturbation of W
λ22 = 10−1 and λ22 = 10−2 , INV fails in the sense that the achieved NMSE is above 0 dB. This happens even
5 Note

c 1 and W
c 2 are not constrained to satisfy A3. Otherwise, the comparison of INV and LS would not give a sense due to
that W

Proposition 5.

11

20
15
10

mean NMSE [dB]

5
0

INV, λ22=10−1

−5

LS, λ2=10−1

−10

INV, λ2=10

−15

LS, λ22=10−2

−20

INV, λ22=10−3

2
2

LS, λ22=10−3

−25

INV, λ22=10−4

−30
−35
−40

Fig. 2.

−2

LS, λ22=10−4
−30

−20

−10

0
10⋅log10(λ21/λ22)

10

20

30

40

[dB]

NMSE averaged over 105 trials where d = 5, m = 2, N = 104 , λ22 is fixed, and λ21 is varied.

c 1 , is relatively “small”. For λ2 = 10−3 and λ2 = 10−4 , the NMSE
if λ21 , which controls the perturbation of W
2
2
of INV decreases with decreasing λ21 . However, the NMSE is lower bounded (does not improve as λ21 → 0).
c 2.
All these results point to the dependency of INV on W
The NMSE of LS depends purely on λ21 . It is always improved with the decreasing value of λ21 (it is
b In this
only limited by the length of data which influences the accuracy of the sample covariance matrix C).
experiment, LS is outperformed by INV only in a few cases, namely, when λ22 = 10−4 and λ21 /λ22 is higher than
c is a sufficiently
−14 dB. INV thus appears to be beneficial compared to LS in situations where the whole W
accurate estimate of W.
2) Varying Dimension: In the situation here, the target component s1 has dimension one, i.e., m = 1, while
the dimension of the mixture d is changed from 2 through 20; N = 104 . The variances λ21 and λ22 are fixed,
namely, λ22 = 10−3 and λ21 is chosen such that 10 log10 λ21 /λ22 corresponds, respectively, to −10, 0, and 10 dB.
The NMSE averaged over 105 trials is shown in Fig. 3.
The NMSE values of both methods are increasing with growing d. In the INV case, the NMSE grows
smoothly until it reaches a certain threshold value of d. The experiments show that this threshold depends on
λ21 and λ22 . Above this threshold, the NMSE of INV abruptly grows. It points to a higher sensitivity of INV to
c when the dimension of data is “high”.
the estimation errors in W
LS yields smooth and monotonic behavior of NMSE for every d. It is outperformed by INV if both λ21 and

12

25
2

20

LS, λ2/λ2 = −10 dB
1 2
INV, λ2/λ2 = 0 dB
1 2
2 2
LS, λ1/λ2 = 0 dB
INV, λ2/λ2 = 10 dB
1 2
LS, λ21/λ22 = 10 dB

15

mean NMSE [dB]

2

INV, λ1/λ2 = −10 dB

10
5
0
−5
−10
−15
−20
−25
2

4

6

8

10

12

14

16

18

20

d
Fig. 3.

NMSE averaged over 105 trials as a function of d = 2, . . . , 20; here m = 1, λ22 = 10−3 , and N = 104 .

λ22 as well as the data dimension d are sufficiently small.
3) Target Component Dimension: The dimension of the mixture d is now put equal to 20, while the dimension
of the target component m is varied from 1 through d − 1. Results for three different choices of λ21 and λ22
are shown in Fig. 4. The scenario with λ21 = λ22 = 10−3 appears to be difficult for both methods as they
do not achieve NMSE below 0 dB. INV also fails when λ22 = 10−3 and λ21 /λ22 corresponds to −10 dB (i.e.,
λ21 = 10−4 ) until m ≤ 17. This is in accordance with the results of the previous example that shows that INV
fails when λ21 , λ22 and d are “too large”. The example here reveals one more detail: INV can benefit from
smaller perturbations of the target component (λ21 = 10−4 ) even if λ22 is larger, but the target dimension must
be large enough with respect to d.
LS performs independently of λ22 , which is confirmed by the cases that are plotted with solid and dashed lines
in Fig. 4: these lines coincide as both correspond to the same λ21 (although different λ22 ). LS is outperformed
c is very small.
by INV when λ21 = λ22 = 10−4 , which, again, occurs when the estimation error of the whole W
B. Underdetermined model
In the example of this subsection, we consider the underdetermined mixture model (33) where m = 1,
d = 2, . . . , 20, and N = 50, . . . , 105 . The goal is to examine the reconstruction of the noise components S2
through (36). H1 is randomly generated. Then, W is such that its rows form a basis of the (d−m)-dimensional
subspace that is orthogonal to H1 plus a random Gaussian perturbation matrix whose elements have the variance

13

20
2

−3

2

2

INV, λ2=10 , λ1/λ2 = −10 dB
LS, λ2=10−3, λ2/λ2 = −10 dB
2
1 2
2
−4 2 2
INV, λ2=10 , λ1/λ2 = 0 dB
LS, λ22=10−4, λ21/λ22 = 0 dB
INV, λ22=10−3, λ21/λ22 = 0 dB
LS, λ22=10−3, λ21/λ22 = 0 dB

mean NMSE [dB]

15

10

5

0

−5

−10

2

4

6

8

10

12

14

16

18

m
Fig. 4.

NMSE averaged over 105 trials where d = 20, N = 104 , and m = 1, . . . , 19.

4
2

λ21 = 10−2

0

λ21 = 10−3
λ21 = 10−4

−2

λ21 = 10−5

2

mean NMSES [dB]

MMSE

λ21 = 10−6

−4
−6
−8
−10
−12
−14
2

4

6

8

10

12

14

16

18

20

d
Fig. 5.

Average NMSES2 as a function of d achieved by (35) in an experiment with the underdetermined model (33); m = 1; N = 104 .

MMSE denotes the NMSE achieved by the optimum minimum mean-squared error solution (41). The signals are generated as random
complex Gaussian i.i.d.

14

2

−3

2

2

mean NMSES [dB]

λ1=10

2

−3

−3.5

−3.5

−4

−4

−4.5

−4.5

−4.5

−5

−5

−5

−5.5

−5.5

−5.5

−6

−6

−6

−6.5

−6.5

−6.5

−7

−7

−7

2

10

3

4

10

10

−7.5

5

10

−5

λ1=10

−3

−7.5

−3
−3.5

LS
LSopt
MMSE

2

10

N
Fig. 6.

−4

λ1=10

3

4

10

10

N

−4

−7.5

5

10

2

10

3

4

10

10

5

10

N

Average NMSES2 as a function of N achieved through (36), (38) and (14); d = 4; λ21 = 10k , k = −3, −4, −5.

values equal to λ21 = 10k , k = −2, −3, . . . , −6, respectively. After applying (36), the evaluation is done using
the normalized mean square distance
NMSES2 =

b 2 k2
kS2 − S
F
.
kS2 k2F

(40)

b 2 with the exact solution of (38), which
Owing to the statement of Proposition 6, it is worth comparing S
will be abbreviated by LSopt, and with the minimum mean square error solution, marked as MMSE, defined
as the minimizer of
min kS2 − QXk2F .

Q∈Cd×d

(41)

The latter gives the minimum achievable value of NMSEs2 by a linear estimator; cf. (38) and (14).
The results averaged over 103 independent trials are shown in Figures 5 and 6. Fig. 5 shows results for
N = 104 . One observation here is that NMSES2 achieved through LS is getting closer to that of MMSE as λ21
approaches zero. Next, NMSES2 improves with growing dimension d, but it appears that it stops improving at
a certain d and grows beyond this threshold value, which depends on λ21 . For example, when λ21 = 10−4 , the
NMSES2 is decaying until d = 8 and grows beyond d ≥ 10.
Fig. 6 shows NMSES2 as a function of N when d = 4. This detailed observation shows that LS approaches
LSopt as N grows and λ21 approaches zero, but does not achieve the performance of MMSE. This is the
fundamental limitation due to the dimension of the separable signal subspace, that is, d − m.

15

VI. P RACTICAL E XAMPLES
A. De-noising of Electrocardiogram
Fig. 7 shows two seconds of a recording from a three channel electrocardiogram (ECG) of a Holter monitor,
which was sampled at 500 Hz. The recording is strongly interfered with a noise signal originating from the
Holter display. The fundamental frequency of the noise is about 37 Hz, and the noise contains several harmonics.
Since the noise is significantly stronger than the ECG components, Principal Component Analysis (PCA)
can be used to find a demixing transform that separates the noise from the mixture. Therefore, we take the
eigenvector corresponding to the highest eigenvalue of the covariance matrix of the recorded data (the principal
vector) as the separating transform. Then, the noise responses on the electrodes are computed using LS and
subtracted from the original noisy recording. This approach is computationally cheaper that doing the whole
PCA and using INV then. According to Proposition 5, both approaches give the same result as PCA yields
component that are exactly orthogonal.
To compare, we repeated the same experiment using the vector obtained through Independent Component
Analysis (ICA). One-unit FastICA [15] with tanh(·) nonlinearity was used to compute the vector separating
the noise component. To avoid the permutation ambiguity, the algorithm was initialized from [1 1 1], because
the noise appears to be uniformly distributed over the electrodes. Also here the approach is faster than doing
the whole orthogonally-constrained ICA (e.g. using Symmetric FastICA) and using INV.
Figures 8 and 9 show the resulting signals where the estimated images of the noise component were removed,
respectively, through PCA and ICA. Both results show very efficient subtraction of the noise. A visual inspection
of the detail in Fig. 8 shows certain residual noise that does not appear in Fig. 9, so the separation through
ICA appears to be more accurate than by PCA. Combining one-unit ICA algorithm with LS, the computational
complexity of the ICA solution is decreased.
B. Blind Separation with Incomplete Demixing Transforms
Proposition 5 points to the fact that if BSS is based on a method that yields (almost) orthogonal estimates of
s1 and s2 , INV and LS are principally not that different. The example of this section demonstrates a situation
where the estimates are significantly nonorthogonal, so INV and LS yield considerably different results.
Recently, a novel approach for blind separation of convolutive mixtures of audio signals has been proposed
in [45], [47]. The idea resides in the application of IVA in the frequency-domain on a constrained subset of
frequencies where the input (mixed) signals are active. This is in contrast with the conventional FrequencyDomain IVA (ICA), which is applied in all frequencies. The motivation behind is threefold: computational
savings, improved accuracy (especially in environments with sparse room impulse responses), and the prediction
of complete demixing transform for separation of future signals whose activity appears in other frequencies.
The method proceeds in three main steps. First, the subset S of p percents of the most active frequency bins
is selected. This can be done through estimating the power spectrum of the signal on a reference microphone
using the coefficients of its short-term Fourier transform. The frequency bins with maximum average magnitude
of Fourier coefficients are selected. Second, an IVA method is applied that estimates the demixing matrices

16

# electrode

3

2

1

0
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time [s]
Fig. 7.

A two-second sample of a three channel electrocardiogram interfered by a noise signal originating from a Holter display.

3.5
3

# electrode

2.5
2
1.5
1
0.5
0
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time [s]
Fig. 8.
LS.

Cleaned data from Fig. 7 after the subtraction of noise responses that were estimated through the main principal component and

17

# electrode

3

2

1

0
0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time [s]

Signal−to−Interference Ratio Improvement [dB]

Fig. 9.

Cleaned data from Fig. 7 after the subtraction of noise images that were estimated using the one-unit FastICA and LS.

35
30
25
20
15

LS: IVA + LASSO
LS: oracle + LASSO
LS: oracle + TDOA
INV: IVA + LASSO
INV: oracle + LASSO
INV: oracle + TDOA

10
5
0
0

10

20

30

40

50

60

70

80

90

100

percentage p [%]
Fig. 10.

Improvement of signal-to-interference ratio as a function of p, i.e., of percents of selected active frequencies in S for the

estimation of incomplete demixing transform. The evaluation was performed on speech signals.

Signal−to−Interference Ratio Improvement [dB]

18

25

20

15

LS: IVA + LASSO
LS: oracle + LASSO
LS: oracle + TDOA
INV: IVA + LASSO
INV: oracle + LASSO
INV: oracle + TDOA

10

5

0
0

10

20

30

40

50

60

70

80

90

100

percentage p [%]
Fig. 11.

Improvement of signal-to-interference ratio as a function of p. The evaluation was performed on white noise signals that

uniformly excite the whole frequency range.

within the subset of the selected frequencies. The subset of the matrices is referred to as Incomplete Demixing
Transform (IDT). Third, the IDT is completed by a given method.
We consider the same experiment with two simultaneously speaking persons and two microphones as in [45].
Signals have 10 seconds in length; the sampling frequency is 16 kHz. The signals are convolved with room
impulse responses (RIR) generated by a simulator and mixed together. Reflection order is set to 1 so that the
RIRs are significantly sparse (the results of this experiment with reflection order 10 when RIRs are no more
sparse are available in [45]). Then, the signals are transformed into the short-term Fourier domain with the
window length of 1024 samples and shift 128. The convolution is, in the frequency domain, approximated by
the set of multiplicative models (1) where d = r = 2 where one model corresponds to one frequency bin; there
are 513 models in total.
As for the second step, the demixing matrices are estimated from the mixed signals using the natural gradient
algorithm for IVA [25] applied to the subset of models (1). To compare, “oracle” demixing matrices are derived
on S using known responses of the speakers. This gives the IDT that is known only on the selected subset S.
The IDT is completed by two alternative methods. The first method, denoted as TDOA, utilizes known timedifferences of arrival of the signals. The unknown demixing matrices are such that their rows correspond to the
null beamformer steering spatial null towards the unwanted speaker. The second approach, denoted as LASSO,
completes the IDT through finding the sparsest representations of incomplete relative transfer functions (RTF)

19

that are derived from the IDT6 . Let q denote an |S| × 1 vector that collects the coefficients of an incomplete
RTF; |S| is the number of elements in S. The completed RTF is obtained as the solution of [46]
arg min khS − qk2 + kFH hk1 ,
h

(42)

where  > 0 controls the time-domain sparsity of the solution, F is the matrix of the DFT, the subscript (·)S
denotes a vector/matrix with selected elements/rows whose indices are in S, and k · k1 denotes the `1 -norm.
Now, it is worth noting that the separated components by the demixing matrices after the completion can be
significantly nonorthogonal. While the IVA applied within S aims to find independent (thus “almost” or fully
orthogonal) components, the method for the IDT completion does not take any regard to the orthogonality7 .
Figures 10 and 11 show results of the experiment from [45] evaluated in terms of the Signal-to-Interference
Ratio Improvement (SIR) after the signals are separated as a function of p (the percentage of frequencies in
S). In Fig. 10, the evaluation is performed with the speech signals, while Fig. 11 shows the results achieved as
if the sources were white Gaussian sequences. The purpose of the latter evaluation is to evaluate the completed
IDT uniformly over the whole frequency range, i.e., also in frequencies that were not excited by the speech
signals. Note that SIR must be evaluated after resolving the scaling ambiguity in each frequency [41]. This
gives us the opportunity to apply either INV or LS.
The results in Figures 10 and 11 point to significant differences between LS and INV in this evaluation. The
results by LS appear to be less biased and stable as compared to those by INV, and can be interpreted in accord
with the theory. In particular, LS shows that oracle+LASSO (oracle IDT completed by LASSO) outperforms
oracle+TDOA for p between 35% and 80%. This gives sense, because LASSO can better exploit the sparsity of
the RIRs generated in this experiment. The results by INV do not reveal this important fact. Next, LS shows in
Fig. 10 that IVA+LASSO can improve the separation of the speech signals when p < 100%. The evaluation on
white noise in Fig. 11 shows that the loss of SIR is not essential until p < 30%. The latter conclusion cannot
be drawn with the results by INV.
VII. C ONCLUSIONS
We have analyzed and compared two estimators of sensor images (responses) of sources that were separated
from a multichannel mixture up to an unknown scaling factor: INV and LS. Simulations and perturbation analysis
have shown pros and cons of the methods, which can be summarized into the following recommendations.
•

LS is more practical in a sense that the whole mixing matrix need not be identified for its use, which is
useful especially in underdetermined scenarios.

•

The advantage of INV resides in the independence on the (estimated) covariance matrix.

•

INV could be beneficial as compared to LS when used with non-orthogonal BSS algorithms, i.e., those not
applying the orthogonal constraint. However, both the target as well as the interference subspaces must be
estimated with a sufficient accuracy.

6 As

pointed in [45], LASSO could be seen as a generalization of TDOA, because impulse responses corresponding to null beamformers

are pure-delay filters, which are perfectly sparse.
7 The

there.

orthogonal constraint cannot be imposed within the frequencies outside of the set S, because signals are not (or purely) active

20

Both approaches have been shown to be equivalent under the orthogonal constraint, so the differences in their
accuracies are less significant when BSS yields signal components that are (almost) orthogonal (e.g. PCA, ICA,
IVA). By contrast, the differences between the reconstructed images of nonorthogonal components can be large,
as demonstrated in the example of Section VI-B.
A PPENDIX : A SYMPTOTIC E XPANSIONS
Computation of (28)
Let E contain first m columns of the d × d identity matrix. It follows that
HE = H1 ,

AE = A1 ,

EH W = W1 ,

EH V = V1 .

To derive an approximate expression for A, we will use the first-order expansion
A = V−1 = (W + Ξ)−1 = (I + HΞ)−1 H

(43)

≈ (I − HΞ)H = H − HΞH.

(44)

Now we apply this approximation and neglect terms of higher than the first order.
2

kH1 W1 − A1 V1 kF = H1 W1 − AEEH V

2
F

≈

H1 W1 − (H − HΞH)EEH (W + Ξ)

2
F

=

2

kH1 W1 − (H1 − HΞH1 )(W1 + Ξ1 )kF ≈
2

kHΞH1 W1 − H1 Ξ1 kF . (45)

Computation of (29)
We start with the first approximation
b H (V1 CV
b H )−1 V1
H1 W1 − CV
1
1

2

=
F

H1 W1 − (C + ∆C)(W1H + ΞH
1 )
−1
· ((W1 + Ξ1 )(C + ∆C)(W1H + ΞH
(W1 + Ξ1 )
1 ))

2

≈
F

H
H1 W1 − (C + ∆C)(W1H + ΞH
1 ) · (W1 CW1 +
−1
W1 ∆CW1H + Ξ1 CW1H + W1 CΞH
(W1 + Ξ1 )
1 )

2
F

. (46)

21

Since W is now the exact inverse of H, it holds that W1 CW1H = Cs1 . By neglecting higher than the first-order
terms and by applying the first-order expansion of the matrix inverse inside the expression,
−1
H
H1 W1 − (C + ∆C)(W1H + ΞH
1 )(I + Cs1 Ξ1 CW1 +
H −1 −1
H
−1
+ C−1
Cs1 (W1 + Ξ1 )
s1 W1 ∆CW1 + Cs1 W1 CΞ1 )

2

≈
F

H
−1
H
H1 W1 − (CW1H + CΞH
1 + ∆CW1 )(I − Cs1 Ξ1 CW1 −
H
H
−1
−1
− C−1
s1 W1 ∆CW1 − Cs1 W1 CΞ1 )Cs1 (W1 + Ξ1 )

2
F

. (47)

Since,
H
H −1
CW1H C−1
s1 = H bdiag(Cs1 , Cs2 )H W1 Cs1

= H1 ,

(48)

the zero order term in (47) vanishes. By neglecting higher than the first-order terms, (29) follows.
ACKNOWLEDGMENTS
This work was supported by The Czech Science Foundation through Project No. 14-11898S and partly by
California Community Foundation through Project No. DA-15-114599.
We thank BTL Medical Technologies CZ for providing us the three-channel ECG recording.
R EFERENCES
[1] H. L. Van Trees, Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory, John Wiley & Sons, Inc.,
2002.
[2] P. Comon and C. Jutten, Handbook of Blind Source Separation: Independent Component Analysis and Applications, Academic Press,
2010.
[3] J.-F. Cardoso, “Blind signal separation: statistical principles”, Proceedings of the IEEE, vol. 90, n. 8, pp. 2009-2026, October 1998.
[4] N. Q. K. Duong, E. Vincent, and R. Gribonval, Underdetermined reverberant audio source separation using a full-rank spatial
covariance model, IEEE Transactions on Audio, Speech, and Language Processing, vol. 18, no. 7, pp. 1830-1840, Sept 2010.
[5] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust and precise method for solving the permutation problem of frequency-domain
blind source separation,” IEEE Trans. Speech Audio Processing, vol. 12, no. 5, pp. 530-538, 2004.
[6] S. Ukai, T. Takatani, H. Saruwatari, K. Shikano, R. Mukai, H. Sawada, “Multistage SIMO-based Blind Source Separation Combining
Frequency-Domain ICA and Time-Domain ICA, IEICE Trans. Fundam., E88-A, no. 3, pp. 642–650, 2005.
[7] K. Matsuoka and S. Nakashima, “Minimal distortion principle for blind source separation,” Proceedings of 3rd International
Conference on Independent Component Analysis and Blind Source Separation (ICA ’01), pp. 722-727, San Diego, Calif, USA,
Dec. 2001.
[8] S. Markovich, S. Gannot and I. Cohen, Multichannel Eigenspace Beamforming in a Reverberant Noisy Environment with Multiple
Interfering Speech Signals,” IEEE Transactions on Audio, Speech and Language Processing, vol. 17, no. 6, pp. 1071–1086, Aug.
2009.
[9] S. Sanei and J. A. Chambers, EEG Signal Processing, Wiley, July 2007.
[10] L. De Lathauwer, B. De Moor, J. Vandewalle, “Fetal Electrocardiogram Extraction by Blind Source Subspace Separation”, IEEE
Transactions on Biomedical Engineering, Special Topic Section on Advances in Statistical Signal Processing for Biomedicine, vol.
47, no. 5, pp. 567–572, May 2000.
[11] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, Wiley-Interscience, New York, 2001.
[12] T. Kim, I. Lee; T.-W. Lee, “Independent Vector Analysis: Definition and Algorithms,” The Fortieth Asilomar Conference on Signals,
Systems and Computers, pp. 1393–1396, 2006.

22

[13] M. Anderson, G. Fu, R. Phlypo, T. Adali, “Independent Vector Analysis: Identification Conditions and Performance Bounds,” IEEE
Transactions on Signal Processing, vol. 62, no. 17, pp. 4399–4410, 2014.
[14] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization,” Nature, vol. 401, pp. 788–791, Oct.
1999.
[15] A. Hyvärinen, “Fast and Robust Fixed-Point Algorithms for Independent Component Analysis”, IEEE Transactions on Neural
Networks, vol. 10, no. 3, pp. 626–634, 1999.
[16] J.-F. Cardoso and A. Souloumiac, “Blind Beamforming from non-Gaussian Signals,” IEE Proc.-F, vol. 140, no. 6, pp. 362-370, Dec.
1993.
[17] S. A. Cruces-Alvarez, A. Cichocki, and S. Amari, “From Blind Signal Extraction to Blind Instantaneous Signal Separation: Criteria,
Algorithms, and Stability,” IEEE Transactions on Neural Networks, vol. 15, no. 4, July 2004.
[18] D. Lahat, J.-F. Cardoso, and H. Messer, “Blind Separation of Multidimensional Components via Subspace Decomposition: Performance
Analysis,” IEEE Transactions on Signal Processing, vol. 62, no. 11, pp. 2894–2905, June 2014.
[19] A. Hyvärinen and U. Köster, “FastISA: A fast fixed-point algorithm for independent subspace analysis,” In Proc. European Symposium
on Artificial Neural Networks, Bruges, Belgium, 2006.
[20] J. Cardoso, “Multidimensional independent component analysis,” in Proceedings of the 1998 IEEE International Conference on
Acoustics, Speech and Signal Processing, vol. 4, pp. 1941–1944, 12-15 May 1998.
[21] D. Lahat, J.-F. Cardoso, and H. Messer, “Second-order multidimensional ICA: performance analysis,” IEEE Transactions on Signal
Processing, vol. 60, no. 9, pp. 4598–4610, Sept. 2012.
[22] Z. Koldovský and P. Tichavský, “Time-domain blind separation of audio sources on the basis of a complete ICA decomposition of
an observation space”, IEEE Trans. on Speech, Audio and Language Processing, vol. 19, no. 2, pp. 406–416, Feb. 2011.
[23] J. Capon, “High-resolution frequency-wavenumber spectrum analysis,” Proc. of IEEE, vol. 57, no. 8, pp. 1408–1418, Aug. 1969.
[24] P. Tichavský and Z. Koldovský, “Optimal Pairing of Signal Components Separated by Blind Techniques”, IEEE Signal Processing
Letters, vol. 11, no. 2, pp. 119–122, 2004.
[25] T. Kim, H. T. Attias, S.-Y. Lee, T-W. Lee, “Blind Source Separation Exploiting Higher-Order Frequency Dependencies,” IEEE
Transactions on Audio, Speech, and Language Processing, vol. 15, no. 1, Jan. 2007.
[26] J. F. Cardoso, M. Le Jeune, J. Delabrouille, M. Betoule and G. Patanchon, “Component Separation With Flexible ModelsApplication
to Multichannel Astrophysical Observations,” IEEE Journal of Selected Topics in Signal Processing, vol. 2, no. 5, pp. 735–746, Oct.
2008.
[27] The Planck collaboration, “Planck 2013 results. XII. Diffuse component separation,” Astronomy and Astrophysics, vol. 571, A12,
Nov. 2014,
[28] J.-F. Cardoso, “On the performance of orthogonal source separation algorithms”, Proc. EUSIPCO, pp. 776-779, Edinburgh, September
1994.
[29] K. Matsuoka, “Elimination of filtering indeterminacy in blind source separation,” Neurocomputing, vol. 71, pp. 2113-2126, 2008.
[30] P. Tichavský, Z. Koldovský, and E. Oja, “Performance Analysis of the FastICA Algorithm and Cramér-Rao Bounds for Linear
Independent Component Analysis”, IEEE Trans. on Signal Processing, Vol. 54, No. 4, April 2006.
[31] S. Makino, Te-Won Lee, and H. Sawada, Blind Speech Separation, Springer, Sept. 2007.
[32] Parra, L., and Spence, C.: “Convolutive Blind Separation of Non-Stationary Sources”, IEEE Trans. on Speech and Audio Processing,
Vol. 8, No. 3, pp. 320-327, May 2000.
[33] F. Nesta and M. Matassoni, “Blind source extraction for robust speech recognition in multisource noisy environments,” Journal
Computer Speech and Language, vol. 27, no. 3, pp. 730–725, May 2013.
[34] J. Benesty, J. Chen, and E. Habets, Speech Enhancement in the STFT Domain, Springer Briefs in Electrical and Computer Engineering,
2011.
[35] F. Abrard, Y. Deville, “A time-frequency blind signal separation method applicable to underdetermined mixtures of dependent sources,”
Signal Processing, vol. 85, issue 7, pp. 1389–1403, July 2005.
[36] T-W. Lee, M. S. Lewicki, M. Girolami, T. J. Sejnowski, “Blind Source Separation of More Sources Than Mixtures Using Overcomplete
Representations,” IEEE Signal Processing Letters, vol. 6, no. 4, 1999.
[37] L. Griffiths and C. Jim, “An alternative approach to linearly constrained adaptive beamforming,” IEEE Trans. Antennas Propag., vol.
30, no. 1, pp. 27–34, Jan. 1982.
[38] J. Even, C. Ishi, H. Saruwatari, and N. Hagita, “Close speaker cancellation for suppression of non-stationary background noise

23

for hands-free speech interface,” Proc. of the 11th Annual Conference of the International Speech Communication Association
(Interspeech 2010), pp. 977-980, Makuhari, Chiba, Japan, September 26-30, 2010.
[39] O. Hoshuyama, A. Sugiyama, A. Hirano, “A robust adaptive beamformer for microphone arrays with a blocking matrix using
constrained adaptive filters,” IEEE Transactions on Signal Processing, vol. 47, no. 10, pp. 2677–2684, Oct. 1999.
[40] Y. Takahashi, T. Takatani, K. Osako, H. Saruwatari, K. Shikano, “Blind Spatial Subtraction Array for Speech Enhancement in Noisy
Environment,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 17, no. 4, pp. 650–664, May 2009.
[41] Vincent, E., Gribonval, R., and C. Févotte, “Performance Measurement in Blind Audio Source Separation,” IEEE Trans. on Speech
and Audio Processing, Vol 14, No. 4, pp. 1462–1469, July 2006.
[42] S. Gannot, D. Burshtein, and E. Weinstein, “Signal enhancement using beamforming and nonstationarity with applications to speech,”
IEEE Trans. on Signal Processing, vol. 49, no. 8, pp. 1614–1626, Aug. 2001.
[43] S. Gannot and I. Cohen, “Speech Enhancement Based on the General Transfer Function GSC and Postfiltering,” IEEE Trans. on
Speech and Audio Processing, vol. 12, No. 6, pp. 561–571, Nov. 2004.
[44] Z. Koldovský, J. Málek, and S. Gannot, “Spatial Source Subtraction Based on Incomplete Measurements of Relative Transfer Function,”
IEEE/ACM Trans. on Speech, Audio and Language Processing, vol. 23, no. 8, pp. 1335–1347, Aug. 2015.
[45] Z. Koldovský, F. Nesta, P. Tichavský, and N. Ono, “Frequency-Domain Blind Speech Separation Using Incomplete De-Mixing
Transform,” The 24th European Signal Processing Conference (EUSIPCO 2016), pp. 1663–1667, Budapest, Hungary, Sept. 2016.
[46] R. Tibshirani, “Regression shrinkage and selection via the lasso,” J. Royal. Statist. Soc. B., vol. 58, pp. 267–288, 1996.
[47] J. Janský, Z. Koldovský and N. Ono, “A Computationally Cheaper Method for Blind Speech Separation Based On AuxIVA and
Incomplete Demixing Transform,” The 15th International Workshop on Acoustic Signal Enhancement (IWAENC)‘, Xian, China, Sept.
2016.

Zbyněk Koldovský (M’04-SM’15) was born in Jablonec nad Nisou, Czech Republic, in 1979. He received the M.S.
degree and Ph.D. degree in mathematical modeling from Faculty of Nuclear Sciences and Physical Engineering
at the Czech Technical University in Prague in 2002 and 2006, respectively. He was also with the Institute of
Information Theory and Automation of the Academy of Sciences of the Czech Republic from 2002 to 2016.
Currently, he is an associate professor at the Institute of Information Technology and Electronics, Technical
University of Liberec, and the leader of Acoustic Signal Analysis and Processing (A.S.A.P.) Group. He is the Vicedean for Science, Research and Doctoral Studies at the Faculty of Mechatronics, Informatics and Interdisciplinary
Studies. His main research interests are focused on audio signal processing, blind source separation, independent component analysis, and
sparse representations.
Zbyněk Koldovský has served as a general co-chair of the 12th Conference on Latent Variable Analysis and Signal Separation (LVA/ICA
2015) in Liberec, Czech Republic.

24

Francesco Nesta received the Laurea degree in computer engineering from Politecnico di Bari, Bari, Italy, in
September 2005 and the Ph.D. degree in information and communication technology from University of Trento,
Trento, Italy, in April 2010, with research on blind source separation and localization in adverse environments.
He has been conducting his research at Bruno Kessler Foundation IRST, Povo di Trento, from 2006 to 2012. He
was a Visiting Researcher from September 2008 to April 2009 with the Center for Signal and Image Processing
Department, Georgia Institute of Technology, Atlanta. His major interests include statistical signal processing,
blind source separation, speech enhancement, adaptive filtering, acoustic echo cancellation, semi-blind source
separation and multiple acoustic source localization. He is currently working at Conexant System, Irvine (CA, USA) on the development
of audio enhancement algorithms for far-field applications.
Dr. Nesta serves as a reviewer for several journals such as the IEEE T RANSACTION ON AUDIO , S PEECH , AND L ANGUAGE P ROCESSING,
Elsevier Signal Processing Journal, Elsevier Computer Speech and Language, and in several conferences and workshops in the field of
acoustic signal processing. He has served as Co-Chair in the third community-based Signal Separation Evaluation Campaign (SiSEC 2011)
and as organizer of the 2nd CHIME challenge.

