➠➡
FREQUENCY DOMAIN BLIND SOURCE SEPARATION USING
SMALL AND LARGE SPACING SENSOR PAIRS
Ryo Mukai Hiroshi Sawada Shoko Araki Shoji Makino
NTT Communication Science Laboratories, NTT Corporation
2–4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619–0237, Japan{ryo,sawada,shoko,maki}@cslab.kecl.ntt.co.jp
ABSTRACTThis paper presents a method for solving the permutation
problem of frequency domain blind source separation (BSS)
when the number of source signals is large, and the potential
source locations are omnidirectional. We propose a combination of small and large spacing sensor pairs with various
axis directions in order to obtain proper geometrical information for solving the permutation problem. Experimental
results show that the proposed method can separate a mixture of six speech signals that come from various directions,even when two of them come from the same direction.
1. INTRODUCTION
Blind source separation (BSS) is a technique for estimating original source signals using only observed mixtures.
When the source signals are s i (t)(i = 1, ..., N ), the signals observed by sensor j are x j (t)(j = 1, ..., M ), andthe separated signals are y k (t)(k = 1, ..., N ), the BSS
Nmodel can be described as: x j (t) =i=1 (hji ∗ si )(t),Myk (t) = j=1 (wkj ∗ xj )(t), where hji is the impulse response from source i to sensor j, w kj are the separating
filters, and ∗ denotes the convolution operator.
Figure 1 shows a flow of the BSS in frequency domain.
A convolutive mixture in the time domain is converted into
multiple instantaneous mixtures in the frequency domain.
Therefore, we can apply an ordinary independent component analysis (ICA) algorithm [1] in the frequency domain
to solve a BSS problem in a reverberant environment. Using a short-time discrete Fourier transform, the model isapproximated as: X(ω, n) = H(ω)S(ω, n), where, ω isthe angular frequency, and n represents the frame index.
The separating process can be formulated in each frequency
bin as: Y(ω, n) = W(ω)X(ω, n), where S(ω, n) =[S1 (ω, n), ..., SN (ω, n)]T is the source signal in frequency
bin ω, X(ω, n) = [X1 (ω, n), ..., XM (ω, n)]T denotes theobserved signals, Y(ω, n) = [Y 1 (ω, n), ..., YN (ω, n)]T isthe estimated source signal, and W(ω) represents the separating matrix. W(ω) is determined so that Y i (ω, n) andYj (ω, n) become mutually independent.
The ICA solution suffers permutation and scaling ambiguities. This is due to the fact that if W(ω) is a solution, then D(ω)P(ω)W(ω) is also a solution, where D(ω)
is a diagonal complex valued scaling matrix, and P(ω) is
an arbitrary permutation matrix. We thus have to solve thepermutation and scaling problems to reconstruct separated
signals in the time domain.

0-7803-8251-X/04/$17.00 ©2004 IEEE

V-1

Sourcesignals

s1

DFT

ICA ω

Permutation
problem

Scalingproblem IDFT

time

s2

time

freq

freq

W (ω)

P (ω)

D(ω)

Permutation
misalignment
Multiple instantaneous mixtures

Convolutive mixtures

Fig. 1. Flow of frequency domain BSS
There is a simple and reasonable solution for the scaling problem: D(ω) = diag{[P(ω)W(ω)] −1 }, which is obtained by the minimal distortion principle (MDP) [2], andwe can use it. On the other hand, the permutation problem
is complicated, especially when the number of source signals is large. Time domain BSS does not suffer the permutation problem, however it takes much computational time
as compared with frequency domain BSS [3]. Therefore we
adopt the frequency domain approach.
Many methods have been proposed for solving the permutation problem, and the use of geometric information,such as beam patterns [4, 5, 6], direction of arrival (DOA)and source locations[7], is effective approach. We have proposed a robust method by combining the correlation based
method [8] and the DOA based method [4, 5], which almost
completely solves the problem for 2-source cases [9]. However it is insufficient when the number of signals is large
or when the signals come from a similar direction. In thispaper, we propose a method for obtaining proper geometric information for solving the permutation problem in suchcases.
2. GEOMETRIC INFORMATION FOR SOLVINGPERMUTATION PROBLEM
2.1. Invariant in ICA solution
If a separating matrix W is calculated successfully
and it extracts source signals with scaling ambiguity,
D(ω)W(ω)H(ω) = I holds (except for singular frequencybins). Because of the scaling ambiguity, we cannot obtain
H simply from the ICA solution. However, the ratio of elements in the same column H ji /Hj  i is invariable in relation
to D, and given byHji[W−1 D−1 ]ji[W−1 ]ji
=
=
,
(1)
−1
−1
Hj  i[W D ]j  i[W−1 ]j  i

ISCAS 2004

➡

➡
where [·]ji denotes ji-th element of the matrix. We canestimate several types of geometric information related tosource signals by using this invariant. The estimated information is utilized for solving the permutation problem.

ri,jj  = 1.4
ri,jj  = 0.71
ri,jj  = 1.6
r  = 0.63
ri,jj  = 2.0 r  = 0.5i,jj
i,jj

12.2. DOA estimation with ICA solution
0.5

We can estimate the DOA of source signals by using theabove invariant [10]. With a farfield model, a frequency
response is formulated as:

z(m)
0

2

Interpretation of the ICA solution by a nearfield modelyields other geometric information. When we adopt thenearfield model, including the attenuation of the wave,
Hji (ω) is formulated as:−1
1
ejωc (qi −pj )(5)Hji (ω) =qi − pj 
where qi represents the location of source i. By taking theratio of (5) for a pair of sensors j and j  we obtain:qi − pj   jωc−1 (qi −pj −qi −pj )Hji /Hj  i =e. (6)qi − pj 
By using the modulus of (6) and (1), we have:qi − pj    [W−1 ]ji .
(7)=qi − pj 
[W−1 ]j  i 

By solving (6) for q i , we have a sphere whose center O i,jj
and radius Ri,jj are given by:1(pj  − pj ),Oi,jj = pj − 2(8)ri,jj − 1
ri,jj
(pj  − pj ),Ri,jj =  2(9)ri,jj − 1
where ri,jj = |[W−1 ]ji /[W−1 ]j  i |. Thus, we can estimate
a sphere (Ôi,jj , R̂i,jj ) on which qi exists by using the result of ICA W and the locations of the sensors p j and pj  .
Figure 2 shows an example of the spheres determined by (7)for various ratios r i,jj . This procedure is valid for sensor
pairs with a large spacing.
The models (2) and (5) are simple approximation without the multi-path propagation and reverberation, however

V-2

1.5

1

0.5

0

−0.5

−1

y(m)



[W−1]ji 
[W−1]j i 

−1.5

−2

0

x(m)

−0.5
−1

Fig. 2. Example of spheres determined by eq.(7) (p j =[0, 0.3, 0], pj  = [0, −0.3, 0])
Nearly verticalto sensor pair axis

Estimated DOA (deg)

180

140

120

100

80

Sources

60

40

20

0

0

100

200

300

400

Sensors

500

600

700

Nearly horizontal
to sensor pair axis

180

160

800

900

160

140

100

Sensors

80

60

40

20

0

1000

Sources

120

0

100

200

Frequency bin

300

400

500

600

700

800

900

1000

Frequency bin

Fig. 3. Source locations and estimated DOAs
6|∆θ̂/∆ arg(Ĥ)|

2.3. Estimation of sphere with ICA solution



−1

−1

= ejωc pj −pj  cos θi,jj ,(3)where θi,jj is the direction of source i relative to the sensor
pair j and j  . By using the argument of (3) and (1), we canestimate:arg(Hji /Hj  i )θ̂i,jj = cos−1 −1ωc (pj − pj  )arg([W−1 ]ji /[W−1 ]j  i ).
(4)= cos−1
ωc−1 (pj − pj  )This procedure is valid for sensor pairs with a small spacing.






∆ 
ri,jj  =qi = [x, y, z]

Estimated DOA (deg)= ejωc−1 aTi (pj −pj  )

0.5

−0.5

−1 T

Hji (ω) = ejωc ai pj ,(2)where c is the speed of wave propagation, a i is a unit vector
that points to the direction of source i, and p j represents alocation of sensor j. According to this model, we have
Hji /Hj  i

1

pj pj 5f=500 Hz43

f=1000 Hz

2

f=2000 Hz

1

f=4000 Hz
0

20

40

60

80

100

θ̂(deg.)

120

140

160

180

Fig. 4. Sensitivity of DOA estimation
we can obtain information for classifying signals by usingthem.
3. SENSITIVITY AND AMBIGUITY IN SOURCELOCATION ESTIMATION
3.1. Sensitivity of DOA estimation
DOA estimation is sensitive to source locations. Figure 3
shows examples of DOA estimation with two different
source locations. When the source signals are almost in
front of a sensor pair, their directions can be estimated robustly. However, when the signals are nearly horizontal tothe axis of the pair, the estimated directions tend to havelarge errors. This can be explained as follows.
When we denote an error in calculated arg(H ji /Hj  i )as ∆ arg(Ĥ), and an error in θ̂i,jj as ∆θ̂, the ratio|∆θ̂/∆ arg(Ĥ)| can be approximated by the partial derivative of (4):
 

 ∆θ̂
 

1

 

(10)

≈
.
 ∆ arg(Ĥ)   ωc−1 pj − pj   sin(θ̂i,jj ) 

➡

➡
S

DOA of S1 and S2 estimated bysmall spacing sensor pairs

S2

Unreliable
cone

S1

Reliablecones

Estimated
source locations(Ô2, R̂2)
(Ô1, R̂1)

Fig. 5. Combination of small spacing sensor pairs with different axes
Figure 4 shows examples of this value for several frequency
bins. We can see that ∆ arg( Ĥ) causes a large error in theestimated DOA when the direction is near the axis of thesensor pair. Therefore, we consider the estimated DOA to
be unreliable in such cases.

Spheres estimated bylarge spacing sensor pairFig. 6. Combination of small spacing sensor pairs and alarge spacing sensor pair3.2. Ambiguity of estimated DOA

estimate DOAs by using the rows of W −1 (ω) corresponding to the small spacing microphone pairs (1-3, 2-4, 1-2 and
2-3). Figure 8 shows a histogram of estimated DOAs. We
can find five clusters in this histogram, and one cluster is
twice the size of the others. This implies that these are six
source signals, and two of them come from the same direction (about 150 ◦ ). We can solve permutation for four
sources by using this information (Fig. 9).
Then, we apply the estimation of spheres to the signals that belong to the large cluster by using the rows of
W−1 (ω) corresponding to the large spacing microphonepairs (7-5, 7-8, 6-5 and 6-8). Figure 10 shows estimatedradiuses for S4 and S5 regarding the microphone pair 7-5.
Finally, we can classify the signals into six clusters.
Unfortunately, the classification by the estimated location tends to be inconsistent especially in a reverberant environment. In many frequency bins, several signals are assigned to the same cluster, and such classification is inconsistent. We solve the permutation only for frequency bins
with a consistent classification, and we employ a correlation based method [9] for the rest. The correlation based
method solves the permutation so that the inter-frequency
correlation for neighboring or harmonic frequency bins becomes maximized. In addition, we use the spectral smoothing method proposed in [11] to make separating filters in thetime domain from the result of ICA W(ω).
The performance is measured from the signal-toinference ratio (SIR). The portion of y k (t) that comes fromMsi (t) is calculated by yki (t) = j=1 (wkj ∗ hji ∗ si )(t). If
we solve the permutation problem so that s i (t) is output toyi (t), the SIR for yk (t) is defined as: 
SIRk = 10 log[ t ykk (t)2 / t ( i=k yki (t))2 ] (dB).

Another problem regarding DOA estimation is ambiguity.
When we use only one pair of sensors or a linear array,we estimate a cone rather than a direction. If we assume
a plane on which sources exist, the cone is reduced to two
half-lines. However, the ambiguity of two directions that aresymmetrical with respect to the axis of the sensor pair still
remains.This is a fatal problem when the source locations
are omnidirectional.
When the spacing between sensors is larger than half awavelength, spatial aliasing causes another ambiguity, but
we do not consider this here.
3.3. Resolving sensitivity and ambiguity
A nonlinear arrangement of sensors is suitable for resolving
both sensitivity and ambiguity. We propose a combinationof small and large spacing sensor pairs that have various
axis directions.
By using the DOA estimation described in Sec.2.2 withthe small spacing sensor pairs that have different axis directions, we can estimate cones which have various vertex
angles for one source direction. Because of the sensitivity
explained in Sec.3.1, we assume that obtuse cones are reliable, and acute cones are unreliable. Then, we can determine a bearing line pointing to a source direction by using
the reliable cones (Fig. 5).
Even when some signals come from the same or a similar direction, we can distinguish between them by using theinformation provided by the large spacing sensor pair described in Sec.2.3. The source locations can be estimatedby combining the estimated direction and spheres (Fig. 6).
Then, we can classify separated signals in the frequency domain according to the estimated source locations.
4. EXPERIMENTS

We measured SIRs for three permutation solving strategies: only the correlation based method (“C”), the estimatedDOAs and correlation (“D+C”), and the combination of
estimated DOAs, spheres and correlation (“D+S+C”, proposed method). We also measured input SIRs by using the
mixture observed by microphone 1 for the reference (“Input
SIR”). The results are summarized in Table 2.
Our proposed method succeeded in separating six

We carried out experiments for 6 sources and 8 microphones using speech signals convolved with impulse responses measured in a room. The room layout is shownin Fig. 7. Other conditions are summarized in Table 1. The
experimental procedure is as follows.
First, we apply ICA to xj (t)(j = 1, ..., 8), and calculate
separating matrix W(ω) for each frequency bin. Then we

V-3

➡

➠

s2

-30 deg

2 cm

s1

mic.3

355 cmmic.5

s3

90 deg

mic.6

30 cm
60 cm mic.7
120 cms4
180 cm

mic.8

S6
100

200

300

400

500

600

700

800

900

1000

5

Estimated radius (m)Room height: 250 cm

8 kHz
6s
2048 point (256 ms)
512 point (64 ms)
Infomax (complex valued)

4
3.5
3

S5

2.5
2
1.5
1

S4

0.5
0

0

100

200

300

400

500

600

700

800

900

1000

Frequency bin

Fig. 10. Estimated radiuses for S 4 and S5

500

Table 2. Experimental results (dB)

500

Number of Signals

−100

4.5

Table 1. Experimental conditions

450

400

400
350

300

300

Input SIR

250

200

C
D+C
D+S+C

200
150

100

100
50
0

S1

−50

Fig. 9. Permutation solved by using DOAs

Fig. 7. Room layout

0

S2

0

Frequency bin

Microphones (omnidirectional, height: 135 cm)
Loudspeakers (height: 135 cm)

Sampling rate
Data length
Frame length
Frame shift
ICA algorithm

S3

50

0

s6
s5

100

−150

-150 deg
150 deg

S4 S5

150

mic.1
mic.2 4 cm mic.4

Estimated DOA (deg)445 cmReverberation time: 130 ms
225 cm30 deg

−150

-150

−100

-100

−50

-50

0

0

50

50

Estimated DOA (deg)

100

100

150

150

Fig. 8. Histogram of estimated DOAs obtained by using
small spacing microphone pairs
speech signals. We can see that the discrimination obtainedby using estimated spheres is effective in improving the separation performance for signals coming from the same direction.
5. CONCLUSION
We proposed the combination of small and large spacing
microphone pairs with various axis directions in order to obtain proper geometrical information for solving the permutation problem in frequency domain BSS. In experiments,
our method succeeded in the separation of six speech signals, even when two come from the same direction. The
computation time was about 1 min. for 6 s. data.
6. REFERENCES[1] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, John Wiley & Sons, 2001.
[2] K. Matsuoka and S. Nakashima, “Minimal distortion principle for blind source separation,” in Proc. of Intl. Workshop
on Independent Component Analysis and Blind Signal Separation (ICA’01), 2001, pp. 722–727.
[3] K. Matsuoka, Y. Ohba, Y. Toyota, and S. Nakashima, “Blindseparation for convolutive mixture of many voices,” in Proc.
IWAENC 2003, 2003, pp. 279–282.
[4] S. Kurita, H. Saruwatari, S. Kajita, K. Takeda, and F. Itakura,“Evaluation of blind signal separation method using direc-

V-4

SIR1

SIR2

SIR3

SIR4

SIR5

SIR6

-8.3
4.4
4.5
12.3

-6.8
2.6
10.8
5.6

-7.8
4.0
14.4
14.5

-7.7
9.2
4.5
7.6

-6.7
3.6
5.4
8.9

-5.2
-2.0
8.8
10.8

ave.
-7.1
3.7
8.1
10.0

tivity pattern under reverberant conditions,” in Proc. of
ICASSP’00, 2000, pp. 3140–3143.
[5] M. Z. Ikram and D. R. Morgan, “A beamforming approach topermutation alignment for multichannel frequency-domain
blind speech separation,” in Proc. of ICASSP’02, 2002,
vol. 1, pp. 881–884.
[6] L. C. Parra and C. V. Alvino, “Geometric source separation: Merging convolutive source separation with geometricbeamforming,” IEEE Trans. on Speech and Audio Processing, vol. 10, no. 6, pp. 352–362, Sep. 2002.
[7] V. C. Soon, L. Tong, Y. F. Huang, and R. Liu, “A robust
method for wideband signal separation,” in Proc. of ISCAS’93, 1993, pp. 703–706.
[8] F. Asano, S. Ikeda, M. Ogawa, H. Asoh, and N. Kitawaki,“A combined approach of array processing and independent
component analysis for blind separation of acoustic signals,”in Proc. ICASSP 2001, May 2001, pp. 2729–2732.
[9] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust
and precise method for solving the permutation problem offrequency-domain blind source separation,” in Proc. Intl.
Symp. on Independent Component Analysis and Blind Signal Separation (ICA2003), 2003, pp. 505–510.
[10] H. Sawada, R. Mukai, and S. Makino, “Direction of arrival estimation for multiple source signals using independent component analysis,” in Proc. of ISSPA’03, 2003, vol. 2,
pp. 411–414.
[11] H. Sawada, R. Mukai, S. de la Kethulle de Ryhove, S. Araki,
and S. Makino, “Spectral smoothing for frequency-domain
blind source separation,” in Proc. IWAENC 2003, 2003, pp.
311–314.