Blind extraction of a dominant source from mixtures of manysources using ICA and time-frequency masking
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
Email: {sawada,shoko,ryo,maki}@cslab.kecl.ntt.co.jp

Abstract— This paper presents a method for enhancing a target sourceof interest and suppressing other interference sources. The target source
is assumed to be close to sensors, to have dominant power at these sensors,
and to have non-Gaussianity. The enhancement is performed blindly, i.e.
without knowing the total number of sources or information about eachsource, such as position and active time. We consider a general case where
the number of sources is larger than the number of sensors. We employ
a two-stage process where independent component analysis (ICA) is first
employed in each frequency bin and time-frequency masking is then used
to improve the performance further. We propose a new sophisticatedmethod for selecting the target source frequency components, and also a
new criterion for specifying time-frequency masks. Experimental results
for simulated cocktail party situations in a room (reverberation time was
130 ms) are presented to show the effectiveness and characteristics ofthe proposed method.

I. I NTRODUCTION
The technique for estimating individual source components fromtheir mixtures at sensors is known as blind source separation (BSS)[1], [2]. With some applications such as brain imaging or wirelesscommunications, it makes sense to extract as many source components as possible, because many sources are equally important.
However, with audio applications such as speech enhancement, thesources do not necessarily have equal significance. We often want to
extract only one source that is close to sensors, has dominant power,
and/or has interesting features.
This paper presents a method for extracting a source signalof interest and suppressing other interference sources blindly. Let
us formulate the task. Suppose that a target source s1 and otherinterference sources s2 , . . . , sN are convolutively mixed and observed
at M sensors
xj (t) =

N 
k=1

l

hjk (l) sk (t − l), j = 1, . . . , M ,(1)where hjk (l) represents the impulse response from source k to sensor
j. The goal is to have an output signal y1 (t) that is close to thecomponent of s1 measured at a selected sensor J:
xJ 1 (t) =

N



l

hJ 1 (l) s1 (t − l).

(2)Note that xj (t) = k=1 xjk (t). The task should be performed onlywith the M observed signals. The number of sources N is unknown
and may be larger than M .
The first problem is how to extract the target source s1 blindly.
Even if N could be larger than M , independent component analysis(ICA) [2] with an N = M assumption produces M components thatmaximize an ICA criterion such as non-Gaussianity. We assume thatthe target source s1 is non-Gaussian, close to sensors, and dominant
in the mixtures. Therefore, we expect that one of the M components
corresponds to s1 whose ICA criterion is high.
We employ ICA in the time-frequency domain. The reason is that itis efficient and also fits time-frequency masking, which is discussed in
the next paragraph. An additional operation that should be performed

0-7803-8834-8/05/$20.00 ©2005 IEEE.

is the selection of the s1 component in every frequency bin. Thisis considered to be the permutation problem of frequency-domain
BSS. It has been reported that the selection of a component withmaximum kurtosis works well when the target is speech and theinterferences are babble sources [3]. However, this does not always
work well for a case where the interferences are also speech. Thus,
we exploit the information of basis vectors (8) produced by ICA. Our
previously reported methods estimate the directions [4], [5] and/or the
distances [5] of the sources from the basis vectors, and then cluster theestimations to solve the permutation problem. However, the systemneeds to know the locations of sensors to estimate such geometric
information about the sources. In Sec. II-C, we propose a new method
for solving the permutation problem. With this approach, we do notneed to know the sensor locations, simply the maximum distance
between a sensor and any other sensor. This relaxation makes it easy
to use a non-uniform arrangement of sensors, and also eliminates theneed for sensor calibration.
The next issue is that some interference still remains in theextracted frequency components when N > M . Post filtering [3],[6] can be used to reduce such residual interference. However, itneeds additional adaptation where the step size should be controlled
based on the short-term power of the target. Another approachis time-frequency masking [7]–[11], which is efficient for sources
with sparseness in the time-frequency domain, such as speech. The
performance of time-frequency masking depends on how well we canspecify the time-frequency slots where the target source is active. A
simple way to specify such slots is to calculate the phase and/or
amplitude difference between the observations of different sensors[7], [8]. Another recently proposed approach involves calculating
the power ratio between an input and outputs of a spatial filter(beamformer [9], [11] or ICA [11]). However, such a power-based
criterion depends on the scaling ambiguity of ICA or beamformer
outputs. We propose a new criterion for specifying masks in Sec. IID. It is based on the cosine distance between a sample vector and thebasis vector corresponding to the target. The closeness is calculated
in a spatially whitened space where the target basis vector is expected
to be almost orthogonal to those of interferences. Therefore, the newcriterion does not suffer from the problem of scaling ambiguity.
The next section describes our proposed method. Section III shows
experimental results, and Section IV concludes this paper.
II. T HE PROPOSED METHOD
A. Frequency domain operations
Figure 1 shows the flow of the method discussed here. First, timedomain signals xj (t) sampled at frequency fs are converted intofrequency-domain time-series signals xj (f, τ ) with an L-point shorttime Fourier transform (STFT):

5882

xj (f, τ ) ←L/2−1r=−L/2xj (τ + r) win(r) e−2πf r ,

(3)

STFT

ICA

Time-frequency
masking

Permutation

Fig. 1.

Flow of proposed method

where f ∈ {0, L1 fs , . . . , L−1fs } is a frequency, win(r) is aLwindow that tapers smoothly to zero at each end, such as a Hanning), and τ is a new index representing time.
window 12 (1 + cos 2πr
L
The remaining operations are performed in the frequency domain.
There are two advantages to this. First, the convolutive mixtures (1)can be approximated as instantaneous mixtures at each frequency:Nxj (f, τ ) ≈k=1

hjk (f )sk (f, τ ),xJ 1 (f, τ ) ≈ hJ 1 (f )s1 (f, τ ),sensorFig. 2.

(5)where J should be the same for all frequency bins f . The second
advantage is that the sparseness of a source signal becomes prominent
in the time-frequency domain if the source is colored and nonstationary such as speech. The possibility of sk (f, τ ) being close
to zero is much higher than that of sk (t).
Through several operations, which will be discussed in the following subsections, we have an output ỹ1 (f, τ ), which should be closeto (5) in each frequency bin. At the end of the flow, we have an
output y1 (t) by an inverse STFT (ISTFT) :1
L·win(r)sourceSince s1 is assumed to be a dominant non-Gaussian source, it isstrongly expected that one of y1 , . . . , yM corresponds to s1 and thus
one of a1 , . . . , aM corresponds to h1 .
C. Permutation
The next operation is to find i for each frequency f such thatai (f ) corresponds to h1 (f ). As shown in [4], integrating the basis
vector ai (f ) and signal envelope |yi (f, τ )| information solves the
permutation problem robustly and precisely, and we also employ thisapproach here. In the rest of this subsection, we discuss a new method
for exploiting the basis vector information.
The new method involves normalizing all basis vectors ai (f ), i =fs such that
1, . . . , M , for all frequency bins f = 0, L1 fs , . . . , L−1
Lthey form clusters, each of which corresponds to an individual source.
The normalization is performed by selecting a reference sensor J andcalculatingỹ1 (f, τ ) e 2πf r .

arg[aji (f )/aJ i (f )]āji (f ) ← |aji (f )| exp 
4f c−1 dmax

1 f , ..., L−1 f }f ∈{0, Lss
L

B. Independent component analysis (ICA)
Let us have a vector notation of the mixing model (4) :x(f, τ ) ≈

N

k=1

hk (f )sk (f, τ ),(6)

where x = [x1 , . . . , xM ] is a sample vector and hk =[h1k , . . . , hM k ]T is the vector of frequency responses from source
sk to all sensors. Independent component analysis (ICA) is used asa first step to identify the vector h1 of a dominant source s1 .
Even though the number of independent components N may belarger than the number of sensors M , we employ ICA by assumingthat N is equal to M :(7)where y = [y1 , . . . , yM ]T is a vector of independent components
and W = [w1 , . . . , wM ]H is an M × M separation matrix. In
the experiments shown in Sec. III, we calculated W by using acomplex-valued version of FastICA [2], and improved it further byusing InfoMax [12] combined with the natural gradient [13] whose
nonlinear function is based on the polar coordinate [14].
Then, we calculate the inverse of W to obtain basis vectors[a1 , · · · , aM ] = W−1 , ai = [a1i , . . . , aM i ]T .

āi (f ) ← āi (f ) / ||āi (f )||

M

i=1

ai (f )yi (f, τ ).

(9)

(11)

for āi (f ) = [ā1i (f ), . . . , āM i (f )] .
The next step is to find clusters C1 , . . . , CM formed by normalized
vectors āi (f ). The centroid ck of a cluster Ck is calculated by
ck ←



ā∈Ck

ā/|Ck |, ck ← ck /||ck ||,where |Ck | is the number of vectors in Ck . The clustering criterion is
to minimize the total sum J of the squared distances between cluster
members and their centroid
J =

M

k=1

Jk , Jk =



ā∈Ck

||ā − ck ||2 .

(12)This minimization can be performed efficiently with the k-meansclustering algorithm [15].
This paragraph explains the reason why normalized basis vectorsāi (f ) form a cluster for a source. Let us approximate the multi-path
mixing model (1) by using a direct-path (nearfield) model (Fig. 2)
hjk (f ) ≈−1

x(f, τ ) =

(10)

T

(8)

By multiplying both sides of (7) by W , the sample vector x(f, τ )is represented by a linear combination of basis vectors a1 , . . . , aM :where c is the propagation velocity and dmax is the maximum distance between the reference sensor J and a sensor ∀ j ∈ {1, . . . , M }.
The rationale of this operation will be explained afterwards. Then,
we apply unit-norm normalizationTy(f, τ ) = W(f ) x(f, τ ),Direct-path (nearfield) model(4)

where hjk (f ) is the frequency response from source k to sensor
j, and sk (f, τ ) is a frequency-domain time-series signal of sk (t)obtained by the same operation as (3). The frequency-domain counterpart of (2) isy1 (τ + r) ←

ISTFT



q(f )
exp  2πf c−1 (djk − dJ k ) ,
djk

(13)

where djk > 0 is the distance between source k and sensor j.
We assume that the phase 2πf c−1 (djk − dJ k ) depends on the
distance normalized with the distance to the reference sensor J. This5883

assumption makes the phase zero at the reference sensor J. We also
assume that the attenuation q(f )/djk depends on both the distanceand a frequency-dependent constant q(f ) > 0. By considering the
permutation and scaling ambiguity of ICA, a basis vector and its
elements are represented asai ≈ αi hk ,aji ≈ αi hjk ,Fig. 3.

Angle θ1 calculated in whitened space(14)

where αi represents the scaling ambiguity, and index k, which may
be different from index i, represents the permutation ambiguity.
Substituting (13) and (14) into (10) and (11) yieldsāji (f ) ≈



1
djk D

exp 



π (djk − dJ k )
, D=
2
dmax


M

i=1

1
,
dik 2

which is independent of frequency, and dependent only on thepositions of the sources and sensors. From the fact that maxj,k |djk −dJ k | ≤ dmax , an inequality−π/2 ≤ arg[āji (f )] ≤ π/2
holds. This property is important for the distance measure (12), since|ā − ā | increases monotonically as | arg(ā) − arg(ā )| increases.
Once we have found M clusters C1 , . . . , CM , we need to identify
a cluster that corresponds to the target source s1 . We decide thata cluster CK with the minimum variance K = argmink Jk /|Ck |corresponds to s1 . The rationale behind this is that the mixing model(13) is more valid for s1 than for the other sources. The directpath components of impulse responses hj1 are distinct since s1 isassumed to be close to the sensors. Finally, the output index i foreach frequency f is selected byI(f ) = argmini ||āi (f ) − cK ||2 .

Fig. 4.

where 0 ≤ M(f, τ ) ≤ 1 is a mask specified for each time-frequency
slot (f, τ ). We specify masks based on the angle θ1 (f, τ ) between
a1 (f ) and x(f, τ ) calculated in the space transformed by a whitening
matrix V(f ) = R−1/2 , R = x(τ ) x(τ )H τ . Let z(f, τ ) =V(f )x(f, τ ) be whitened samples and b1 (f ) = V(f )a1 (f ) be thebasis vector in the whitened space. The angle is calculated byθ1 (f, τ ) = arccos

a1 (f )y1 (f, τ ) = (α1 a1 (f )) (y1 (f, τ )/α1 ),for any non-zero complex scalar α1 . This is easily solved by

1
,
(16)
1 + eg(θ1 −θT )where θT and g are parameters specifying the transition point and its
steepness, respectively. As θT becomes smaller, the residual powerthat appears in ỹ1 decreases but the musical noise in y1 increases.
The effectiveness of the above operation depends on the sparsenessof sources. If we assume that the possibility of sk (f, τ ) being close
to zero is very high, (6) can be approximated asx(f, τ ) ≈ hk (f )sk (f, τ ), k ∈ {1, . . . , N },where J is the index of the sensor specified in (5). The reason is asfollows. The goal in each frequency bin is to make y1 (f, τ ) as close
to xJ 1 (f, τ ) defined in (5) as possible. And we can derive relations
xJ 1 (f, τ ) ≈ hJ 1 (f )s1 (f, τ ) ≈ aJ 1 (f )y1 (f, τ ).

Suppose that the permutation ambiguity of ICA is solved at thisstage. The extraction of s1 by ICA (7) is represented by=

w1H h1 s1 (τ ) +N
k=2

w1H hk sk (τ ).

If N ≤ M , w1 satisfies w1H hk = 0, ∀ k ∈ {2, . . . , N } and makes
the second term zero. However, we assume that the number ofsources N is generally larger than M . In this case, there exists a
H∀
set K ⊆ {2, . . . , N } such that w
1 hk = 0, k ∈ K. Thus, y1 (τ )Hwhcontains an unwanted residualk sk (τ ). The purpose
1k∈K
of time-frequency masking is to obtain
another
output ỹ1 (τ ) that

H
w
h
s
contains less power of the residualkk (τ ) than y1 (τ ).
1k∈KTime-frequency masking is performed byỹ1 (f, τ ) = M(f, τ ) · y1 (f, τ ),z(f, τ ) ≈ V(f )h1 (f )s1 (f, τ ) ≈ V(f )a1 (f )y1 (f, τ )(18)

V(f )hk (f )sk (f, τ ).
k=2

(19)

z(f, τ ) ≈D. Time-frequency masking(17)

where k depends on each time-frequency slot (f, τ ). Let us consider
the whitened-space counterpart of (17), while distinguishing between
cases where s1 is the only active source (18) and other cases (19):from (5), the h1 term in (6) and the a1 term in (9).

w1H x(τ )(15)

for each time-frequency slot (Fig. 3). Then, we calculate a mask byusing a logistic function (Fig. 4)

y1 (f, τ ) ← aJ 1 (f )y1 (f, τ ),=|bH
1 (f ) · z(f, τ )|||b1 (f )|| · ||z(f, τ )||M(θ1 (f, τ )) =This means that basis vector aI(f ) (f ) corresponds to h1 (f ).
After we align the index as a1 (f ) ← aI(f ) (f ) and y1 (f, τ ) ←yI(f ) (f, τ ), we solve the scaling ambiguity in (9) :y1 (τ )

Masking functions with three sets of parameters (θT , g)N

If the number of sources N is equal to or less than the numberof sensors M , vectors Vh1 , . . . , VhN in the whitened space areorthogonal to each other. Even if N > M , the vector b1 =Va1 of a dominant source s1 , which points in almost the same
direction as Vh1 , tends to have large angles with the other vectors
Vh2 , . . . , VhN . Figure 3 shows such a case. Therefore, calculating
the angle (15) provides information about whether or not s1 is the
only active source at a time-frequency slot (f, τ ).
III. E XPERIMENTS
We performed experiments to enhance a dominant speech that wasclose to microphones. We measured impulse responses hjk (l) under
the conditions shown in Fig. 5. The speaker positions simulated a
cocktail party situation. Mixtures at the microphones were made byconvolving the impulse responses and 6-second English speechessampled at 8 kHz. Microphone arrangement was 3-dimensional

5884

TABLE I
AVERAGE SIR IMPROVEMENT FOR EACH POSITION (dB)

445 cm25cm

b120
Target source
candidates

90cm

Target position
InputSIR
Only ICAICA and T-F masking (0.375π, 40)
ICA and T-F masking (0.333π, 20)
ICA and T-F masking (0.25π, 20)

80cm
270cm

355 cm

70cm

a120

80cm
c120 (120 cm height)
c170 (170 cm height)

2.5cm

Mic. 3Mic. 1

4cm

2.5cm
2.3cm

Mic. 4

c170
−0.0
13.0
16.9
18.3
21.0

IV. C ONCLUSION
Interference
speech sources

3.5cm Mic. 2• Microphones arranged
3-dimensionally

80cm

150cm

80cm

25cm
25cm

• Reverberation time: RT60 = 130 ms• Room height: 250 cm• Loudspeakers and microphones at height of 120cm except c170

Fig. 5.

a120 b120 c120
1.3
1.5
1.9
11.7 11.8
9.0
15.4 14.6 12.5
16.8 15.8 14.1
19.5 18.2 16.9

Experimental conditions

and non-uniform. The system knew only the maximum distance(4 cm) between the reference microphone (Mic. 1) and others. For
each setup, we selected one of the four speakers (a120, b120,
c120, c170) as a dominant target source, and the others were kept
silent. The six speakers away from the microphones were used asinterferences for every setup. The frame size L of STFT (3) was
1024 (128 ms). The computational time was around 12 seconds
for 6-second speech mixtures. The program was coded in Matlab
and run on Athlon 64 FX-53. The performance was evaluated interms of the signal-to-interference ratio (SIR) improvement, which is
OutputSIR − InputSIR. These two types of SIRs are defined by
InputSIR = 10 log 10

 J1

|x (t)|2 t
| k=1 xJ k (t)|2 t(dB),

|y (t)|2 t
| k=1 y1k (t)|2 t(dB),

OutputSIR = 10 log10

 11where xJ k (t) is defined in (2), and y1k (t)the component of sk
isNthat appears at output y1 (t), i.e. y1 (t) = k=1 y1k (t).
Experiments were conducted with 16 combinations of 7 speeches
for each target position. Table I shows the average SIR improvements
obtained only with ICA, and by the combination of ICA and timefrequency (T-F) masking. The SIR improvements depend on the target
position. Positions a120 and b120 were fairly good for enhancement.
This is because the interferences came from different directions. If
we consider the speaker arrangement 2-dimensionally, positions c120
and c170 seems to be a hard position as many interferences camefrom similar directions. However, the result for position c170 wasvery good. This is because the height of c170 was different fromthose of interferences, and the 3-dimensionally arranged microphones
enable the system to exploit the height difference.
We used three sets of parameters for function (16) specifying amask for each time-frequency slot. The shapes of these functions
are shown in Fig. 4. Table I shows that a smaller θT resulted ingreater SIR improvements by T-F masking. However, some soundswith a small θT were unnatural. We observed that in many cases
parameter (θT , g) = (0.333π, 20) produced natural sounds with
sufficient interference suppression. Some sound examples can befound on our web site [16].

We have presented a method for extracting a dominant target
source and suppressing other interferences. The process of ICA
and following permutation alignment extracts the target source, andestimates the corresponding basis vector. The new method for permutation alignment makes it easy to use a 3-dimensional non-uniform
arrangement of sensors without exact measurement or calibration.
Time-frequency masking in the second stage reduces the power ofthe residuals caused by ICA when N > M . It exploits the sparsenessof sources. We have proposed a new criterion for specifying masks.
It is based on the angle between the target basis vector and a sample
vector, and gives information about whether or not the target sourceis active. The experiments showed good results for extracting adominant source out from six interferences mixed in a real room.
R EFERENCES[1] S. Haykin, Ed., Unsupervised Adaptive Filtering (Volume I: Blind Source
Separation). John Wiley & Sons, 2000.
[2] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis. John Wiley & Sons, 2001.
[3] S. Y. Low, R. Togneri, and S. Nordholm, “Spatio-temporal processing
for distant speech recognition,” in Proc. ICASSP 2004, vol. I, May 2004,
pp. 1001–1004.
[4] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust and precise
method for solving the permutation problem of frequency-domain blindsource separation,” IEEE Trans. Speech Audio Processing, vol. 12, no. 5,
pp. 530–538, Sept. 2004.
[5] R. Mukai, H. Sawada, S. Araki, and S. Makino, “Frequency domain
blind source separation using small and large spacing sensor pairs,” inProc. ISCAS 2004, vol. V, May 2004, pp. 1–4.
[6] R. Mukai, S. Araki, H. Sawada, and S. Makino, “Removal of residualcrosstalk components in blind source separation using LMS filters,” inProc. NNSP 2002, Sept. 2002, pp. 435–444.
[7] M. Aoki, M. Okamoto, S. Aoki, H. Matsui, T. Sakurai, and Y. Kaneda,“Sound source segregation based on estimating incident angle of eachfrequency component of input signals acquired by multiple microphones,” Acoustical Science and Technology, vol. 22, no. 2, pp. 149–157,
2001.
[8] S. Rickard, R. Balan, and J. Rosca, “Real-time time-frequency based
blind source separation,” in Proc. ICA2001, Dec. 2001, pp. 651–656.
[9] N. Roman and D. Wang, “Binaural sound segregation for multisource
reverberant environments,” in Proc. ICASSP 2004, vol. II, May 2004,
pp. 373–376.
[10] S. Araki, S. Makino, A. Blin, R. Mukai, and H. Sawada, “Underdetermined blind separation for speech in real environments with sparsenessand ICA,” in Proc. ICASSP 2004, vol. III, May 2004, pp. 881–884.
[11] D. Kolossa and R. Orglmeister, “Nonlinear postprocessing for blindspeech separation,” in Proc. ICA 2004 (LNCS 3195), Sept. 2004, pp.
832–839.
[12] A. Bell and T. Sejnowski, “An information-maximization approach toblind separation and blind deconvolution,” Neural Computation, vol. 7,
no. 6, pp. 1129–1159, 1995.
[13] S. Amari, “Natural gradient works efficiently in learning,” NeuralComputation, vol. 10, no. 2, pp. 251–276, 1998.
[14] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Polar coordinate based
nonlinear function for frequency domain blind source separation,” IEICE
Trans. Fundamentals, vol. E86-A, no. 3, pp. 590–596, Mar. 2003.
[15] R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification, 2nd ed.
Wiley Interscience, 2000.
[16] [Online]. http://www.kecl.ntt.co.jp/icl/signal/sawada/demo/dominant/5885