SOLVING THE PERMUTATION PROBLEM OF FREQUENCY-DOMAIN BSS
WHEN SPATIAL ALIASING OCCURS WITH WIDE SENSOR SPACING
Hiroshi Sawada

Shoko Araki

Ryo Mukai

Shoji Makino

NTT Communication Science Laboratories, NTT Corporation
2-4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japan

ABSTRACTThis paper describes a method for solving the permutation problem
of frequency-domain blind source separation (BSS). The method analyzes the mixing system information estimated with independent
component analysis (ICA). When we use widely spaced sensors or
increase the sampling rate, spatial aliasing may occur for high frequencies due to the possibility of multiple cycles in the sensor spacing. In such cases, the estimated information would imply multiple possibilities for a source location. This causes some difficultywhen analyzing the information. We propose a new method designed to overcome this difficulty. This method first estimates themodel parameters for the mixing system at low frequencies where
spatial aliasing does not occur, and then refines the estimations byusing data at all frequencies. This refinement leads to precise parameter estimation and therefore precise permutation alignment. Experimental results show the effectiveness of the new method.
1. INTRODUCTION
The technique for estimating individual source components fromtheir mixtures at multiple sensors is known as blind source separation (BSS) [1]. With acoustical applications of BSS, such as solving
a cocktail party problem, signals are generally mixed in a convolutive manner with reverberations. Let s1 , . . . , sN be source signals
and x1 , . . . , xM be sensor observations. The convolutive mixture
model is formulated asxj (t) =

N


k=1

hjk (l) sk (t − l), j = 1, . . . , M ,(1)lwhere t represents time and hjk (l) represents the impulse response
from source k to sensor j. In a practical room situation, impulse responses hjk (l) can have thousands of taps even with an 8 kHz sampling rate. This makes the convolutive BSS problem very difficultcompared with the BSS of simple instantaneous mixtures.
An efficient and practical approach for such convolutive mixtures is frequency-domain BSS [2–9], where we apply a short-time
Fourier transform (STFT) to the sensor observations xj (t). In thefrequency domain, the convolutive model (1) can be approximatedas an instantaneous mixture model at each frequency:
xj (f, t) =N

hjk (f )sk (f, t), j = 1, . . . , M,

(2)

k=1

where f represents frequency, t is now down-sampled with the distance of the frame shift, hjk (f ) is the frequency response from source
k to sensor j, and sk (f, t) is a frequency-domain representation of asource signal sk (t).
Independent component analysis (ICA) [10] is a major statistical
tool for BSS. In the frequency-domain approach, ICA is employed in

1­4244­0469­X/06/$20.00 ©2006 IEEE

STFT

ICA

ISTFT

inverse

Permutation& Scaling

Fig. 1. Flow of frequency-domain BSS

each frequency bin with the instantaneous mixture model (2). This
makes the convergence of ICA relatively fast compared with convolutive ICA where the mixture model (1) is explicitly assumed. However, the ICA solution has permutation ambiguity. Even if we change
the order of the separated signals at the ICA output, it is still an ICA
solution. This causes the permutation problem of frequency-domain
BSS. We need to align the permutation of each frequency bin so thatthe frequency components of the same source are grouped together.
Various methods have been proposed for solving the permutation problem: 1) making the separation matrices smooth [2, 3], 2)
maximizing the correlation of separated signal envelopes [4], 3) analyzing the directivity patterns calculated from the separation matrices [5, 6] and 4) analyzing the mixing system information estimatedwith ICA [7, 8]. The third and fourth methods utilize similar kinds ofinformation, such as estimated directions of sources. However, thefourth method is more general than the third one, as it can easily beapplied to a situation where there are more than two sources [7, 8].
We have experimentally shown that the fourth (or the third) method
provides a robust solution for the permutation problem [9]. When
we employ such a method, we prefer the sensor spacings to be nolarger than half the minimum wavelength of interest. We typicallyuse a 4 cm sensor spacing for an 8 kHz sampling rate to satisfy thiscondition. If sensor spacing is wider, spatial aliasing might occur athigh frequencies [11], and the ICA solution in such a frequency bin
implies multiple possibilities for a source location.
In this paper, we propose a new method for dealing with a situation where the distance between sensors is larger than half the wavelength. Although a wider sensor spacing is discussed in [6], it isbased on the third method, which is hard to generalize for more than
two sources. With the new proposed method working effectively, we
can use widely spaced sensors to achieve better separation for low
frequencies or we can increase a sampling rate, for example up to
16 kHz, to obtain better speech recognition accuracy for separated
signals.
2. FREQUENCY-DOMAIN BSS
This section presents an overview of frequency-domain BSS. Figure 1 shows the flow that we consider in this paper. First, the sen-V ­ 77ICASSP 2006

Table 1. Experimental setup A

sourceSource directions
Sensor spacing
Source distance from sensors
Sampling rate
Reverberation time
Frame size of STFT
Source signal
Propagation velocity

sensorFig. 2. Direct-path model

sor observations (1) are converted into frequency-domain time-series
signals (2) by an STFT. Let us rewrite (2) in a vector notation:
x(f, t) =N

hk (f )sk (f, t),

(3)

k=1

where x = [x1 , . . . , xM ]T is the vector of observed signals andhk = [h1k , . . . , hM k ]T is the vector of frequency responses fromsource sk to all sensors.
Then, complex-valued instantaneous ICA [10] is applied to the
mixtures x(f, t) to obtain separated frequency components:y(f, t) = W(f ) x(f, t),(4)where y = [y1 , . . . , yN ]T is the vector of separated frequency components and W = [w1 , . . . , wN ]H is an N × M separation matrix.
With the method considered here, the inverse of separation matrix W (or the Moore-Penrose pseudoinverse if W is not square) is
calculated in each frequency bin. The inverse is represented as[a1 , · · · , aN ] = W−1 , ai = [a1i , . . . , aM i ]T ,(5)which we call basis vectors obtained by ICA, because the mixture
x(f, t) is represented with a linear combination of basis vectors bymultiplying W−1 and (4):x(f, t) =Nai (f )yi (f, t).

(6)

i=1

By comparing (6) and (3), we see that a basis vector ai (f ) represents the same information as hk (f ) up to permutation and scaling
ambiguity if ICA works well. The use of different subscripts, i andk, indicates the permutation ambiguity.
Next, the permutation ambiguity is aligned. This paper focuseson the fourth method discussed in the introduction, which analyzes
the mixing system information represented by the basis vectors ai (f ).
Section 3 will discuss the way to analyze the information and thenalign the permutation.
Then, the scaling ambiguity of ICA is aligned byyi (f, t) ← aJ i (f )yi (f, t),(7)where J is the index of a reference sensor (see [8] for the rationaleof this operation). Finally, time-domain output signals yi (t) are obtained from separated frequency components yi (f, t) by an inverse
STFT (ISTFT).

3. PERMUTATION ALIGNMENT

70◦ and 150◦ (2 sources)
dmax = 20 cm (2 sensors)
120 cm
16 kHzRT60 = 130 ms
128 msspeech of 3 s
c = 340 m/s

3.1. Assumption and basic idea
We assume a simple direct-path model (Fig. 2) for the mixing system, even though in reality signals are mixed in a multi-path model(1). This simplified model is expressed in the frequency domain:hjk (f ) = λjk · e−j 2πf τjk ,(8)

where τjk and λjk ≥ 0 are the time delay and attenuation fromsource k to sensor j, respectively. Since we cannot distinguish thephase (or amplitude) of sk (f, t) and hjk (f ), these two parameters
can be considered to be relative. Thus, without loss of generality, we
normalize them by
τjk = (djk − dJ k )/c,

M

(9)

λ2 = 1,j=1 jk(10)
where djk is the distance from source k to sensor j (Fig. 2), and cis the propagation velocity of the signal. Normalization (9) makesτJ k = 0, i.e. the relative time delay is zero at the reference sensor J.
As explained in Sec. 2, a basis vector ai (f ) represents the same
information as hk (f ) up to permutation and scaling ambiguity. Following (9) and (10), the scaling ambiguity is aligned by the operationai −j arg(aJi )
(11)
e
ai ←
||ai ||
which makes arg(aJ i ) = 0 and ||ai || = 1. Now, the task for the
permutation problem is to determine a permutation Πf that relates
the subscript i and k with i = Πf (k), and to estimate parametersτjk , λjk that make the model (8) match the aji (f ) element of thebasis vector. This can be formulated so as to find Πf , τjk and λjkthat minimize the cost function:
J =

N 
M


|aji (f ) − λjk · e−j 2πf τjk |2 , i = Πf (k), (12)f ∈F k=1 j=1

where F is the set of frequencies that we have to consider.
To make the discussion here intuitively understandable, we performed an experiment with setup A shown in Table 1. This was asimple M = N = 2 case, but the sensor spacing was 20 cm, which
induced spatial aliasing for a 16 kHz sampling rate. Figure 3 shows
the argument of a21 and a22 after the normalization (11) where we
set J = 1 as a reference sensor. The arguments of a1i (f ) are notshown because they are all zero. The relative time delay τ21 andτ22 can be estimated from these data. However, we see some circular jumps at high frequencies, which are caused by spatial aliasing.
They complicate the estimation of τ21 and τ22 .
3.2. For frequencies without spatial aliasing

In this section, we propose a new method for solving the permutation problem. It is based on an analysis of the mixing system information estimated by ICA and represented in basis vectors (5). The
method can handle a situation where spatial aliasing occurs becauseof a wide sensor spacing or a high sampling rate.

Let us first consider the lower frequency range
FL = {f : −π < 2πf τjk < π, ∀ j, k}(13)

where we can guarantee that spatial aliasing does not occur. Let
dmax be the maximum distance between the reference sensor J andV ­ 78

0.85

Fig. 3. Arguments of a21 and a22 before permutation alignment

Fig. 4. Arguments of ā21 and ā22 after permutations are aligned
only for frequency range FL = {f : 0 < f < 850 Hz}

any other sensor. Then the relative time delay is bounded by
max |τjk | ≤ dmax /c

(14)

jk

and therefore FL can be defined as

c
}.
(15)
2dmax
For frequency range FL , appropriate Πf for (12) can be obtained by the method shown in our previous work [8], which further
normalizes the basis vectors (11) to remove frequency dependence:
FL = {f : 0 < f <arg[aji (f )]āji (f ) ← |aji (f )| exp j
4f c−1 dmax



.

(16)

The rationale of dividing the argument by 4f c−1 dmax is discussed
in [8]. With this operation, the cost function (12) is converted into
J¯ =

M
N
 


|āji (f ) − h̄jk |2 , i = Πf (k)(17)

f ∈FL k=1 j=1

where

π c · τjk]
2 dmaxis a frequency-normalized model. In a vector notation,h̄jk = λjk · exp[−j

J¯ =

N
 

||āi (f ) − h̄k ||2 , i = Πf (k),

(18)

(19)

f ∈FL k=1

h̄k ←



N


||āΠ(k) (f ) − h̄k ||2 ,of spatial aliasing. This means that the frequency normalization (16)does not work for these higher frequencies, and therefore we need a
new method to minimize the cost function (12).
3.3. For frequencies where spatial aliasing may occur

where āi = [ā1i , . . . , āM i ]T and h̄k = [h̄1k , . . . , h̄M k ]T . Here we
see that frequency-normalized basis vectors āi (f ) and the frequencynormalized model h̄k do not depend on frequency. Therefore, J¯ canbe minimized efficiently by a clustering algorithm that iterates thefollowing two updates until convergence:Πf ← argminΠ

Fig. 5. Arguments of a21 and a22 after permutation alignment using
model parameters estimated with data in the low frequency range
FL . Because τ21 and τ22 are not precisely estimated, there are somepermutation errors at high frequencies.

This subsection proposes a new method for deciding the permutation
Πf for frequencies where spatial aliasing may occur. The model
parameters τjk , λjk can be extracted from the frequency-normalized
model h̄jk (18) calculated by (21):
2 dmaxτjk = −
(23)
arg(h̄jk ), λjk = |h̄jk | .
π cThus, a simple way is to use these extracted parameters for the costfunction (12), and decide the permutations Πf by(20)

Πf = argminΠ

k=1āΠf (k) (f ), h̄k ← h̄k /||h̄k ||.

|ajΠ(k) (f ) − λjk · e−j 2πf τjk |2 .

(24)

k=1 j=1

(21)

f ∈FL

The first update (20) optimizes the permutation Πf for each frequency with the current model h̄k . The second update (21) calculates the most probable model h̄k with the current permutations.
This set of updates is very similar to that of the k-means algorithm[12]. After the algorithm has converged, we update the subscript ofthe basis vectors by
ak (f ) ← aΠf (k) (f ), k = 1, . . . , N.

N
M



(22)

Figure 4 shows the arguments of ā21 and ā22 calculated by operation (16) in the setup A experiment. For frequency range FL ,the permutations Πf were aligned and updated by (22). We see twoclusters whose centroids are the two lines represented by arg(h̄21 )and arg(h̄22 ). For frequencies higher than 850 Hz, we see the effectHowever, τjk and λjk estimated only with frequencies in FL may
not be very accurate. Figure 5 shows arg(a21 ) and arg(a22 ) after
permutations were calculated by (24) using model parameters extracted by (23). We see some estimation error for τ21 and τ22 .
A better way is to refine the model parameters τjk and λjk by
the gradient of the cost function J (12):∂J∂J, λjk ← λjk − µ,(25)τjk ← τjk − µ
∂τjk
∂λjk
where µ is a step size parameter, and

∂J
∝ λjk
f · imag[ajΠ(k) (f ) · ej 2πf τjk ] , (26)∂τjkV ­ 79

f ∈F

∂J
∂λjk

∝





λjk − real[ajΠ(k) (f ) · ej 2πf τjk ] , (27)f ∈F

Table 2. Comparison of average signal-to-interference ratio (SIR)
obtained with three different methods for permutation alignment.
Setup A
Setup B

Fig. 6. Argument of a21 and a22 after permutation alignment using
model parameters refined with data at all frequencies. Now τ21 andτ22 are precisely estimated, and permutations are aligned correctly.
445 cmMethod I
11.5 dB
9.0 dB

Method II
15.4 dB
11.2 dB

Method III
17.2 dB
13.8 dB

tations were aligned for all frequency bins with τjk and λjk refinedwith data at all frequencies (corresponding to Fig. 6). We see thatmethod III provides superior results to the other methods. Since separated frequency components generated by ICA were the same for all
methods, the differences among these average SIRs were due solely
to the preciseness of the permutation alignment.
5. CONCLUSION

Reverberation time: RT60 = 130 ms
LoudspeakerWe have proposed a new method for the permutation problem offrequency-domain BSS, which analyzes mixing system information
estimated with ICA. The method can handle a situation where thesensor spacing is larger than half the minimum wavelength and thus
spatial aliasing occurs. Experimental results clearly show the effectiveness of the proposed method.

120cm
Microphones
30cm (omnidirectional)

355 cmInput SIR
0.0 dB−3.1 dB

Loudspeaker
Loudspeaker

6. REFERENCES

Room height: 250 cm
Height of microphones and loudspeakers: 135 cm

Fig. 7. Experimental setup B

are the partial derivatives of J with respect to τjk and λjk . The
operations imag[·] and real[·] extract the imaginary and real parts ofa complex number, respectively. We can iteratively update Πf by(24) and (τjk , λjk ) by (25) to obtain better estimations of the modeland consequently better permutations. Note that the structure that
iterates (24) and (25) has the same structure as (20) and (21). Figure
6 shows arg(a21 ) and arg(a22 ) after τjk and Πf were refined by(24) and (25). We see that τ21 and τ22 were precisely estimated andthe permutations were aligned correctly even for high frequencies.
4. EXPERIMENTS
To see the effectiveness of the new method, we conducted experiments to separate 3-second English speeches blindly. We used two
setups. Setup A (Table 1) was used to illustrate the process of themethod with Figs. 3-6. Setup B (Fig. 7) was used to examine thevalidity of the method for a more complicated case. One attractive
feature of the method shown in Sec. 3 is that the system does nothave to know the sensor geometry, but simply the maximum distance dmax from a reference sensor J to the other sensors. In setup
B, we assigned the center microphone as the reference sensor, andprovided the system with the information dmax = 30 cm.
Table 2 shows the BSS results measured with the average signalto-interference ratio (SIR) for 16 combinations of speeches. The
SIR was calculated as the ratio of the power of a target component
and interference components [9]. In method I, permutations werealigned only for a low frequency range FL (corresponding to Fig. 4).
In method II, permutations were aligned for all frequency bins, butwith model parameters τjk and λjk estimated with data in the low
frequency range FL (corresponding to Fig. 5). In method III, permu-[1] S. Haykin, Ed., Unsupervised Adaptive Filtering (Volume I:
Blind Source Separation), John Wiley & Sons, 2000.
[2] P. Smaragdis, “Blind separation of convolved mixtures in thefrequency domain,” Neurocomputing, vol. 22, pp. 21–34, 1998.
[3] L. Parra and C. Spence, “Convolutive blind separation of nonstationary sources,” IEEE Trans. Speech Audio Processing,vol. 8, no. 3, pp. 320–327, May 2000.
[4] N. Murata, S. Ikeda, and A. Ziehe, “An approach to blindsource separation based on temporal structure of speech signals,” Neurocomputing, vol. 41, no. 1-4, pp. 1–24, Oct. 2001.
[5] H. Saruwatari, S. Kurita, K. Takeda, F. Itakura, T. Nishikawa,
and K. Shikano, “Blind source separation combining independent component analysis and beamforming,,” EURASIP Journal on Applied Signal Processing, vol. 2003, no. 11, pp. 1135–
1146, Nov. 2003.
[6] M. Z. Ikram and D. R. Morgan, “Permutation inconsistency
in blind speech separation: Investigation and solutions,” IEEE
Trans. Speech Audio Processing, vol. 13, no. 1, pp. 1–13, Jan.
2005.
[7] R. Mukai, H. Sawada, S. Araki, and S. Makino, “Frequency
domain blind source separation for many speech signals,” inProc. ICA 2004 (LNCS 3195), pp. 461–469. Springer-Verlag,
Sept. 2004.
[8] H. Sawada, S. Araki, R. Mukai, and S. Makino, “Blind extraction of a dominant source signal from mixtures of manysources,” in Proc. ICASSP 2005, Mar. 2005, vol. III, pp. 61–
64.
[9] 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, Sept.
2004.
[10] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, John Wiley & Sons, 2001.
[11] D. H. Johnson and D. E. Dudgeon, Array Signal Processing:Concepts and Techniques, Prentice-Hall, 1993.
[12] R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification,
Wiley Interscience, 2nd edition, 2000.

V ­ 80