ARTICLE IN PRESS

Signal Processing 87 (2007) 1833–1847
www.elsevier.com/locate/sigpro

Underdetermined blind sparse source separation for
arbitrarily arranged multiple sensors
Shoko Arakia,b,, Hiroshi Sawadaa, Ryo Mukaia, Shoji Makinoa,ba

NTT Communication Science Laboratories, NTT Corporation, 2– 4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japanb
Graduate School of Information Science and Technology, Hokkaido University, Kita 14, Nishi 9,
Kita-ku, Sapporo-shi, Hokkaido 060-0814, Japan
Received 31 July 2006; received in revised form 7 February 2007; accepted 9 February 2007
Available online 1 March 2007

Abstract
This paper presents a new method for blind sparse source separation. Some sparse source separation methods, which
rely on source sparseness and an anechoic mixing model, have already been proposed. These methods utilize level ratios
and phase differences between sensor observations as their features, and they separate signals by classifying them.
However, some of the features cannot form clusters with a well-known clustering algorithm, e.g., the k-means. Moreover,most previous methods utilize a linear sensor array (or only two sensors), and therefore they cannot separate symmetrically
positioned sources. To overcome such problems, we propose a new feature that can be clustered by the k-means algorithm
and that can be easily applied to more than three sensors arranged non-linearly. We have obtained promising results fortwo- and three-dimensionally distributed speech separation with non-linear/non-uniform sensor arrays in a real room evenin underdetermined situations. We also investigate the way in which the performance of such methods is affected by room
reverberation, which may cause the sparseness and anechoic assumptions to collapse.
r 2007 Elsevier B.V. All rights reserved.
Keywords: Blind source separation; Sparseness; Clustering; Normalization; Binary mask; Speech separation; Reverberation

1. Introduction
Blind source separation (BSS) [1] is an approachfor estimating source signals that uses only themixed signal information observed at each sensor.
The BSS technique for speech dealt with in thispaper has many applications including hands-freeteleconference systems and automatic conference
minute generators.
Corresponding author. NTT Communication Science
Laboratories, NTT Corporation, 2–4 Hikaridai, Seika-cho,
Soraku-gun, Kyoto 619-0237, Japan. Tel.: +81 774 93 5319;
fax: +81 774 93 5158.
E-mail address: shoko@cslab.kecl.ntt.co.jp (S. Araki).

Two approaches have been widely studied andemployed to solve the BSS problem; one is based on
independent component analysis (ICA) (e.g., [2])and the other relies on the sparseness of source
signals (e.g., [3]). Recently, many ICA methods havebeen proposed even for the convolutive BSS
problem [2,4–10]. ICA works well even in areverberant condition when the number of sources
N is less than or equal to the number of sensors M.
On the other hand, the sparseness-based approaches
are attractive because they can handle the underdetermined problem, i.e., N4M.
The sparseness-based approaches can be divided
into two main categories. One method is based on

0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved.
doi:10.1016/j.sigpro.2007.02.003

ARTICLE IN PRESS
1834

S. Araki et al. / Signal Processing 87 (2007) 1833–1847

MAP estimation, where the sources are estimatedafter mixing matrix estimation [11–17], and theother extracts each signals with time-frequency
binary masks [3,18–20]. The former method includes mixing matrix estimation and l 1 -norm minimization in the frequency domain (i.e., for complex
numbers), both of which still present difﬁculties [16].
The latter, the binary mask approach, has theadvantage of being implemented in real time [21]. In
this paper we focus on the binary mask approach.
In the binary mask approach, we assume thatsignals are sufﬁciently sparse, and therefore, we canassume that at most one source is dominant at eachtime–frequency slot. If this assumption holds, ahistogram of the level and frequency normalized
phase differences between two sensor observations
has N clusters [3,18,20]. Because an individual
cluster in the histogram corresponds to an individual source, we can separate each signal by selecting
the observation signal at time–frequency points ineach cluster with a binary mask. The best-known
approach may be the Degenerate Unmixing Estimation Technique (DUET) [3,18,21].
Previously, such clustering was performed manually [3,18], by using kernel density estimation [20],or with an ML-based gradient method [21]. On theother hand, if clustering could be performed with awell-known algorithm such as the k-means clustering or hierarchical clustering [22], the clustering will
be automated and simpliﬁed. To employ a widely
utilized clustering algorithm such as the k-means,
we should be careful about the variances of multiple
variables, in this case the level ratios and phasedifferences. However, frequency normalization of
the phase difference, which is important in terms of
avoiding the permutation problem among frequencies [16,17], sometimes makes the phase differencemuch smaller than the level ratio as shown in
Section 3.2. Such different variances between thefeatures make clustering with the k-means difﬁcult.
This is the prime motivation for this work.
Our second motive is to employ more than three
sensors arranged two- or three-dimensionally, which
could have a non-linear/non-uniform alignment.
Only a few authors have generalized [16,17,23] amethod for more than two sensors. Authors of [23]used up to eight sensors, however, their sensors werestill linearly arranged. The paper [24] has alreadytried a multichannel DUET (DESPRIT) by combining the sparse assumption and the Estimation of
Signal Parameters via Rotational Invariance Technique (ESPRIT); however, their method still limits

the array shape: a linear array or two sets ofcongruent arrays. A two-sensor system and a linear
sensor array limits the separation ability on a twodimensional half-plane, e.g., the previous methods
cannot separate sources placed in a mirror image
arrangement. To allow the free location of sources,
we need more than three sensors arranged two- or
three-dimensionally.
Based on these two motivations, we propose a
new binary mask approach MENUET (MultiplesENsor dUET), which employs the well-known
k-means clustering algorithm. As a feature, ourmethod utilizes the level ratios and phase differences
between multiple observations. To realize level ratio
and phase difference variances of a comparable
level, we propose a way of weighting the phase term
for successful clustering. Moreover, our proposed
method does not require sensor location information. This allows us to employ freely arranged
multiple sensors easily. Therefore, the proposed
method can separate signals that are distributed
two- or three-dimensionally. Our previous paper,[16], utilized a two-dimensional sensor array to test
the MAP approach proposed in [16]. However, thatwork did not employ the frequency normalization,and therefore, suffered from the abovementioned
permutation problem. On the other hand, in thispaper, we employ appropriate frequency normalization for the k-means algorithm. Moreover, we
also apply our proposed method to a threedimensional sensor array, and describe the result.
An additional contribution of this paper is that it
undertakes an investigation of the separation
performance in real world acoustic environments.
Both our proposed method and previous methods
employ assumptions of source sparseness andanechoic mixing (i.e., a simple attenuation anddelay model for a room impulse response). Suchassumptions can easily be affected by reverberation.
We show how the performance is affected when the
problem does not satisfy the assumed conditions.
The organization of this paper is as follows.
Section 2 presents the basic framework of the binary
mask-based BSS method. In Section 3, we describe
some features for clustering, and test how eachfeature will be clustered by the k-means clustering
algorithm. In Section 4, we propose a novel method
MENUET, which includes the estimation of geometric features from multiple sensor observations.
Our proposed feature is suitable for k-meansclustering. Section 5 reports some experimental
results obtained with non-linearly arranged sensors

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

1835

2. Binary mask approach to BSS

where f 2 f0; ð1=TÞf s ; . . . ; ððT  1Þ=TÞf s g is a frequency, winðrÞ is a window that tapers smoothly tozero at each end, t is a new index representing time,and S is the window shift size. As the windowwinðrÞ, in this paper, we utilized a Hanning window
1
2 ð1  cosð2pr=TÞÞ (r ¼ 0; . . . ; T  1).
There are two advantages to working in thetime–frequency domain. First, convolutive mixtures(1) can be approximated as instantaneous mixturesat each frequency:

2.1. Problem description

xj ðf ; tÞ in underdetermined scenarios. Even when thesources and sensors were distributed two- or threedimensionally, we obtained good separation resultswith the k-means algorithm for each scenario under
weak reverberant ðRT60 ¼ 128 msÞ conditions. We
also investigated the performance under more
reverberant conditions ðRT60 ¼ 300 msÞ. The ﬁnal
section concludes this paper.

NX

hjk ðf Þsk ðf ; tÞ

(3)

k¼1

Suppose that sources s1 ; . . . ; sN are convolutively
mixed and observed at M sensors
xj ðtÞ ¼

N X
X
k¼1

hjk ðlÞsk ðt  lÞ;

j ¼ 1; . . . ; M,or in a vector notation,
xðf ; tÞ 

(1)

l

(4)

where hjk ðf Þ is the frequency response from source k
to sensor j, and sk ðf ; tÞ is a frequency–domain timeseries signal of sk ðtÞ obtained by the same operationas (2), x ¼ ½x1 ; . . . ; xM T , and hk ¼ ½h1k ; . . . ; hMk T isa mixing vector that consists of the frequency
responses from source sk to all sensors. The second
advantage is that the sparseness of a source signal
becomes prominent in the time–frequency domain[12,19], if the source is colored and non-stationary
such as speech. The possibility of sk ðf ; tÞ being close
to zero is much higher than that of sk ðtÞ. When thesignals are sufﬁciently sparse in the time–frequency
domain, we can assume that the sources rarely
overlap and, (3) and (4), respectively, can beapproximated as

2.2. Separation procedures
Step 1. Signal transformation to the time– frequency domain: Fig. 1 shows the ﬂow of the binary
mask approach. The binary mask approach usuallyemploys a time–frequency domain representation.
First, time–domain signals xj ðtÞ sampled at frequency f s are converted into frequency–domain
time-series signals xj ðf ; tÞ with a T-point short-time
Fourier transform (STFT):
T=21
X

hk ðf Þsk ðf ; tÞ,

k¼1

where hjk ðlÞ represents the impulse response fromsource k to sensor j. In this paper, we focus
particularly on a situation where the number ofsources N can exceed the number of sensors M(N4M). We assume that N and M are known, andthat the sensor spacing is small enough to avoid thespatial aliasing problem. The goal is to obtainseparated signals yk ðtÞ that are estimations of sk
solely from M observations.

xj ðf ; tÞNX

xj ðr þ tSÞ winðrÞe|2pfr ,xðf ; tÞ  hk ðf Þsk ðf ; tÞ;

STFT

(2)

x(f,t)

yk(f,t)Featureextraction

Mk(f,t)

ISTFT

xM(t)

Θ(f,t)
Clustering

9

k 2 f1; . . . ; Ng,

ð5Þ

k 2 f1; . . . ; Ng,

ð6Þ

where sk ðf ; tÞ is a dominant source at the time–
frequency point ðf ; tÞ. For instance this is approximately true for speech signals [3,15]. Fig. 2(a) shows
example spectra of three speech sources, in which
we can see their temporal/frequency sparseness.

r¼T=2x1(t)9xj ðf ; tÞ  hjk ðf Þsk ðf ; tÞ;

Ck

Mask
designFig. 1. Basic scheme of binary mask approach.

yk(t)ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

1836

Sources
2.0
S1

Observations
2.0
50

100

150

200 X1

2.0
1.0
0.5
50

100

150

200 X2Frequency [kHz]Frequency [kHz]

S2

1.0
0.5

2.0
S3

1.0
0.5
50

50

100
150
Time (frame index)

Y1

1.0
0.5
50

100

150

200

2.0
1.0
0.5

Y2
50

100

150

1.0
0.5
Frequency [kHz]Frequency [kHz]

200

2.0

200

2.0
M3

100
150
Time (frame index)

Separated signals

Masks

M2

200

200

2.0
M1

150

1.0
0.5

1.0
0.5
50

100

2.0

50

100

150

200

50

100

150

200

100
150
Time (frame index)

200

2.0
1.0
0.5

2.0
Y3

1.0
0.5
50

100
150
Time (frame index)

1.0
0.5

200

50

Fig. 2. Example spectra of (a) speech sources, (b) observations, (c) masks and (d) separated signals (N ¼ 3, M ¼ 2).

Step 2. Feature extraction: If the sources sk ðf ; tÞare sufﬁciently sparse, separation can be realized bygathering the time–frequency points ðf ; tÞ where
only one signal sk is estimated to be dominant. To
estimate such time–frequency points, some featuresHðf ; tÞ are calculated by using the frequency–domain observation signals xðf ; tÞ. Here, Hðf ; tÞ is avector that consists of certain geometric features.
Generally, the level ratios and phase differences
between observations are utilized for Hðf ; tÞ. Previously employed features are discussed in Section
3, and our newly proposed feature is introduced in
Section 4.
Step 3. Clustering: Then the features Hðf ; tÞare grouped into N clusters C 1 ; . . . ; C N , where N is
the number of possible sources. Formerly, suchclustering was undertaken manually [3,18],
with a kernel density estimation [20] or with an
ML-based gradient method [21]. On the other
hand, if we can employ a standard clustering
algorithm such as the k-means algorithm or
hierarchical clustering [22], the clustering procedure
will be automated and simpliﬁed. In this work weutilize the k-means clustering algorithm [22] with agiven source number N. The clustering criterion is
to minimize the total sum J of the Euclidean
distances (ED) between cluster members and theircentroids ck :

J¼

M
X
k¼1

Jk ;

Jk ¼

X
Hðf ;tÞ2C kkHðf ; tÞ  ck k2 .

(7)

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

After setting appropriate initial centroids ck
ðk ¼ 1; . . . ; NÞ, this J can be minimized by thefollowing iterative updates:
C k ¼ fHðf ; tÞ j k ¼ argmin kHðf ; tÞ  ck k2 g,

ð8Þ

k

ck

E½Hðf ; tÞH2C k ,ð9Þwhere E½H2C k is a mean operator for the membersof a cluster C k . The cluster members are determinedby (8). If the feature Hðf ; tÞ is properly chosen, theneach cluster corresponds to an individual source.
Here, it should be noted that the k-meansclustering utilizes the ED kHðf ; tÞ  ck k2 , not theMahalanobis distance (MD) ðHðf ; tÞ  ck ÞT R1kðHðf ; tÞ  ck Þ, where Rk is the covariance matrix ofcluster k. That is, k-means assumes clusters of amultivariate isotropic variance Rk ¼ I for all k,
where I denotes an identity matrix.
Step 4. Mask design: Next, the separated signals
yk ðf ; tÞ are estimated based on the clustering result.
We design a time–frequency domain binary maskthat extracts the time–frequency points of eachcluster(1; Hðf ; tÞ 2 C k ;M k ðf ; tÞ ¼
ð10Þ
0 otherwise.
Example of binary mask spectra is shown in Fig. 2(c).
Then, applying the binary masks (Fig. 2(c)) to one ofthe observations (Fig. 2(b)) xJ ðf ; tÞ, we obtainseparated signals (Fig. 2(d)):
yk ðf ; tÞ ¼ M k ðf ; tÞxJ ðf ; tÞ,
where J is a selected sensor index.
Step 5. Separated signal reconstruction: At the endof the ﬂow (Fig. 1), we obtain outputs yk ðtÞ byemploying an inverse STFT (ISTFT) and theoverlap-and-add method [25],yk ðtÞ ¼

S 1
1Xymþl ðtÞ;A l¼0 kð11Þ

where A ¼ 12T=S is a constant for the Hanningwindow case,8Pyðf ; mÞe|2pfr ;>< f 2f0;T1 f s ;...;T1
T f sgymðmSptpmS þ T  1Þ;
k ðtÞ ¼>:0 (otherwise);and r ¼ t  mS.

1837

3. Discussion of features
Because the binary mask approach depends
strongly on the clustering of the feature vectors
Hðf ; tÞ, the selection of an appropriate feature
vector Hðf ; tÞ is essential to this approach. In thissection, we provide examples of the features Hðf ; tÞ
including the previously utilized feature. We also
test how each feature will be clustered by the kmeans algorithm.
3.1. Features
Most previous methods utilized the level ratio
and/or phase difference between observations as
their features Hðf ; tÞ. The previously proposedfeatures can be summarized as
jx2 ðf ; tÞj
x2 ðf ; tÞ T
; arg
(12)
Hðf ; tÞ ¼
jx1 ðf ; tÞj
x1 ðf ; tÞand some examples are shown in Table 1.
Such features (12) represent geometric information on sources and sensors, if the sources are
sufﬁciently sparse. Let us assume that the mixing
process is expressed byhjk ðf Þ  ljk exp½|2pf tjk ,(13)where ljk X0 and tjk are the attenuation and thetime delay from source k to sensor j. If there is noreverberation (i.e., an anechoic situation), ljk and tjk
are determined solely by the geometric distributionof the sources and sensors. If the sources are sparse(5), the feature vector (12) becomesTl2kHðf ; tÞ ¼
; 2pf ðt2k  t1k Þ ; 9 k.
(14)l1k
We can see that Hðf ; tÞ contains geometric information on the dominant source sk at each timefrequency point ðf ; tÞ.
To avoid frequency dependence in the phasedifference (14), some authors have employed afrequency normalization that involves dividing thephase difference by 2pf or 2pfc1 d where c is the
propagation velocity and d is the sensor spacing (seeTable 1). The former is utilized in [3,18] and thelatter gives the directions of arrival (DOA) of
sources if the sensor spacing d is given correctly[26]. If we do not use such frequency normalization,
we have to solve the permutation problem amongfrequencies after clustering the features [16,17].
Moreover, frequency normalization makes it

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

1838

Table 1
Typical features and their separation performance (SIR improvement in dB) with the k-means algorithm
Feature Hðf ; tÞ
(A)

h

(B)

h

ðAÞ0

h

(D)
(E)
(F)
(G)

[18]

h
iiT
jx2 ðf ;tÞjx2 ðf ;tÞ1
jx1 ðf ;tÞj ; 2pfc1 d arg x1 ðf ;tÞ
"
#
h
i T
jx2 ðf ;tÞjx2 ðf ;tÞ
1
1

;
arg
jx2 ðf ;tÞj 2pf
jx ðf ;tÞjx ðf ;tÞ
1

(C)

iiT

jx2 ðf ;tÞj 1x2 ðf ;tÞjx1 ðf ;tÞj ; 2pf arg x1 ðf ;tÞ

hjx1 ðf ;tÞjx2 ðf ;tÞ1
2pf arg x1 ðf ;tÞ

Opt.(ED)

Opt. (MD)

1.9

8.3

14.0

5.7

14.1

14.0

1.8

7.9

14.0

10.5

14.0

14.0

11.6

14.0

14.0

5.2

7.9

14.3

12.4

14.1

14.2

12.2

14.1

14.1

1

i

h
i
1
arg xx2 ðfðf ;tÞ
2pfc1 d
1 ;tÞ
h
h
iiT
jx1 ðf ;tÞj jx2 ðf ;tÞj 1x2 ðf ;tÞAðf ;tÞ ; Aðf ;tÞ ; 2pf arg x1 ðf ;tÞ
hhiiT
jx1 ðf ;tÞj jx2 ðf ;tÞjx2 ðf ;tÞ1
Aðf ;tÞ ; Aðf ;tÞ ; 2pfc1 d arg x1 ðf ;tÞ
hi
arg½xj ðf ;tÞ=xJ ðf ;tÞȲj ðf ; tÞ ¼ jxj ðf ; tÞj exp |aj fHðf ; tÞ[3]

k-means[26]

H̄ðf ; tÞ=kH̄ðf ; tÞk‘‘opt.’’ shows the performance with the known centroid and two distance measures: the Euclidean distance (ED) (8) and the Mahalanobis
distance (MD). N ¼ 3, M ¼ 2. The performance difference between features C and D was caused by the convergence criteria for theqﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
PM2
k-means. Aðf ; tÞ ¼
j¼1 jxj ðf ; tÞj . aj : A weight parameter introduced in Section 4.1.

possible to apply the method to short data without
signiﬁcant performance degradation [27].
3.2. Clustering result with k-means algorithm
Previously, features Hðf ; tÞ were clustered manually [3,18], or with an ML-based gradient method[21]. In contrast, in this subsection, we attempt toemploy the well-known k-means clustering algorithm [22], which can both automate and simplify
the clustering. We show that the previously utilizedfeature (A) cannot be clustered well by the k-means
algorithm and provide possible reasons for this.
Table 1 shows the separation performance (thesignal to interference ratio (SIR) improvement: see
Appendix A) when we cluster each feature with the
k-means algorithm. In Table 1, we also show the
optimal results with known centroid values, which
are calculated with known sources and impulse
responses (unblind). Here, we utilized two omnidirectional microphones with a 4 cm spacing for threespeech sources set at 30 , 70 and 135 , where thedistance between the microphones and sources was
50 cm and the room reverberation time was 128 ms.
We investigated eight combinations of speakers andaveraged the separation results. From Table 1, it
can be seen that all features perform similarly when
the centroids are known and MD-based clustering isused. Therefore, the separation problem amounts toﬁnding a feature that leads to accurate centroid

estimates blindly. However, we can also see thatsome features (A), ðAÞ0 , (B) and (E) cannot achieve
good separation performance with the k-means.
There are two main reasons for this. The ﬁrst
is related to the outliers of the level ratio jx2 j=jx1 j.
Fig. 3 shows examples. We can see several large
values in the level ratio of (A), although we used
omnidirectional microphones where jx2 j=jx1 j  1.
Due to the outliers in the level ratio, the phase of(A) cannot be clustered (Fig. 3 (A)), although thephase terms themselves can be clustered (Table 1(C)). This is the reason for the poor performancewith (A) and (B). Such outliers occur at too manytime-frequency points for them to be removed
without degrading the separated sound quality.
We found that when we normalize the level ratios
as seen in features (E) and (F), they become p1 andprevent such outliers (Fig. 3 (F)). Therefore,features (E) and (F) provide better performance
than (A) and (B). However, the performance with(E) is still insufﬁcient.
This suggests a second reason: namely that thephase term of feature (A) is too small. This is more
important and more fatal than the ﬁrst reason. For
multivariate clustering with the k-means algorithm,
the level ratios and phase differences should have
similar variances. This is because the k-means
assumes distributions of isotropic variance. However, the phase term of feature (A) is far smaller(Fig. 3 (A)) than the level ratio. The poor

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

1839

histogram (log scale)

105
104
103
102

histogram
8000

6000

4000

0

2000

101
2

phase difference (x 10-3)

1.5
1
0.5
0
-0.5
-1
-1.5
-2
20

40

60
level ratio

80

100

120

1500
1000

histogram
2000

1500

0

0
3

500

500

1000

histogram

2000

phase difference

2
1
0
-1
-2
-3
0

0.2

0.4
0.6
level ratio

0.8

1

Fig. 3. Example histograms with features (A) and (F). For each feature, top: histogram of the level ratio term of each feature, bottom left:
the contour plot of the two-dimensional histogram, bottom right: histogram of the phase difference term of each feature. In the contour
plots,  denotes the cluster centroids ck obtained by the k-means algorithm. In (F), we plot only two components½jx1 ðf ; tÞj=Aðf ; tÞ; 1=2pfc1 d arg½x2 ðf ; tÞ=x1 ðf ; tÞT .

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

1840

Table 2
Performance (SIR improvement in dB and SDR in dB: see
Appendix A for the deﬁnitions) and computation time with the kmeans algorithm and the GMM ﬁtting
k-means

GMM

Feature SIR
imprv.

SDR

time (s) SIR
imprv.

SDR

time (s)

(A)
ðAÞ0

1.9
1.8

7.5
7.5

5.0
7.1

11.8
11.3

9.4
9.8

26.3
26.3

(E)
(F)

5.2
12.4

6.0
10.3

4.6
2.5

12.8
14.2

8.8
9.7

33.7
26.5

N ¼ 3, M ¼ 2. GMM ﬁtting needed 6 Gaussians.

performance with (A) and (E) result from thisimbalance between the level ratio and phasedifference terms. With features (B) and (F), where
the phase is divided by 2pfc1 d, the phase difference
becomes larger (see Fig. 3 (F)). Therefore, feature(F), where both the level ratios and phase difference
are normalized appropriately, achieves good performance (see Table 1) with the k-means algorithm.
It should be noted that if we can handle thedifferent variance with, for example the Gaussian
mixture model (GMM) ﬁtting, features (A) and (E)can also work as shown in Table 2. In our
experiments, the GMM ﬁtting with only N ¼ 3
Gaussians did not work. We needed 6 Gaussians forthe successful ﬁtting, and therefore, utilized theposteriors of 3 dominant Gaussians as separation
masks. This selection of the number of Gaussians
required a lot of trial and error. Furthermore, it
should be noted that we had to set appropriate
initial values for the mean and variance of theGaussians carefully for the GMM ﬁtting. If theselection is successful, the GMM ﬁtting works. If
not, it does not work at all. Table 2 shows that theappropriate GMM can achieve reasonable performance even for features (A) and (E). From Table 2,
we can say that the reason for the poor performancewith the k-means for features (A) and (E) arises
from the different variances of the level ratio andthe phase difference.
We would like to note that the GMM ﬁtting
needs sufﬁcient computational time. Table 2 also
gives the calculation time with the k-means and
GMM ﬁtting. Here, we separated mixtures of 5 s
with AthlonXP 3200þ and MATLAB 6.5. The
calculation time was measured with the cputime
command of MATLAB. The computation time
provided in Table 2 is the averaged result for eight

speaker combinations. We can see from Table 2 thatthe k-means with feature (F) achieves sufﬁciently
high performance with a shorter computationaltime than realized with the GMM ﬁtting. The small
computational cost of the k-means is attractive.
In conclusion, we found that feature (F) provides
more accurate clustering result than other features
when either the GMM ﬁtting or the k-meansclustering are employed. Moreover, clustering
is successfully and effectively executed with the
k-means, when we use normalized level ratios andphase differences such as feature (F) (see Fig. 3 (F)and Table 1 (F)).
4. Proposed new feature
Based on the clustering results described in theprevious section, we propose a new feature that can
be clustered by the k-means algorithm. We also
extend the feature to a multiple sensor version,
where we can utilize more than three sensors
arranged non-linearly to separate two- or threedimensionally located sources. As our method canbe considered as an extension of the DUET, we callour method Multiple sENsor dUET: MENUET.
Our proposed feature has the same property as
feature (F) in Section 3, that is, the level ratios andphase differences are appropriately normalized. Our
new feature also has a parameter that controls theweight for the level ratios and phase differences.
Moreover, our normalization does not require
sensor position information. This allows us to apply
our method to an arbitrarily arranged sensor array.
Because the basic scheme is the same as that inFig. 1, here we focus mainly on our new feature
vector.
4.1. New feature
Our new feature employs the normalized level
ratios and phase differences between multiple
observations:Hðf ; tÞ ¼ ½HL ðf ; tÞ; HP ðf ; tÞT ,

(15)

where
HL ðf ; tÞ ¼


jx1 ðf ; tÞjjxM ðf ; tÞj;...;,Aðf ; tÞAðf ; tÞ

ð16ÞHP ðf ; tÞ





1
x1 ðf ; tÞ1xM ðf ; tÞ
¼
arg
arg
;...;
ð17Þ
a1 f
xJ ðf ; tÞaM fxJ ðf ; tÞ

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

ﬃ
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
PM2
In the above equations, Aðf ; tÞ ¼
jx
ðf
;
tÞj
,J
j¼1 j
is the index of one of the sensors, and aj(j ¼ 1; . . . ; MÞ is a positive weighting constant. Bychanging aj , we can control the weights for the level
ratio and phase difference information of the observed
signals; a larger value puts weight on the level ratio
and a smaller value emphasizes the phase difference.
The normalized level ratio has the property of
0pYLj ðf ; tÞp1, where YLj is the jth component of
HL . This can prevent the outliers discussed in theprevious section.
An appropriate value for the phase weight isaj ¼ a ¼ 4pc1 d max , where c is the propagation
velocity and d max is the maximum distance1 between
sensor J and sensor 8 j 2 f1; . . . ; Mg. Let us provide
the reason for this. Here, we use the mixing model(13) and, without loss of generality, we assume thatthe delay parameter tjk is determined by the path
difference l jk  l Jk :tjk ¼ ðl jk  l Jk Þ=c,(18)

where l jk is the distance from source k to sensor j.
We also use the fact that the maximum distanced max between the sensors is greater than the
maximum path difference:
max jl jk  l Jk jpd max .
j;kUsing these assumptions and the mixing model (13),YPj ðf ; tÞ, which is the jth component of HP ðf ; tÞ,
becomes
xj ðf ; tÞ
2pc1 ðl jk  l Jk Þ 2pc1 d max
1
arg
p
.
¼
aj f
xJ ðf ; tÞ
aj
aj
(19)
Because the level ratio is normalized to have therange 0pHL ðf ; tÞp1, the phase difference YPj ðf ; tÞshould also be normalized so that it has a similar
range. If we allow YPj ðf ; tÞ to have the range
12pHP ðf ; tÞp12 (note that YPj ðf ; tÞ can take anegative value), we have equality in (19) whenaj ¼ a ¼ 4pc1 d max . That is, a ¼ 4pc1 d max realizes
the same range width as that of the level ratio.
4.2. Modified proposed feature
We can modify our proposed new feature (15) byusing the complex representation,
Yj ðf ; tÞ ¼ YLj ðf ; tÞ exp½|YPj ðf ; tÞ,1

(20)

If we do not have an accurate value for d max , we may use arough positive constant, as shown in Section 5.2.4.

1841

where YLj and YPj are the jth components of (16)and (17). This modiﬁcation can also be realizedby [28]
arg½xj ðf ; tÞ=xJ ðf ; tÞȲj ðf ; tÞ ¼ jxj ðf ; tÞj exp |, ð21Þ
aj f
Hðf ; tÞH̄ðf ; tÞ=kH̄ðf ; tÞk,

ð22Þ

where H̄ðf ; tÞ ¼ ½Ȳ1 ðf ; tÞ; . . . ; ȲM ðf ; tÞT . Feature(22) is our modiﬁed feature, where the phasedifference information is held in the argument term(21), and the level ratio is normalized by the vector
norm normalization (22). The weight parameter aj
has the same property as (15); however, a ¼
4c1 d max should be the lower limit for successful
clustering (see Appendix B).
Now the normalized vectors Hðf ; tÞ (22) are
M-dimensional complex vectors, and therefore theclustering of the features will be carried out in an Mdimensional complex space. The unit-norm normalization (22) makes the distance calculation in theclustering (7) easier, because it projects the vector on a hyper unit sphere. If the features Hðf ; tÞand the cluster centroid ck are on the unit sphere, i.e., kHðf ; tÞk ¼ kck k ¼ 1, the square distancekHðf ; tÞ  ck k2 ¼ 2ð1  ReðcH
k Hðf ; tÞÞÞ. That is, theminimization of the distance kHðf ; tÞ  ck k2 isequivalent to the maximization of the real part ofthe inner product cH
k Hðf ; tÞ, whose calculation is lessdemanding in terms of computational complexity.
5. Experiments
5.1. Experimental conditions
We performed experiments with measured impulse responses hjk ðlÞ in a room as shown in Figs. 4
and 5. The room reverberation times RT60 were 128
and 300 ms. We used the same room for bothreverberation times but changed the wall condition.
We also changed the distance R between the sensors
and sources. The distance variations were R ¼ 50,
110, and 170 cm (see Figs. 4 and 5). Mixtures weremade by convolving the measured impulse responses in the room and 5-s English speeches. For
the anechoic test, we simulated the mixture by using
the anechoic model ((13) and (18)) and the mixture
model (1). The sampling rate was 8 kHz. The STFT
frame size T was 512 and the window shift was T=4.
Unless otherwise noted, we utilized modiﬁed
feature (22) with aj ¼ a ¼ 4c1 d max for the features, because the computational cost of distanceARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

1842

for M ¼ 2), and the elevation fk ¼ 0 for all sources
k. Note that these initial values of d j and qk were
not exactly the same as those in each setup.
The separation performance was evaluated interms of the SIR improvement and the signal todistortion ratio (SDR). Their deﬁnitions are foundin Appendix A. We investigated eight combinationsof speakers and averaged the separation results.

445 cmMicrophone (120cm height, omni-directional)
Loudspeaker (120cm height)• Reverberation time: 128 ms or 300ms• Room height: 250cm
S2
120°

210°

4cm

256cm
180cm

355cm

S3

4cm

5.2. Separation results

4cm
S4
315°

R cm

S1
45°

R = 50, 110, 170 (cm)Fig. 4. Room setup ðN ¼ 4; M ¼ 3Þ.
445 cmMicrophone (120cm height, omni-directional)
Loudspeaker (height is written with H)• Reverberation time : 128 ms or 300ms• Room height: 250cm

S3(180°, H120cm)(240°, H149cm)
256cm

S2
(120°, H98cm)
4.0cm
5.0cm

180cm

355cm

S4

5.1cm

S1
R cm
(315°, H170cm)

S5

5.0cm

(45°, H120cm)
5.5cm

R = 50, 110, 170 (cm)Fig. 5. Room setup ðN ¼ 5; M ¼ 4Þ.

calculation is lowered (see Section 4.2). We utilized
the k-means algorithm for the clustering, where thenumber of sources N was given. We set the initialcentroids of the k-means using the far-ﬁeld model
where the frequency response hjk ðf Þ is given as
hjk ðf Þ  exp½|2pfc1 dTj qk , and using the same
normalization as each feature. Here, c is the
propagation velocity of the signals, and the threedimensional vectors dj and qk represent the locationof sensor j and the direction of source k, respectively[29]. The sensor locations dj (j ¼ 1; . . . ; MÞ were on
almost the same scale in each setup, and the initial
directions qk were set so that they were as scatteredas possible. Concretely, we utilized the sensor vector
qk ¼ ½cos yk cos fk ; sin yk cos fk ; sin fk T where theazimuth of kth source was set as yk ¼ P  k(k ¼ 1; . . . ; N, P ¼ 2p=N for MX3 and P ¼ p=N

5.2.1. With two sensors
First, we tested our feature with two sensors
under the condition described in Section 3, andcompared the result with that of previous features.
Table 1 in Section 3.2 shows the result. The
proposed feature (15) corresponds to (F) and themodiﬁed feature (22) is shown as (G). We obtained
better separation performance with our proposedfeatures than with other features (A)–(E). A
comparison of the performance achieved with our
proposed method and with the GMM ﬁtting isshown in Table 2. The comparison was investigatedonly for the two sensor case. Our proposed feature(F) achieves high performance with the k-means
within a shorter computation time than with the
GMM ﬁtting. Moreover, we can see that our
proposed feature (F) is also suitable for the GMM
ﬁtting. A comparison with the MAP approach canbe found in [16]. In [16], it is shown that the
proposed method yields better performance in termsof SIR than the MAP approach. It is also pointed
out that proposed method causes larger non-linear
distortion in its outputs than the MAP approach.
We also compared our proposed method with the
DESPRIT algorithm [24], using a linear array of
three microphones for four sources. It should be
noted that the previous DESPRIT limits its arrayshape to a linear array or two sets of congruent
arrays, as discussed in Section 1. In the experiments,
we did not see big difference in performance
between our MENUET and the DESPRIT. Thatis, the proposed MENUET also works with a linear
array (i.e., one-dimensional array).
Note that two sensors/linear arrays do not workwhen the sources are placed at axisymmetricallocations with respect to the microphone array,because they have the equivalent features in (12).
5.2.2. With two-dimensional three sensors
Here, we show the separation results obtainedwith three sensors arranged two-dimensionally

ARTICLE IN PRESS

5.2.4. Weight parameter aHere, we examine the relationship between thephase weight parameter a and the separation
performance. As mentioned in Section 4.1, whena is large the level ratio is emphasized, and when a issmall the phase difference is emphasized. Fig. 7
shows the relationship when N ¼ 4 and M ¼ 3(Fig. 4) with the proposed feature (15) and themodiﬁed feature (22). Note that a ¼ 2p corresponds
to the previous feature (A) (Table 1).
Fig. 7(a) shows the result with the proposedfeature (15). It achieved high performance when awas sufﬁciently small. This is because the phasedifference between the sensors was more reliable
than the level ratio, due to our microphone setup.
As a became small, the performance saturated.
On the other hand, the performance degraded asa became large. This is caused by the imbalance
between level ratio and phase difference terms,
because the phase term becomes too small whena becomes large.
With modiﬁed feature (22), we obtained the best
performance around a ¼ 4c1 d max (Fig. 7(b)). Thisis because a ¼ 4c1 d max realizes the full use of thephase difference information (Appendix B), which ispreferable for our sensor setting. As a became large,
the performance degraded. When ao4c1 d max theperformance also worsened. It should be remembered that, with the modiﬁed feature, a ¼ 4c1 d max
should be the lower limit (see Section 4.2). When

1843

18
16
14
12
10
8
6
4
2
0

anechoic
128ms, 50cm
128ms, 110cm
128ms, 170cm
300ms, 50cm
300ms, 110cm
300ms, 170cm

SIR

[dB]

5.2.3. With four sensors
We also applied our method to a three-dimensional sensor array arranged non-uniformly (Fig. 5).
Here, the system knew only the maximum distanced max (5.5 cm) between the reference microphone andthe others. To avoid the spatial aliasing problem, we
utilized frequency bins up to 3100 Hz in this setup.
Fig. 6(c) shows the separation result when N ¼ 5
and M ¼ 4. Fig. 6(c) shows that our proposed newfeature can be applied to such three-dimensional
microphone array systems.

SDR

18
16
14
12
10
8
6
4
2
0

anechoic
128ms,50cm
128ms,110cm
128ms,170cm
300ms,50cm
300ms,110cm
300ms,170cm

SIR

[dB]

(Fig. 4). Note that sources were also distributed
two-dimensionally.
Fig. 6(a) shows the separation result when N ¼ 4
and M ¼ 3. We can see that our proposed method
achieved good separation performance with the
non-linear sensor arrangement. We also evaluatedthe performance for N ¼ 5; M ¼ 3, where thesource positions were 45 , 120 , 210 , 280 and345 , and obtained good performance (Fig. 6(b)).

[dB]

S. Araki et al. / Signal Processing 87 (2007) 1833–1847

SDR

18
16
14
12
10
8
6
4
2
0

anechoic
128ms,50cm
128ms,110cm
128ms,170cm
300ms,50cm
300ms,110cm
300ms,170cm

SIR

SDR

Fig. 6. Average SIR improvement and SDR for each condition.
The error bar shows the standard deviation for all outputs andcombinations: (a) N ¼ 4; M ¼ 3 (average input SIR  4:8[dB]); (b) N ¼ 5; M ¼ 3 (average input SIR  6:0 [dB]); (c) N ¼
5; M ¼ 4 (average input SIR  6:3 [dB]).

ao4c1 d max , the distance measure (7) for theclustering is not evaluated correctly (see AppendixB), and therefore, the clustering stage failed and theperformance worsened.
We can also see from Fig. 7 that both proposedfeatures (15) and (22) achieved good performanceover a wide a range. This means that we do not need
the exact maximum sensor spacing d max . This allows
us to utilize an arbitrarily arranged sensor array,
although similar distances between pairs of sensors

ARTICLE IN PRESS1844

S. Araki et al. / Signal Processing 87 (2007) 1833–1847
Table 3
Average clarity index [dB]

RT60 ¼ 128 ms
RT60 ¼ 300 ms

R ¼ 50 cm

R ¼ 110 cm

R ¼ 170 cm

45.5
40.5

40.8
34.9

36.6
32.7

are preferable so that the k-means can use all sensor
information optimally.

sparseness increases the approximated W-disjoint
orthogonality rk increases, and vice versa. As seen in
Fig. 8, the sparseness decreases with increases inboth the reverberation time and distance R. That is,
the sparseness decreases when the contribution of
the direct sound is small (see Table 3). In addition,
we can see that an increase in the number of sources
also reduces the sparseness.
It is also important to mention non-linear
distortion in separated signals. There is non-linear
distortion (musical noise) in the outputs with our
method, just as there is in the outputs with theprevious binary mask approaches. The results ofsubjective tests with 10 listeners can be found in [28].
Some sound examples can be found at [31].

5.3. Discussion

6. Conclusion

Fig. 6 also shows the performance tendency inreverberant conditions. The performance degradesas the reverberation time RT60 becomes long.
Moreover, performance degradation was observed
as the distance R became large. This is because,
under long reverberation and/or large distance R
conditions, the direct sound contribution to theimpulse responses becomes smaller, and the source
sparseness (5) and anechoic assumptions (13)cannot hold.
We assessed the way in which reverberation timeand distance R affect the source sparseness (5) andanechoic assumptions (13). Table 3 shows theaverage clarity index C ([30] and Appendix C),
which explains the ratio between direct sound andreverberant sound. Small (large) C means thereverberant sound (direct sound) is large. We cansee from Table 3 that the clarity index C becomes
small as the reverberation time and distance R
increase. That is, when the reverberation time islong and distance R is large, the anechoic assumption (13) seems to become corrupted. For thesparseness measure, we employed the approximate
W-disjoint orthogonality rk ([3] and Appendix C).
Fig. 8 shows the approximate W-disjoint orthogonality under some reverberant conditions. As theWe proposed a novel sparse source separation
method (MENUET) based on the normalizationand clustering of the level ratios and phasedifferences between multiple observations. Our
proposed features can effectively employ the level
ratios and phase differences, and are clustered easily
by the well-known k-means algorithm. It should be
noted that the k-means is optimal when the clusters
are Gaussian; however, this is not always true evenfor our proposed feature (F) (see Fig. 3 (F)).
However, as shown in this paper, the proposedfeature with the k-means achieved sufﬁciently
high performance. Moreover, our feature makes it
easy to employ multiple sensors arranged in a nonlinear/non-uniform way. We obtained promisingexperimental results in a room with weak reverberation even under underdetermined conditions.
Although we provided results solely for underdetermined cases in this paper, our proposed
method can also be applied to (over-) determinedcases [28].
We also reported the separation performance
under some reverberant conditions, where thesparseness and anechoic assumptions were deteriorating. From the results, we saw that the direct andreverberant ratio is important for the current sparse

Fig. 7. Relationship between a and separation performance when
N ¼ 4, M ¼ 3, R ¼ 50 cm, and RT60 ¼ 128 ms. (a) with feature(15) and (b) with modiﬁed feature (22). In (a), ‘‘P’’ denotes theperformance with H ¼ HP , and ‘‘L’’ means with H ¼ HL .

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847

100
900ms(N=4)
128ms, 50cm(N=4)

80r_k [%]

300ms, 170cm(N=4)

700ms(N=5)
128ms, 50cm(N=5)

60300ms, 170cm(N=5)

50
40
30
20
−10

0

10
z [dB]

20

smaller than that of the level term. On the other
hand, when a in (21) is too small, the phase changes
too fast and causes a kind of aliasing problem.
Moreover, it is important that the distance measure(7) of the clustering holds the condition: jY  Y0 j
increases monotonically as j arg½Y  arg½Y0 j increases. However, if the phase term is larger thanp=2, such a monotonic increase cannot hold. That isthe phase term should have the following relationship:pp
 p arg½Yp .
2
230

Fig. 8. Source sparseness rk ðzÞ for some reverberant conditions.

source separation. The sparse source separation inreverberant conditions is still an open problem.

(B.1)

Let us model the mixing process as (13) and,without loss of generality, we assume that the delay
parameter tjk is determined by the path differencel jk  l Jk :tjk ¼ ðl jk  l Jk Þ=c,

Appendix A. Performance measures

1845

(B.2)

where l jk is the distance from source k to sensor j.
This assumption makes tJk ¼ 0. Substituting themixing model (B.2) and (13), and the sparsenessassumption (5) into (21) and (22) yields

The SIR improvement was calculated by
OutputSIRi  InputSIRi ,whereP
jxJi ðtÞj2
InputSIRi ¼ 10 log10 P Pt
ðdBÞ,
ðA:1Þ
2
tj
kai xJk ðtÞjPjy ðtÞj2
OutputSIRi ¼ 10 log10 P Pt ii
ðdBÞ, ðA:2Þ
2
tj
kai yik ðtÞjP
where xJk ðtÞ ¼ l hJk ðlÞsk ðt  lÞ and yik ðtÞ is thecomponent
P of sk that appears at output yi ðtÞ:yi ðtÞ ¼ N
k¼1 yik ðtÞ.
The SDR is employed to evaluate the soundquality:P2
t jxJi ðtÞjSDRi ¼ 10 log10 PðdBÞ,2
t jxJi ðtÞ  byii ðt  DÞj(A.3)
where b and D are parameters used to compensate
for the amplitude and phase difference between xJi
and yii .
Appendix B. Value of a for modiﬁed feature (21)In this section, we show the required condition
for the phase weight parameter a for modiﬁed
feature (21). Because the modiﬁed feature (22) is acomplex vector, we have to consider the phase termwhen we perform clustering. When a in (21) is toolarge, the variance of the phase term becomesljk
2pc1 ðl jk  l Jk Þ
exp |
,
(B.3)
aj
Dk
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
PM 2 ﬃ
where Dk ¼
j¼1 ljk .
From the condition (B.1) and Eq. (B.3), the lower
limit of a is given asHðf ; tÞ 

 1
2pc ðl jk  l Jk Þj arg½Yj ¼ 

aj
 1

2pc d max  p
p ,
p
 2
aj
aj X4c1 d max .

ðB:4Þ
ðB:5Þ

In (B.4), we used the fact that maxj;k jl jk  l Jk j
pd max .
From (B.5), we can conclude that the phaseparameter a ¼ 4c1 d max should be the minimum
value to maintain the relationship (B.1). In addition(B.1) has equality when a ¼ 4c1 d max , which meansthat the phase difference information is mostscattered. That is, the weight with a ¼ 4c1 d max
allows us to make full use of the phase difference
information. This is a preferable property for small
sensor array systems (e.g., see Section 5), wherephase differences between sensors are more reliable
than level ratios for clustering.

ARTICLE IN PRESS1846

S. Araki et al. / Signal Processing 87 (2007) 1833–1847

Appendix C. Measures for reverberation andsparseness assessments
The clarity index [30]
R 80 ms 2h ðtÞ dt
C ¼ 10 log10 R01
ðdBÞ2
80 ms h ðtÞ dt
explains the ratio between direct and reverberant
sound. Small (large) C means the reverberant sound(direct sound) is large.
The sparseness measure, that is the approximate
W-disjoint orthogonality, is deﬁned as [3]P2
ðf ;tÞ kFðk;zÞ ðf ; tÞsk ðf ; tÞkrk ðzÞ ¼
 100ð%Þ,(C.1)
P2
ðf ;tÞ ksk ðf ; tÞk
where Fðk;zÞ is a time-frequency binary mask thathas a parameter z

1 20 logðjsk ðf ; tÞj=jy^ k ðf ; tÞjÞ4z;
Fðk;zÞ ðf ; tÞ ¼
0 otherwise(C.2)
PN
and y^ k ðf ; tÞ ¼ STFT½ i¼1;iak si ðtÞ (sum of interference components). The approximate W-disjoint
orthogonality rk ðzÞ indicates the percentage of theenergy of source k for time–frequency points where
it dominates the other sources by z dB. A larger
approximate W-disjoint orthogonality rk ðzÞ means
more sparseness, and vice versa.
References[1] S. Haykin (Ed.), Unsupervised Adaptive Filtering (Volume
I: Blind Source Separation), Wiley, New York, 2000.
[2] A. Hyvärinen, J. Karhunen, E. Oja, Independent Component Analysis, Wiley, New York, 2001.
[3] Ö. Yılmaz, S. Rickard, Blind separation of speech mixtures
via time–frequency masking, IEEE Transactions on SP 52(7) (2004) 1830–1847.
[4] H. Buchner, R. Aichner, W. Kellermann, Blind source
separation for convolutive mixtures: a uniﬁed treatment, in:Y. Huang, J. Benesty (Eds.), Audio Signal Processing for
Next-Generation Multimedia Communication Systems,
Kluwer Academic Publishers, Dardrecht, 2004, pp. 255–293.
[5] H. Sawada, R. Mukai, S. Araki, S. Makino, Frequencydomain blind source separation, in: J. Benesty, S. Makino,
J. Chen (Eds.), Speech Enhancement, Springer, Berlin, 2005,
pp. 299–327.
[6] S. Amari, S. Douglas, A. Cichocki, H. Yang, Multichannel
blind deconvolution and equalization using the natural
gradient, in: Proceedings of IEEE Workshop on Signal
Processing Advances in Wireless Communications, 1997,
pp. 101–104.
[7] P. Smaragdis, Blind separation of convolved mixtures in thefrequency domain, Neurocomputing 22 (1998) 21–34.

[8] L. Parra, C. Spence, Convolutive blind separation of nonstationary sources, IEEE Trans. Speech Audio Process. 8 (3)(2000) 320–327.
[9] J. Anemüller, B. Kollmeier, Amplitude modulation decorrelation for convolutive blind source separation, in: Proceedings of the ICA 2000, 2000, pp. 215–220.
[10] S. Araki, R. Mukai, S. Makino, T. Nishikawa, H.
Saruwatari, The fundamental limitation of frequency
domain blind source separation for convolutive mixturesof speech, IEEE Trans. Speech Audio Process. 11 (2) (2003)
109–116.
[11] F. Theis, E. Lang, C. Puntonet, A geometric algorithm forovercomplete linear ICA, Neurocomputing 56 (2004)
381–398.
[12] P. Boﬁll, M. Zibulevsky, Blind separation of more sources
than mixtures using sparsity of their short-time Fouriertransform, in: Proceedings of the ICA2000, 2000, pp. 87–92.
[13] L. Vielva, D. Erdogmus, C. Pantaleon, I. Santamaria,
J. Pereda, J.C. Principe, Underdetermined blind source
separation in a time-varying environment, in: Proceedings ofthe ICASSP 2002, 2002, pp. 3049–3052.
[14] P. Boﬁll, Underdetermined blind separation of delayedsound sources in the frequency domain, Neurocomputing 55(2003) 627–641.
[15] A. Blin, S. Araki, S. Makino, Underdetermined blindseparation of convolutive mixtures of speech using time–
frequency mask and mixing matrix estimation, IEICE Trans.
Fundam. E88-A (7) (2005) 1693–1700.
[16] S. Winter, W. Kellermann, H. Sawada, S. Makino, MAPbased underdetermined blind source separation of convolutive mixtures by hierarchical clustering and l1-norm
minimization, EURASIP J. Adv. Signal Process. (2007),
Article ID 24717.
[17] J.M. Peterson, S. Kadambe, A probabilistic approach forblind source separation of underdetermined convolutive
mixtures, in: Proceedings of the ICASSP 2003, vol. VI, 2003,
pp. 581–584.
[18] A. Jourjine, S. Rickard, Ö. Yılmaz, Blind separation of
disjoint orthogonal signals: demixing N sources from 2
mixtures, in: Proceedings of the ICASSP 2000, vol. 12, 2000,
pp. 2985–2988.
[19] M. Aoki, M. Okamoto, S. Aoki, H. Matsui, T. Sakurai, Y.
Kaneda, Sound source segregation based on estimating
incident angle of each frequency component of input signals
acquired by multiple microphones, Acoust. Sci. Technol. 22(2) (2001) 149–157.
[20] N. Roman, D. Wang, G.J. Brown, Speech segregation basedon sound localization, J. Acousit. Soc. Am. 114 (4) (2003)
2236–2252.
[21] S. Rickard, R. Balan, J. Rosca, Real-time time–frequency
based blind source separation, in: Proceedings of the ICA
2001, 2001, pp. 651–656.
[22] R.O. Duda, P.E. Hart, D.G. Stork, Pattern Classiﬁcation,
second ed., Wiley Interscience, New York, 2000.
[23] R. Balan, J. Rosca, S. Rickard, Non-square blind source
separation under coherent noise by beamforming andtime–frequency masking, in: Proceedings of the ICA 2003,
2003, pp. 313–318.
[24] T. Melia, S. Rickard, C. Fearon, Histogram-based blindsource separation of more sources than sensors using a
DUET-ESPRIT technique, in: Proceedings of the EUSIPCO
2005, 2005.

ARTICLE IN PRESS
S. Araki et al. / Signal Processing 87 (2007) 1833–1847[25] S. Araki, S. Makino, H. Sawada, R. Mukai, Reducing musical noise by a ﬁne-shift overlap-add method
applied to source separation using a time–frequency mask,
in: Proceedings of the ICASSP 2005, vol. III, 2005,
pp. 81–84.
[26] S. Araki, S. Makino, A. Blin, R. Mukai, H. Sawada,
Underdetermined blind separation for speech in realenvironments with sparseness and ICA, in: Proceedings ofthe ICASSP 2004, vol. III, 2004, pp. 881–884.
[27] S. Araki, H. Sawada, R. Mukai, S. Makino, Normalized observation vector clustering approach for sparsesource separation, in: Proceedings of the EUSIPCO 2006,
2006.

1847

[28] S. Araki, H. Sawada, R. Mukai, S. Makino, A novel blindsource separation method with observation vector clustering, in: Proceedings of the 2005 International Workshop on
Acoustic Echo and Noise Control (IWAENC 2005), 2005,
pp. 117–120.
[29] S. Araki, H. Sawada, R. Mukai, S. Makino, DOA
estimation for multiple sparse sources with normalized
observation vector clustering, in: Proceedings of the ICASSP
2006, vol. 5, 2006, pp. 33–36.
[30] ISO 3382: Acoustics-measurement of the reverberation time ofrooms with reference to other acoustical parameters (1997).
[31] hhttp://www.kecl.ntt.co.jp/icl/signal/araki/xcluster_ﬁne.
htmli.