IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 2005

1640

INVITED PAPER

Special Section on Multi-channel Acoustic Signal Processing

Blind Source Separation of Convolutive Mixtures of Speech in
Frequency Domain
Shoji MAKINO†a) , Hiroshi SAWADA† , Ryo MUKAI† , Members, and Shoko ARAKI† , Nonmember

SUMMARY
This paper overviews a total solution for frequencydomain blind source separation (BSS) of convolutive mixtures of audio
signals, especially speech. Frequency-domain BSS performs independent
component analysis (ICA) in each frequency bin, and this is more eﬃcient than time-domain BSS. We describe a sophisticated total solution
for frequency-domain BSS, including permutation, scaling, circularity, andcomplex activation function solutions. Experimental results of 2 × 2, 3 × 3,
4 × 4, 6 × 8, and 2 × 2 (moving sources), (#sources × #microphones) in aroom are promising.
key words: blind source separation, convolutive mixtures, independent
component analysis, frequency-domain BSS, microphone array, adaptive
beamformer

1.

Introduction

Blind source separation (BSS) [1]–[3] is an approach forestimating source signals by using only the information of
mixed signals observed at each input channel. The estimation is performed blindly, i.e., without possessing information on each source, such as its location and active time.
Typical examples of such source signals include mixturesof simultaneous speech signals that have been picked up by
several microphones. Its potential audio signal applications
include speech enhancement for speech recognition, teleconferences, and hearing aids. In such applications, signals
are mixed in a convolutive manner with reverberations. This
makes the BSS problem diﬃcult. We need very long finite
impulse response (FIR) filters (e.g., around a thousand taps
for 8 kHz sampling) to separate the acoustic signals mixed
under such conditions.
Independent component analysis (ICA) [4], [5] is a major statistical tool for dealing with the BSS problem. If signals are mixed instantaneously, we can directly employ an
instantaneous ICA algorithm to separate them. However,signals are mixed in a convolutive manner in the applications mentioned above. Therefore, we need to extend the
ICA/BSS technique so that it is applicable to convolutive
mixtures.
The first approach is time-domain BSS, where ICA is
directly extended to the convolutive mixture model [6]–[11].
This approach is theoretically sound and achieves good separation once an algorithm converges, since the algorithm
correctly evaluates the independence of separated signals.
Manuscript received February 16, 2005.
Final manuscript received March 9, 2005.
†
The authors are with NTT Communication Science Laboratories, NTT Corporation, Kyoto-fu, 619-0237 Japan.
a) E-mail: maki@cslab.kecl.ntt.co.jp
DOI: 10.1093/ietfec/e88–a.7.1640

However, an ICA algorithm for convolutive mixtures is notas simple as an ICA algorithm for instantaneous mixtures,and is computationally expensive for long FIR filters because it includes convolution operations.
The second approach is frequency-domain BSS, wherecomplex-valued ICA for instantaneous mixtures is employed in each frequency bin [12]–[29]. The merit of thisapproach is that the ICA algorithm remains simple andcan be performed separately at each frequency. Also, any
complex-valued instantaneous ICA algorithm can be employed with this approach. The computational time for BSS
can be reduced by employing a fast algorithm such as FastICA [30], [31], and/or by performing parallel computation
for multiple frequency bins. However, the permutation ambiguity of the ICA solution becomes a serious problem. We
need to align the permutation in each frequency bin so that aseparated signal in the time domain contains frequency components from the same source. This problem is well known
as the permutation problem of frequency-domain BSS [12]–[21], [25]–[27], which is the main focus of this paper. Another problem relates to the circularity eﬀect of discrete frequency representation. Frequency responses calculated in
the frequency domain assume a periodic time-domain filter
for their implementation. However, such a periodic filter
is unrealistic, and we usually use its one-period realization
for the separation filter. Therefore, the frequency responses
should be smoothed so that the one-period realization does
not rely on the circularity eﬀect [18], [29]. This paper also
discusses this problem.
The third approach uses both the time and frequency
domains. In some time-domain BSS methods, convolutions in the time domain are speeded up by the overlap-savemethod in the frequency domain [10], [32]. Furthermore, insome methods [33]–[35], filter coeﬃcients are updated in
the frequency domain while nonlinear functions for evaluating independence are applied in the time domain. The permutation problem does not occur in either case since the independence of separated signals is evaluated in the time domain. Nor does the circularity problem occur when there isan appropriate constraint for filter coeﬃcients [36] by suchmeans as rectangular windowing. However, the algorithm
moves back and forth between the two domains at every
iteration, spending non-negligible time on discrete Fouriertransforms (DFTs) and inverse DFTs. Therefore, we consider that the permutation and circularity problems are inevitable if we hope to benefit from the merits of frequencydomain BSS.

c 2005 The Institute of Electronics, Information and Communication Engineers
Copyright MAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1641

Fig. 2
Fig. 1This paper deals with the second approach, i.e.,frequency-domain BSS. We begin by formulating the BSS
problem for convolutive mixtures in Sect. 2. Section 3 provides an overview of frequency-domain BSS. We thenpresent several important techniques that enable this approach to achieve eﬀective separation of many sources
mixed in a reverberant environment. Section 4 discusses
complex-valued ICA for instantaneous mixtures. Understanding the separation mechanism of BSS in Sect. 5 greatly
helps us to cope with the problem. Section 7 presents amethod for solving the permutation problem, which is themost important technique for frequency-domain BSS. To
solve the permutation problem, information on source location is very useful. This can be estimated from ICA solutions as shown in Sect. 6. The key point with respect tosource localization is that the estimation of the mixing system is easily obtained. This is because the ICA algorithm isjust for instantaneous mixtures, and therefore it is straightforward to calculate the (pseudo)-inverse of a separation matrix, which corresponds to the mixing system. This fact also
makes it easy to solve the scaling ambiguity as shown in
Sect. 8. Section 9 discusses a spectral smoothing techniquedesigned to solve the circularity problem. The experimental results shown in Sect. 10 are very promising. Section 11
concludes this paper.
2.

BSS for Convolutive Mixtures

In the case of audio source separation, several sensor microphones are placed in diﬀerent positions so that each records
a mixture of the original source signals at a slightly diﬀerent
time and level. In the real world, where the source signals
are speech and the mixing system is a room, the signals thatare picked up by the microphones are aﬀected by reverberation. Suppose that N source signals si (t) are mixed andobserved at M sensors
x j (t) =

N 

i=1

Task of blind source separation of speech signals.

BSS system configuration.

h ji (l) si (t − l), j = 1, . . . , M,(1)lwhere h ji (l) represents the impulse response from source i tosensor j. We assume that the number of sources N is known
or can be estimated in some way (e.g., by [37]), and theFig. 3

BSS for convolutive mixtures.

number of sensors M is more than or equal to N (N ≤ M).
The separation system typically consists of a set of FIR
filters wi j (l) of length L to produce N separated signals
yi (t) =

L−1
M 


wi j (l) x j (t − l), i = 1, . . . , N(2)

j=1 l=0

at the outputs. The separation filters are estimated so that the
separated signals become mutually independent. The separation filters wi j (l) should be obtained blindly, i.e., without
knowing si (t) or h ji (l).
A two-input, two-output convolutive BSS problem, i.e.,
N = M = 2, is shown in Figs. 1 and 2. It is assumed that thesource signals s1 and s2 are mutually independent. This assumption usually holds for sounds in the real world. There
are two microphones which pick up the mixed speech. Only
the observed signals x1 and x2 are available, and they arecorrelated. The goal is to adapt the separation systems wi j
and to extract y1 and y2 so that they are mutually independent. With this operation, we can obtain s1 and s2 in the
output y1 and y2 . No information is needed on the source
positions or period of source existence/absence. Nor is any
information on the mixing systems h ji (l) required. Thus,
this task is called blind source separation.
Figure 3 shows a block diagram of BSS. The ideal goalof BSS is to separate and deconvolve the mixtures x j (t), andto obtain a delayed version of source si (t) at each output
i. However, this is very diﬃcult if si (t) is a colored signal,
which is the case when separating natural sounds such as
speech [9]. A practical alternative goal [8], [11] is to obtainthe convolved version of a source si (t) measured at a sensor
Ji :


L 
(3)
h Ji i (l) si t − − l ,yi (t) =2lwhere the sensor index Ji can be selected according to each

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 2005

1642

[a1 , · · · , aN ] = W+ ,ai = [a1i , . . . , a Mi ]T ,Fig. 4is very useful for source localization and scaling alignment,
as described in Sects. 6 and 8, respectively. It should be
noted that it is not diﬃcult to make W invertible by using
an appropriate ICA procedure (for an example, see Sect. 4).
By multiplying both sides of (7) by W+ , the sensor sample
vector x(τ) is represented by a linear combination of basis
vectors a1 , . . . , aN :Flow of frequency-domain BSS.

output i. The way used to attain this goal will be discussed
in Sect. 8.
3.

Overview of Frequency-Domain Approach

Figure 4 shows the flow of frequency-domain BSS. Timedomain signals x j (t) sampled at frequency f s are converted
into frequency-domain time-series signals x j ( f, τ) with an
L-point short-time Fourier transform (STFT):
2 −1

L

x j ( f, τ) =x j (τ + r) win(r) e− 2π f r,

(4)

r=− L2where f ∈ {0, L1 f s , . . . , L−1
L f s } is a frequency, win(r) is awindow that tapers smoothly to zero at each end, such as a
Hanning window 12 (1 + cos 2πr
L ), and τ is a new index representing time.
The remaining operations are performed in the frequency domain. The advantage is that the convolutive mixtures in (1) can be approximated as instantaneous mixturesin each frequency bin:x j ( f, τ) =Nh ji ( f )si ( f, τ),

where h ji ( f ) is the frequency response from source i to sensor j, and si ( f, τ) is a frequency-domain time-series signalof si (t) obtained by the same operation as (4). The vector
notation of the mixing model (5) isNhi ( f )si ( f, τ),

(6)

i=1

where x = [x1 , . . . , x M ]T is a sensor sample vector and hi =[h1i , . . . , h Mi ]T is the vector of the frequency responses fromsource si to all M sensors.
To obtain the frequency responses wi j ( f ) of separation
filters wi j (l) in (2), complex-valued ICAy( f, τ) = W( f )x( f, τ)x( f, τ) =(7)

is solved, where y = [y1 , . . . , yN ]T is a vector of separated
signals, W = [w1 , . . . , wN ]H is an N × M separation matrix,
wi = [wi1 , . . . , wiM ]H , and wi j = [W]i j . The details of the
ICA algorithm are discussed in Sect. 4.
Calculating the Moore-Penrose pseudoinverse W+ ( reduced to the inverse W−1 if N = M) of W asNai ( f )yi ( f, τ).

(10)

i=1

It is well-known that an ICA solution (7) has permutation and scaling ambiguities: even if we permute the rowsof W( f ) or multiply a row by a constant, it is still an ICA
solution. In matrix notation,W( f ) ← Λ( f )P( f )W( f )(11)is also an ICA solution for any permutation P( f ) and diagonal Λ( f ) matrix. Permutation alignment is to decide P( f )so that a time-domain separated signal contains frequency
components from the same source. Section 7 presents amethod for solving this problem. Scaling alignment is todecide Λ( f ) so that a time-domain separated signal satisfies
the goal (3), as discussed in Sect. 8.
Then, we perform spectral smoothing so that a timedomain separation filter tapers smoothly to zero at each end.
This is typically achieved by multiplying the time-domain
filter by a Hanning window, which is equivalent to smoothing the frequency-domain separation matrices asW( f ) ←(5)i=1x( f, τ) =

(8)
(9)

1
[ W( f −∆ f ) + 2W( f ) + W( f +∆ f ) ],4

where ∆ f = Lfs is the diﬀerence from the adjacent frequency. However, this smoothing changes the ICA solution
and causes an error. Section 9 discusses the error and how
to minimize it.
Finally, separation filters wi j (l) are obtained by applying inverse DFT to wi j ( f ) = [W( f )]i j :L
wi j ( f )e 2π f (l− 2 ) ,wi j (l) =f ∈{0, L1 f s , ..., L−1
L fs }

where l = 0, . . . , L−1. The reason for using e 2π f (l− 2 ) insteadof e 2π f l is to make the separation filter wi j (l) causal. Then,the separated signals yi (t) are produced by (2).
L

4.

Complex-Valued ICA

This section discusses how to solve the ICA Eq. (7). Oneof the advantages of frequency-domain BSS is that we can
employ any ICA algorithm for instantaneous mixtures, suchas the information maximization approach (InfoMax) [38]combined with the natural gradient [39], FastICA [30],JADE [40], or an algorithm based on the non-stationarity of

MAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1643

signals [41]. Here, we explain a procedure that was shown
to be eﬃcient by the experiments described in Sect. 10. The
procedure consists of the following three steps:
1. Dimension reduction and whitening by eigenvalue decomposition
2. ICA by a unitary matrix (FastICA)
3. ICA by InfoMax combined with the natural gradient
The first step performs a linear transformationz(τ) = Vx(τ)
for M-dimensional sensor observations x(τ) such that thedimension of z(τ) is reduced (if necessary) to the numberof sources N and z(τ) is spatially whitened (sphered), i.e.,z(τ)z(τ)H τ = I, where I is the N × N identity matrix. The
linear transformation V is typically obtained by eigenvalue
decomposition. Let λ1 ≥ . . . ≥ λ M be sorted eigenvalues of the spatial correlation matrix R = x(τ)x(τ)H τ ande1 , . . . , e M be their corresponding eigenvectors. Then, thelinear transformation is
V = D−1/2 EH ,where D = diag(λ1 , . . . , λN ) is the diagonal matrix of the Nlargest eigenvalues, E = [e1 , . . . , eN ] is the matrix of theircorresponding eigenvectors, and ei = [e1i , . . . , eMi ]T .
This step has practical importance for the following
two reasons. First, the outputs y(τ) of ICA (7) adhere to
the signal subspace that is identified by the N eigenvectors e1 , . . . , eN . This means that the following ICA algorithm does not pursue its solution in the noise subspace,which consequently stabilizes the algorithm and also hasa noise/reverberation reduction eﬀect [18]. A geometricalinterpretation of the dimension reduction is given in [28].
Second, the whitening zzH τ = I is necessary for FastICA,and also provides an eﬃcient convergence for InfoMax even
if the step size is constant over all frequency bins.
The second step performs ICA in a constrained form:y(τ) = Bz(τ),
where B is an N × N unitary matrix: BBH = I. This is performed by a complex-valued version of FastICA [30], [31].
It is very eﬃcient because a fairly good solution can be obtained with only several iterations. The eﬃciency comesfrom the fact that z is whitened and B is unitary. However,there remains room for improving the solution by using another ICA algorithm. One of the reasons is that the output yof FastICA is whitened y( f, τ)y( f, τ)H τ = I and therefore
uncorrelated, whereas original sources s1 ( f, τ), . . . , sN ( f, τ)
are not always completely uncorrelated with a limited number of samples.
The third step improves the ICA solution obtained so
far as an initial valuey(τ) = Wx(τ) = BVx(τ)
by employing another ICA algorithm that does not have theunitary constraint. Based on the use of InfoMax combined

Fig. 5 Directivity patterns (a) obtained by BSS (T R =0 ms), (b) obtained
by BSS (T R =300 ms), (c) obtained by ABF (T R =0 ms), and (d) obtained by
ABF (T R =300 ms).

with the natural gradient, a separation matrix W is gradually
improved by the learning rule:W ← W + µ [I − Φ(y(τ))y(τ)H τ ] W,

(12)

where µ is a step-size parameter. Φ(y) = [Φ(y1 ), . . . ,Φ(yN )]T is an element-wise nonlinear function defined byΦ(yi ) = −∂log p(yi ),
∂yi

(13)

where p(yi ) is the probability density function (pdf) of
a complex-valued signal yi = |yi | e ·arg(yi ) . Since yi is afrequency-domain signal whose phase can be shifted arbitrarily by shifting the STFT window position (4), a feasible assumption is that the pdf is independent of the phasep(yi ) = α · p(|yi |), where α is a constant. This assumption
reduces (13) toΦ(yi ) = ϕ(|yi |) e ·arg(yi ) ,∂log p(|yi |).
ϕ(|yi |) = −∂|yi |(14)(15)

If we assume the Laplacian distribution p(|yi |) = 12 e−|yi | ,which is typical for speech modeling, we have ϕ(|yi |) = 1
and thus a simple nonlinear function
Φ(yi ) = e ·arg(yi ) .
A nonlinear function of the form (14) has a better convergence property [22] than one where the nonlinearity
is applied separately to the real and imaginary parts of acomplex-valued signal yi .
5.

Separation Mechanism of BSS

The mechanism of BSS based on ICA has been shown tobe equivalent to that of an adaptive microphone array system, i.e., N sets of adaptive beamformers (ABFs) with anadaptive null directivity aimed in the direction of unnecessary sounds [23], [24]. From the equivalence between BSS
and ABF, it becomes clear that the physical behavior of BSS
reduces the jammer signal by making a spatial null toward

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 20051644

the jammer, and extract the target.
The separation performance of BSS is compared withthat of ABF. Figure 5 shows the directivity patterns obtained
by BSS and ABF. In Fig. 5, (a) and (b) show directivity patterns by W obtained by BSS, and (c) and (d) show directivity patterns by W obtained by ABF. When T R = 0, a sharp
spatial null is obtained by both BSS and ABF [see Figs. 5(a)and (c)]. When T R = 300 ms, the directivity pattern becomes
duller for both BSS and ABF [see Figs. 5(b) and (d)].
BSS can be regarded as an intelligent version of ABF
in the sense that it can adapt without any information on thesource positions or period of source existence/absence [1].
6.

Source Localization

This section presents a source localization method by analyzing the ICA solution (7) or equivalently (10). The information on source locations can be used to solve the permutation problem, as described in the next section. Many
source localization methods have been proposed. A widelyused method is MUSIC (MUltiple SIgnal Classification)[42], which employs subspace analysis with second-order
statistics. The ICA-based method, on the other hand, employs higher-order statistics (or multiple second-order statistics based on non-stationarity). In this sense, the ICAbased method has certain advantages over the subspacebased method [43].
The source localization technique that employs ICA
is a by-product of research on frequency-domain BSS.
Direction-of-arrival (DOA) estimation methods [19]–[21]have been proposed that are based on beamforming theory[44]. They calculate directivity patterns as shown in Fig. 5from the separation matrix W, and then search the null directions, which correspond to the directions of sources [24].
However, it is simpler and more eﬀective to estimate the directions directly from the basis vectors ai , which are givenby the pseudoinverse of W. The source localization method[25]–[27], [43] presented in this section is based on this idea.
Such an idea was taken for granted in research on blind identification [45], [46], where the mixing system is estimateddirectly.
6.1 Basic Theory of Nearfield Model
Let us assume a mixing model that is suitable for source
localization. Although the mixing model (1) in the time domain is a multi-path mixing model, we approximate the frequency response h ji ( f ) in (5) with a nearfield (direct-path)model (Fig. 6):
h ji ( f ) ≈1−1e 2π f c (||qi −p j ||−||qi ||) ,||qi − p j ||(16)

where p j and qi are 3-dimensional vectors representing thelocations of sensor j and source i, respectively, and c is the
propagation velocity of the signals. We assume that the amplitude is attenuated based on the distance ||qi −p j ||. We also
assume that the phase depends on the diﬀerence between theFig. 6

Nearfield (direct-path) model.

distances ||qi − p j || − ||qi || from the source to the sensor andto the origin o = [0, 0, 0]T . This makes the phase zero atthe origin. If the phase 2π f c−1 (||qi − p j || − ||qi ||) is outside
the range (−π, π), this model suﬀers from spatial aliasing.
Therefore, the model is feasible as long as the condition




c
f <
 2 · ( ||qi − p j || − ||qi || ) is satisfied.
The ICA-based source localization discussed in thissection estimates the location qi of source i from information on sensor locations p j and the separation matrix W( f )obtained by ICA (7). Let us assume here that the decomposition (10) of observations x( f, τ) has been obtained ineach frequency bin by the pseudoinverse of W( f ). By comparing (6) and (10), we observe the following fact. If the
ICA algorithm works well and the outputs y1 , . . . , yN are
the estimation of the sources s1 , . . . , sN , then the basis vectors a1 , . . . , aN are also estimations of the mixing vectors
h1 , . . . , hN up to the permutation and scaling ambiguity.
Following the model (16), the ratio between two elements a ji , a j i of the same basis vector ai provides the key
equation for source localization:a ji
αi h ji=aj i
αi h j i||qi − p j || 2π f c−1 (||qi −p j ||−||qi −p j ||)e=,||qi − p j ||(17)

where the scaling ambiguity αi is canceled out by calculating the ratio. The permutation ambiguity still remains.
However, if we estimate the location qi for all i = 1, . . . , N,
the set of all estimated locations does not depend on the permutation.
With respect to the phase diﬀerences, the set of vectors
qi in the argument of (17),||qi − p j || − ||qi − p j || =arg(a ji /a j i ),2π f c−1(18)defines a surface where the diﬀerence between the distancesfrom p j and p j is constant. The surface is one sheet of a
two-sheet hyperboloid.
Alternatively, with respect to the level diﬀerences, theset of vectors qi in the modulus of (17),||qi − p j ||  a ji 
=
(19)
,
||qi − p j ||  a j i MAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1645

Fig. 7 Source localization by intersection of two hyperboloids and asphere.

Fig. 8

Farfield model.

defines a sphere where the ratio of the distances from p j andp j is constant. Therefore, with these two Eqs. (18) and (19),
we can estimate the possible location qi of source si . Such
hyperboloid and sphere are defined by a pair of sensors j
and j . If we select another pair of sensors, a diﬀerent hyperboloid and sphere are obtained. In this way, the location
qi is estimated as the intersection of several hyperboloids
and spheres. An example is shown in Fig. 7.
6.2 DOA Estimation with Farfield Model
Although it is useful to estimate a 3-dimensional location,
calculating the intersections of hyperboloids and spheres iscomputationally demanding. In many cases it is suﬃcient toestimate just the direction-of-arrival (DOA) of source si . If
we assume the source location qi is far from sensors p j andp j , (18) can be approximated as a farfield model (Fig. 8):(p j − p j )Targ(a ji /a j i )qi=,||qi ||
2π f c−1(20)

and the cosine of angle θij j between the two vectors qi andp j − p j can be calculated ascos θij j(p j − p j )T qi=||p j − p j || · ||qi ||arg(a ji /a j i )=.
2π f c−1 ||p j − p j ||(21)

The set of vectors qi that satisfy (20) represents a coneFig. 9 3-dimensional arrangement of eight microphones and three loudspeakers (upper) and DOA estimation results for this case (lower).

[26], which is the asymptotic surface of the corresponding
hyperboloid (18). To estimate the DOA of a source, the intersections of several cones should be obtained. Let us assume that we select u cones whose corresponding sensor
pairs are ( j1 , j1 ), . . . , ( ju , ju ). The set of Eqs. (20) for u sensor pairs is represented as
ri
qi
=
,
(22)
D
||qi || 2π f c−1
where
D = [p j1 −p j1 , . . . , p ju −p ju ]T,
ri = [ arg(a j1 i /a j1 i ), . . . , arg(a ju i /a ju i ) ]T .
In practical situations, there is no exact solution for (22) because the u conditions do not coincide exactly. Therefore,
we typically solve it in the least-square sense by using theMoore-Penrose pseudoinverse [27]:qi
D+ r i=.
||qi || 2π f c−1(23)

If rank(D) ≥ 3, the set of vectors qi that satisfy (23) represents a line in 3-dimensional space, which represents the
DOA of a source i.
The upper photo in Fig. 9 shows the case where
eight microphones and three loudspeakers are arranged 3dimensionally, and the lower plot shows the DOA estimation results for this case. Each point shows a location vector

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 20051646

qi ( f ) that is normalized to unit norm qi ( f ) ← ||qqii (( ff ))|| . The
estimations are obtained for all frequencies f and all output
indexes i. As shown in the plot, they form clusters, each ofwhich corresponds to the location of each source.
If the sensor and source locations are limited to a 2dimensional plane, the dimensionality of location vectors,
such as pi and qi , can be reduced to two. In this case,rank(D) ≥ 2 is suﬃcient to reach a solution in (23). Moreover, the DOA of source i can be represented simply by theangle θi that satisfies
qi = [cos(θi ), sin(θi )]T ,−180◦ < θi ≤ 180◦ .

(24)

Figure 15 shows the case where the sensor and source locations are limited to 2-dimensions. The DOA estimations inthis case are shown in Figs. 16 and 17.
If the sensors are arranged linearly and the potential
source location is in a 2-dimensional half-plane, which is
to one side of the sensor arrangement line, the angle θij j(0◦ ≤ θij j ≤ 180◦ ) by (21) provides suﬃcient informationon the source location. For example, Fig. 13 shows DOA
estimation results for such a case with the conditions shownin Fig. 12.
7.

from the second and the third categories utilize basically thesame information because the separation matrix W( f ) andthe basis vectors a1 ( f ), . . . , aN ( f ) are directly connected by
the pseudoinverse operation (8). However, the information
used in the third category is easier to handle since it directly
represents the mixing system (6). The last category includes
an approach that employs the inter-frequency correlations of
output signal envelopes [16], [17]. This is particularly eﬀective for a non-stationary signal such as speech.
In the next two subsections, we explain the approachesof the third and the fourth categories, respectively. Since
these two approaches have diﬀerent but complementarycharacteristics, integrating them is a good way to pursue abetter solution to the permutation problem [25]. Subsection
7.3 presents a method that eﬀectively integrates the two approaches to solve the permutation problem in a better way.
In the following subsections, let Π f be a permutation corresponding to the inverse P−1 ( f ) of the permutation matrix of(11). The permutation problem can be formulated to obtainΠ f for every frequency f , which is a mapping from source
index k to output index i:i = Π f (k).

Permutation Alignment
7.1 Localization ApproachThis section discusses how to solve the permutation problem. Various methods have already been proposed. Withreference to the ICA Eq. (7) as well as to the decomposition(10) of observations x( f, τ), we classify these methods into
four categories based on the following strategies:
1. Applying an operation to the separation matrix W( f ),2. Utilizing the information on the separation matrix
W( f ) itself,
3. Utilizing the information on the basis vectors
a1 ( f ), . . . , aN ( f ),4. Utilizing the information on the separated signals
y1 ( f, τ), . . . , yN ( f, τ).
The operation of the first strategy basically involves
smoothing the separation matrices in the frequency domain.
This has been realized by reducing the filter length by rectangular windowing in the time domain [10], [13]–[15], or by
averaging the separation matrices with adjacent frequencies[13]. However, this operation makes the separation matrix
W( f ) diﬀerent from the ICA solution (7), which may havea detrimental eﬀect on the separation performance. A possible way to solve this problem is to interleave the ICA update, e.g., (12), and this operation until convergence. In thissense, this strategy is related to the third approach to BSS
discussed in the Introduction.
The second category includes the beamforming approach [19]–[21], where the directivity patterns formed by
the separation matrix are analyzed to identify the DOA ofeach source. The third category includes an approach thatutilizes the results of source localization with the basis vectors [25]–[27], [46]. The theory and operation for source localization were discussed in Sect. 6. These two approaches

The basic idea of this approach is to estimate the locationsof sources and then cluster them to decide the permutation.
ICA-based source localization (Sect. 6) estimates the location qi ( f ) of a source that corresponds to the i-th basis vector ai ( f ) for each frequency f . Let the following function
localize estimate the location in this way:qi ( f ) = localize( f, ai ( f )) .
If just the DOA estimation is adequate, the location vectorqi ( f ) should be normalized to the unit norm qi ( f ) ← ||qqii (( ff ))|| .
If the locations of sensors and sources are limited to a 2dimensional plane, we simply obtain θi ( f ) that satisfies (24)as a DOA estimation.
Then, we employ a clustering algorithm to find N clusters C1 , . . . , C N formed by estimated locations qi ( f ) or θi ( f ).
Each Ck corresponds to the location of source k. Let the following function clustering perform clustering for all of theestimated locations qi ( f ) and return the centroid ck and thevariance σ2k of each cluster Ck :[c1 , σ1 , . . . , cN , σN ]= clustering (∀ f, q1 ( f ), . . . , qN ( f )),
 q
ck =
,
|Ck |
q∈C
k

σ2k =

 ||ck − q||2
q∈Ck

|Ck |

,

where |Ck | is the number of vectors in the cluster. The optimization
criterion for clustering is to minimize the total sum
N
2
σ
of
the variances. This optimization is eﬃciently
k=1 kMAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1647

and v1566, as well as v1562
and v1566, are highly correlated.
1
2
2
Thus, calculating such correlations helps us to align permutations.
A simple criterion for deciding Π f is to maximize thesum of the correlations between neighboring frequencies
within distance δ:N fcor(vΠ(i), vgΠg (i) ),(26)
Π f = argmaxΠ
|g− f |≤δ i=1

Fig. 10

Envelopes of two output signals at diﬀerent frequencies.

performed with the k-means clustering algorithm [47]. Once
we have N clusters, permutations for all frequencies f canbe decided byΠ f = argminΠ

N||ck − qΠ(k) ( f )||2 .

(25)

k=1

The advantage of this source localization approach isthat it is very simple to decide the permutation Π f for eachfrequency once the centroids of N clusters are obtained.
However, the downside of this approach is that the estimatedlocations or DOAs, and thus the permutations Π f are not accurate for some frequencies. Such situations typically happen at low frequencies, where the phase diﬀerence causedby the sensor spacing is very small, as shown in Fig. 13.
7.2 Correlation Approach

where Πg is the permutation at frequency g. This criterion
is based on local information and has a drawback in that
mistakes in a narrow range of frequencies may lead to thecomplete misalignment of the frequencies beyond the range.
To avoid this problem, the method in [17] does not limit
the frequency range in which correlations are calculated. It
decides permutations one by one based on the criterion


N




f
cor vΠ(i),vgΠg (i)  ,Π f = argmaxΠ
i=1

where F is a set of frequencies in which the permutation
is decided. This method assumes high correlations of envelopes even between frequencies that are not close neighbors. This assumption is not satisfied for all pairs of freand v3516
in Fig. 10 do not have a highquencies, e.g., v1566
ii
correlation. Therefore, this method still has the drawbackof permutations possibly being misaligned at many frequencies.
If a source signal has a harmonic structure, as in thecase of speech, there are strong correlations between theenvelopes of a fundamental frequency f and its harmonics
2 f, 3 f, . . .. Therefore, maximizing the correlation among
harmonics is another idea for permutation alignment [25]:Π f = argmaxΠN g∈H( f ) i=1

This subsection presents an approach to permutation alignment based on the inter-frequency correlation of separated
signals. The correlation should be calculated for the amplitude |yi ( f, τ)| or (log-scaled) power |yi ( f, τ)|2 of separated signals. The correlation of raw complex-valued signals yi ( f, τ) would be very low due to the STFT property.
Here, we use the amplitude (so-called envelope)vif (τ) = |yi ( f, τ)|of a separated signal yi ( f, τ). The correlation of two sequences x(τ) and y(τ) is usually calculated by the correlationcoeﬃcientcor(x, y) = (µ x·y − µ x · µy )/(σ x · σy ),where µ x is the mean and σ x is the standard deviation of x.
Based on this definition, cor(x, x) = 1, and cor(x, y) = 0 if x
and y are uncorrelated.
Envelopes have high correlations at neighboring frequencies if separated signals correspond to the same source
signal. Figure 10 shows an example. Two envelopes v1562
1

g∈F

f
cor(vΠ(i), vgΠg (i) ),(27)

where H( f ) provides a set of harmonic frequencies of f .
The permutation accuracy improves if we take the signal’s
harmonic structure into consideration. However, maximizing (26) and (27) simultaneously is not very straightforward
and is computationally expensive.
7.3 Integrated Method
This subsection presents a method that integrates the two approaches discussed in the last two subsections. The intention
behind this integration is to solve the permutation problem
robustly and precisely. Let us review the characteristics ofthe above two approaches.
• robustness: The localization approach is robust sincea misalignment at one frequency does not aﬀect otherfrequencies. The correlation approach is not robust
since a misalignment at one frequency aﬀects the results of other frequencies and may cause consecutive
misalignments.
• preciseness: The localization approach is not precise

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 2005

1648

since the evaluation is based on a direct-path approximation (16) of the mixing system. The correlation approach is precise as long as signals are well separated
by ICA, since the measurement is based on the separated signals themselves.
To benefit from both advantages, namely the robustness of
the localization approach and the preciseness of the correlation approach, the integrated method first decides permutations with the localization approach and then refines thesolution with the correlation approach. An implementationof the integrated method consists of the following four steps[25]:
1. Decide the permutations by the localization approach(25) at certain frequencies where the confidence ofsource localization is suﬃciently high,
2. Decide the permutations based on neighboring correlations (26) as long as the criterion gives a clear-cut
decision,
3. Decide the permutations at certain frequencies where
the correlation among harmonics (27) is suﬃciently
high,
4. Decide the permutations for the remaining frequencies
based on neighboring correlations (26).
The key to the first step is fixing a permutation only if
the confidence of source localization is suﬃciently high. We
assume that the confidence is high if the squared distance between an estimated location and its corresponding centroidis smaller than the variance, i.e., ||ck − qΠ(k) ( f )||2 < σ2k . In
the second step, permutations are decided one by one for thefrequency f where the sum of the correlations with fixed frequencies g ∈ F within distance |g − f | ≤ δ is the maximum.
This is repeated as long as the maximum correlation sum islarger than a threshold thcor . In the third step, the permutations are decided for frequencies f where the sum of thecorrelations among harmonics is larger than a threshold thha .
The last step decides the permutations for the remaining frequencies with the same criterion as the second step.
Let us discuss the advantages of the integrated method.
The main advantage is that it does not cause a large misalignment as long as the permutations fixed by the localization approach are correct. Moreover, the correlation partcompensates for the lack of preciseness of the localization
approach. The correlation part consists of three steps (step
2,3,4) for two reasons. First, the harmonics part works well
if most of the other permutations are fixed. Second, themethod becomes more robust by quitting step 2 if there isno clear-cut decision. With this structure, we can avoid
fixing the permutations for consecutive frequencies without high confidence. As shown in the experimental results(Sect. 10), this integrated method is eﬀective in separating
many sources.
8.

Scaling Alignment

The scaling ambiguity Λ( f ) in (11) is easily solved by calculating the (pseudo)-inverse of a separation matrix W( f ) [8],Fig. 11 Impulse responses uik (l) obtained with periodic filters (above)
and with their one-period realization (below).

[17]. The frequency-domain counterpart of the BSS goal (3)isyi ( f, τ) = h Ji i ( f )si ( f, τ),(28)

where Ji can be selected according to each output i but
should be the same for all frequencies f . Let us assume
that the ICA and the permutation problem have been solved.
Then the ai term in (10) is close to the hi term in (6):hi ( f )si ( f, τ) ≈ ai ( f )yi ( f, τ).

(29)

By substituting (28) into (29), we have the condition for
scaling alignment:hi ( f ) ≈ ai ( f )h Ji i ( f ) ⇔ a Ji i ( f ) ≈ 1.
This condition, i.e., a Ji i ( f ) = 1, is attained byW( f ) ← Λ( f )W( f ),Λ( f ) = diag(a J1 1 ( f ), . . . , a JN N ( f )),where a ji ( f ) = [W+ ( f )] ji is an element of the pseudoinverseof W( f ).
9.

Spectral Smoothing

The frequency-domain BSS described in this paper is influenced by the circularity of discrete frequency representation.
The circularity refers to the fact that frequency responses
sampled at L points with an interval f s /L ( f s : sampling frequency) represent a periodic time-domain signal whose period is L/ f s . Since this filter is unrealistic, we usually use
its one-period realization. However, such one-period filters
may cause a problem. Figure 11 shows impulse responses
from a source sk (t) to an output yi (t) defined by (39). Responses on the left u11 (l) correspond to the extraction of atarget signal, and those on the right u14 (l) correspond to the
suppression of an interference signal. The upper responses
are obtained with infinite-length filters, and the lower ones
with one-period filters. We can see that the one-period filters create spikes, which distort the target signal and degradeMAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1649

the separation performance. Note that these spikes are inevitable in the frequency-domain BSS since we have an ICA
solution in the frequency domain.

wi ( f ) and wi ( f −∆ f ) (or wi ( f +∆ f )). Since these diﬀerences
are usually not very large, the error ei does not seriously affect the separation if we use a Hanning window for spectral
smoothing.

9.1 Windowing
9.2 Minimizing Error by Adjusting Scaling Ambiguity
To solve this problem, we need to control the frequency responses wi j ( f ) so that the corresponding time-domain filter
wi j (l) does not rely on the circularity eﬀect whereby adjacent periods work together to perform some filtering. The
most widely used approach is spectral smoothing, which isrealized by multiplying a window g(l) that tapers smoothly
to zero at each end, such as a Hanning window g(l) =12πl2 (1 + cos L ). This makes the resulting time-domain filter wi j (l) · g(l) fit length L and have small amplitude aroundthe ends [18]. As a result, the frequency responses wi j ( f )are smoothed asw̃i j ( f ) =fs −∆ fg(φ)wi j ( f − φ),

φ=0

where g( f ) is the frequency response of g(l) and ∆ f = Lfs .
If a Hanning window is used, the frequency responses aresmoothed asw̃i j ( f ) =1[wi j ( f −∆ f ) + 2wi j ( f ) + wi j ( f +∆ f )], (30)4since the frequency responses g( f ) of the Hanning window
are g(0) = 12 , g(∆ f ) = g( f s −∆ f ) = 14 , and zero for the otherfrequency bins.
The windowing successfully eliminates the spikes.
However, it changes the frequency response from wi j ( f ) tow̃i j ( f ) and causes an error. Let us evaluate the error foreach row wi ( f ) = [wi1 ( f ), . . . , wiM ( f )]T of the ICA solution
W( f ). The error isei ( f ) = min[w̃i ( f ) − αi wi ( f )]Hw̃i ( f ) wi ( f )wi ( f ),||wi ( f )||2

1 −[ e ( f ) + e+i ( f ) ],
4 i

(32)

where
wi ( f −∆ f )H wi ( f )wi ( f ),e−i ( f ) = wi ( f −∆ f ) −||wi ( f )||2
wi ( f +∆ f )H wi ( f )wi ( f ).
e+i ( f ) = wi ( f +∆ f ) −||wi ( f )||2(35)

We want to find di ( f ) such that the error (31) is minimized.
The scalar di ( f ) should be close to 1 to avoid any greatchange in the predetermined scaling. Thus, an appropriate
total cost to be minimized is||ei ( f )||2
Ji ( f ), Ji ( f ) =+ β|di ( f ) − 1|2 ,J=||wi ( f )||2fwhere β is a parameter indicating the importance of maintaining the predetermined scaling. With the Hanning window, the error after the scaling adjustment is easily calculated by substituting (35) for (32):ei ( f ) =1[ di ( f −∆ f )e−i ( f ) + di ( f +∆ f )e+i ( f ) ],4

(36)

where e−i and e+i are defined in (33) and (34), respectively.
The minimization of the total cost can be performed
iteratively bydi ( f ) = di ( f ) − µ(31)where w̃i ( f ) = [w̃i1 ( f ), . . . , w̃iM ( f )]T and αi is a complexvalued scalar representing the scaling ambiguity of the ICA
solution. The minimization minαi is based on least-squares,
and can be represented by the projection of w̃i to wi . We canevaluate the error for the Hanning window case by substituting (30) for w̃ of (31):ei ( f ) =wi ( f ) ← di ( f )wi ( f ).

∂J∂di ( f )(37)

with a small step-size µ. With the Hanning window, thegradient is

αi= w̃i ( f ) −Even if the error caused by the windowing is not very large,
the separation performance is improved by minimizing theerror [29]. The minimization is performed by adjusting thescaling ambiguity of the ICA solution before the windowing. Let di ( f ) be a complex-valued scalar for the scaling
adjustment:∂J∂Ji ( f −∆ f ) ∂Ji ( f +∆ f ) ∂Ji ( f )=++∂di ( f )∂di ( f )∂di ( f )∂di ( f )H +ei ( f −∆ f ) ei ( f −∆ f ) + ei ( f +∆ f )H e−i ( f +∆ f )=8 · ||wi ( f )||2+2β(di ( f ) − 1).
(38)
With equations from (36) to (38), we can optimize the scalardi ( f ) for the scaling adjustment, and minimize the error
caused by the spectral smoothing (30) with the Hanningwindow.
10.

(33)
(34)

This e−i (or e+i ) represents the diﬀerence between two vectors

Experimental Results

The performance of BSS is evaluated by a signal-tointerference ratio (SIR), which is the power ratio between
the target component and the interference components. Let
uik (l) be the impulse responses from source sk (t) to separated
signal yi (t):IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 2005

1650
Table 1

Separation performance with linear array.

#sources / position
Spectral smoothing
Average SIR at microphones (dB)
Average SIR of output (dB)
Execution time (s)Fig. 12

2/ac
no
yes
0.1
20.1
22.3
5.2
5.2

3/abd
no
yes
-2.9
14.7
17.0
8.0
8.1

4/abcd
no
yes
-4.6
9.3
11.5
12.3
12.4

Experimental conditions with linear array.

Fig. 14 Comparison of diﬀerent methods for solving permutation
problem.

Fig. 13uik (l) =DOA estimations by (21) with four sources.

L−1
M 


wi j (τ)h jk (l − τ).

(39)

j=1 τ=0

Then, the SIR of output i is calculated as| uii (l)si (t − l)|2 t(dB),
SIRi = 10 log10  l| ki l uik (l)sk (t − l)|2 t(40)

where ·t denotes the averaging operator over time t.
10.1 Linear Array
We performed experiments to separate speech signals in anenvironment whose conditions are summarized in Fig. 12.
Our experiments involved two, three and four sources whose
locations are indicated in Table 1. The sensors were arranged linearly, and the number of sensors used was thesame as the number of sources. We used filters of length
L = 2048 because this length provided the best performance
under the conditions. The BSS program was coded in Matlab and run on Athlon XP 3200+.

The results shown in Table 1 are the average SIRs of
output for eight combinations of 7-second speeches. We cansee that the spectral smoothing discussed in Sect. 9 improves
the average SIR for every setup. The short execution time,
as shown in Table 1, enables the BSS system to perform inreal time if the number of source signals is not very large.
Figure 13 shows DOA estimations for mixtures of four
sources obtained with (21). Figure 14 shows SIRs for threeand four sources with the diﬀerent methods for solving the
permutation problem discussed in Sect. 7. Here, “Localization” is the localization (DOA) approach (25) alone, “Correlation” is the correlation approach (26) alone, “Integrated”
is the integrated method, and “Optimal” is the optimal solution obtained by utilizing the si (t) and h ji (l) information.
The performance of “Localization” was stable but insuﬃcient. The performance of “Correlation” was unstable andvery poor in the four-source cases. The “Integrated” method
performed very well and was close to “Optimal.”
10.2 Planar Array
Next, we carried out experiments on separating six sources
with a planar array of eight microphones. The room layout and other experimental conditions are shown in Fig. 15.
All six sources were 6-second speech signals, and two came

MAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1651

Fig. 15

Experimental conditions for planar array case.

Fig. 17 Permutation solved by using estimated DOAs (upper) andspheres (lower).

Fig. 16
pairs.

Histogram of DOAs estimated with small spacing microphonefrom the same direction. The filter length was again L =2048 for an 8-kHz sampling rate.
Let us explain the method for solving the permutation
problem in this situation. First, the source directions were
estimated with small-spacing microphone pairs (1-3, 2-4, 12 and 2-3 shown in the right-top corner of Fig. 15). This was
performed based on (20), (22) and (23). Figure 16 shows ahistogram of the estimated DOAs. There are five clusters inthis histogram, and one cluster is twice the size of the others.
This implies that two sources came from the same direction(about 150◦ ). We solved the permutation problem for theother four sources by using this DOA information as shownon the upper plot of Fig. 17.
Then, to distinguish between the two sources that camefrom the same direction, the spheres of these sources were
estimated with large-spacing microphone pairs (7-5, 7-8, 65 and 6-8 shown in the center of Fig. 15). This was performed based on (19). The lower plot of Fig. 17 shows theradiuses of the spheres estimated with microphone pair 75. Although the radius estimations had large variances, it
provided suﬃcient information to distinguish between thetwo sources. Consequently, the signal components of all
frequencies were classified into six clusters. We decided the
permutation only for frequency bins where the classification
was reliable, as discussed in Sect. 7.3.

To show the eﬀectiveness of this method, we compared
SIRs by three diﬀerent methods for the permutation problem. Table 2 shows the results. The last row, “DOA +Sphere + Correlation,” shows the results obtained with the
integrated method. The two methods for comparison were
“Correlation” where only the correlations (26) were maximized, and “DOA + Correlation” where only the DOA information was used for the source localization step in the integrated method. To see how much the SIRs were improved,
we also measured the SIR of the mixture observed at microphone 1 (“SIR at microphone 1”). The eﬀectiveness ofthe two integrated methods can again be observed. If wecompare the results of “DOA + Correlation” with “DOA +Sphere + Correlation,” the improvement of the latter overthe former is apparent for sources 4 and 5, which came fromthe same direction. This means that the sphere information
was important in terms of distinguishing between sources
coming from the same direction. The BSS program wasagain coded in Matlab and run on Athlon XP 3200+. The
computational time for separating six speeches of 6 seconds
was around one minute.
10.3 Moving Sources
In most realistic applications, the source location may
change. A mixing system is time-varying when source signals move. A naive approach for tracking a time-varying
system is an online algorithm that updates the separation
system sample by sample [48], [49].

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 2005

1652
Table 2

Separation performance with planar array measured by SIR (dB).

SIR at microphone 1
Correlation
DOA + Correlation
DOA + Sphere + Correlation

SIR1 SIR2 SIR3 SIR4 SIR5 SIR6 average-8.3 -6.8 -7.8 -7.7 -6.7 -5.2
-7.1
4.4 2.6 4.0 9.2 3.6 -2.0
3.7
9.6 9.3 14.7 2.7 6.5 14.0
9.4
10.8 10.4 14.5 7.0 11.0 12.2
11.0

Table 3
CommonICA part

Fig. 18

Experimental conditions.

Sampling rate = 8 kHz
Window = hanning
Reverberation time T R =130 ms
Frame length T ICA = 1024 points (128 ms)
Frame shift = 256 points (32 ms)
g = 100.0µ = optimized for block size T b
Number of iterations NI = 100

Average and standard deviation of SIR for fixed sources.

Fig. 20 SIR of blockwise batch algorithm without postprocessing. Target and interference signals moved at 10 s (T b = 1.0 s).

Fig. 19

Layout of room used in experiments. T R = 130 ms.

Indeed, an online algorithm can track a time-varying
system; however, its performance is generally worse thana batch algorithm, which can employ a number of samples,when the system is stationary. Although we are dealing withmoving sources, we do not want to degrade the performancefor fixed sources.
In this section, we describe a real-time BSS method[50] that employs frequency domain ICA with a blockwise
batch algorithm. This algorithm achieves better separation
performance than an online algorithm for fixed source signals.
We measured the BSS performance using ICA. Figure 18 shows the average and standard deviation of SIR for
fixed sources (the target is at A and the interference at C in
Fig. 19). This indicates that the blockwise batch algorithm
outperforms the online algorithm (in which µ is tuned to optimize the performance) when we use the update Eq. (12).
In addition, the deviation of the batch algorithm is smaller
than that of the online algorithm, which is why we adopt theblockwise batch algorithm. We used block size T b = 1.0 s in

the experiments.
We carried out experiments using speech signals
recorded in a room. The reverberation time of the room was
130 ms. We used two omni-directional microphones withan inter-element spacing of 4 cm. The layout of the room isshown in Fig. 19. The target source signal was first locatedat A and then moved to B at a speed of 30 deg/s. The interference signal was located at C and moved to D at a speedof 40 deg/s.
The step size parameter µ in (12) aﬀects the separation
performance of BSS when the block size changes. We carried out preliminary experiments and chose µ to optimize theperformance for each block size. The other conditions are
summarized in Table 3. We measured SIRs with 30 combinations of source signals using three male and three female
speakers, and averaged them.
We investigate the BSS performance for moving
sources using the blockwise batch algorithm. Figure 20
shows the SIR for a moving target (solid line) and that for amoving interference (dotted line). We can see that the SIR
is not degraded even when the target moves. By contrast,
interference movement causes a decline in the SIR.
This can be explained by the directivity pattern of
the separation system obtained by ICA. The solution of
frequency-domain BSS works in the same way as an adap-

MAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1653

tive beamformer, which forms a spatial null toward an interference signal (Fig. 5). Because of this characteristic, BSS
using ICA is robust with a moving target signal but fragile
with a moving interference signal. Taking advantage of thisnature, we can estimate residual crosstalk components evenwhen the interference signal moves by employing postprocessing in the second stage [50].
11.

Conclusion

This paper presented a comprehensive description of
frequency-domain BSS as well as various techniques thatenable frequency-domain BSS to be used for separating
many speech signals mixed in a real-room environment. The
permutation problem has been a major concern with the frequency domain approach. However, with the methods described in Sect. 7, this problem can be solved even in a practical situation. Moreover, the locations of sources can be
estimated by the method described in Sect. 6. This ability is unique to the frequency domain approach, and cannot be seen in time-domain BSS. Our experimental results
show that the separation performance was fairly good andthe computational cost was feasible. These results demonstrate the eﬀectiveness of frequency-domain BSS.
References[1] S. Makino, “Blind source separation of convolutive mixtures ofspeech,” in Adaptive Signal Processing: Applications to Real-World
Problems, eds. J. Benesty and Y. Huang, Springer, 2003.
[2] S. Haykin, ed., Unsupervised Adaptive Filtering (Volume I: Blind
Source Separation), John Wiley & Sons, 2000.
[3] A. Cichocki and S. Amari, Adaptive Blind Signal and Image Processing, John Wiley & Sons, 2002.
[4] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component
Analysis, John Wiley & Sons, 2001.
[5] T.W. Lee, Independent Component Analysis—Theory and Applications, Kluwer Academic Publishers, 1998.
[6] S. Amari, S. Douglas, A. Cichocki, and H. Yang, “Multichannel
blind deconvolution and equalization using the natural gradient,”
Proc. IEEE Workshop on Signal Processing Advances in WirelessCommunications, pp.101–104, April 1997.
[7] M. Kawamoto, K. Matsuoka, and N. Ohnishi, “A method of blindseparation for convolved non-stationary signals,” Neurocomputing,vol.22, pp.157–171, 1998.
[8] K. Matsuoka and S. Nakashima, “Minimal distortion principle forblind source separation,” Proc. ICA 2001, pp.722–727, Dec. 2001.
[9] S.C. Douglas and X. Sun, “Convolutive blind separation of speech
mixtures using the natural gradient,” Speech Commun., vol.39,
pp.65–78, 2003.
[10] H. Buchner, R. Aichner, and W. Kellermann, “Blind source separation for convolutive mixtures: A unified treatment,” in Audio Signal
Processing for Next-Generation Multimedia Communication Systems, eds. Y. Huang and J. Benesty, Kluwer Academic Publishers,
pp.255–293, 2004,[11] T. Takatani, T. Nishikawa, H. Saruwatari, and K. Shikano, “Highfidelity blind separation of acoustic signals using SIMO-modelbased independent component analysis,” IEICE Trans. Fundamentals, vol.E87-A, no.8, pp.2063–2072, Aug. 2004.
[12] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Frequency-domain
blind source separation,” in Speech Enhancement, eds. J. Benesty,
S. Makino, and J. Chen, Springer, 2005.

[13] P. Smaragdis, “Blind separation of convolved mixtures in the frequency domain,” Neurocomputing, vol.22, pp.21–34, 1998.
[14] L. Parra and C. Spence, “Convolutive blind separation of nonstationary sources,” IEEE Trans. Speech Audio Process., vol.8, no.3,
pp.320–327, May 2000.
[15] L. Schobben and W. Sommen, “A frequency domain blind signal
separation method based on decorrelation,” IEEE Trans. Signal Process., vol.50, no.8, pp.1855–1865, Aug. 2002.
[16] J. Anemüller and B. Kollmeier, “Amplitude modulation decorrelation for convolutive blind source separation,” Proc. ICA 2000,
pp.215–220, June 2000.
[17] N. Murata, S. Ikeda, and A. Ziehe, “An approach to blind source
separation based on temporal structure of speech signals,” Neurocomputing, vol.41, no.1-4, pp.1–24, Oct. 2001.
[18] F. Asano, S. Ikeda, M. Ogawa, H. Asoh, and N. Kitawaki, “Combined approach of array processing and independent component
analysis for blind separation of acoustic signals,” IEEE Trans.
Speech Audio Process., vol.11, no.3, pp.204–215, May 2003.
[19] S. Kurita, H. Saruwatari, S. Kajita, K. Takeda, and F. Itakura, “Evaluation of blind signal separation method using directivity patternunder reverberant conditions,” Proc. ICASSP 2000, pp.3140–3143,
June 2000.
[20] H. Saruwatari, S. Kurita, K. Takeda, F. Itakura, T. Nishikawa, andK. Shikano, “Blind source separation combining independent component analysis and beamforming,” EURASIP Journal on Applied
Signal Processing, no.11, pp.1135–1146, 2003.
[21] M.Z. Ikram and D.R. Morgan, “A beamforming approach to permutation alignment for multichannel frequency-domain blind speech
separation,” Proc. ICASSP 2002, pp.881–884, May 2002.
[22] 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,
March 2003.
[23] S. Araki, R. Mukai, S. Makino, T. Nishikawa, and H. Saruwatari,“The fundamental limitation of frequency domain blind source separation for convolutive mixtures of speech,” IEEE Trans. Speech Audio Process., vol.11, no.2, pp.109–116, March 2003.
[24] S. Araki, S. Makino, Y. Hinamoto, R. Mukai, T. Nishikawa,and H. Saruwatari, “Equivalence between frequency domain blindsource separation and frequency domain adaptive beamforming for
convolutive mixtures,” EURASIP Journal on Applied Signal Processing, vol.2003, no.11, pp.1157–1166, 2003.
[25] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust andprecise method for solving the permutation problem of frequencydomain blind source separation,” IEEE Trans. Speech Audio Process., vol.12, pp.530–538, Sept. 2004.
[26] R. Mukai, H. Sawada, S. Araki, and S. Makino, “Frequency domain
blind source separation using small and large spacing sensor pairs,”
Proc. ISCAS 2004, vol.V, pp.1–4, May 2004.
[27] R. Mukai, H. Sawada, S. Araki, and S. Makino, “Frequency domain
blind source separation for many speech signals,” Proc. ICA 2004
(LNCS 3195), pp.461–469, Sept. 2004.
[28] S. Winter, H. Sawada, and S. Makino, “Geometrical understandingof the PCA subspace method for overdetermined blind source separation,” Proc. ICASSP 2003, pp.769–772, April 2003.
[29] H. Sawada, R. Mukai, S. de la Kethulle, S. Araki, and S. Makino,“Spectral smoothing for frequency-domain blind source separation,”
Proc. IWAENC2003, pp.311–314, Sept. 2003.
[30] A. Hyvärinen, “Fast and robust fixed-point algorithm for independent component analysis,” IEEE Trans. Neural Netw., vol.10, no.3,
pp.626–634, 1999.
[31] E. Bingham and A. Hyvärinen, “A fast fixed-point algorithm for independent component analysis of complex valued signals,” Int. J.
Neural Systems, vol.10, no.1, pp.1–8, Feb. 2000.
[32] M. Joho and P. Schniter, “Frequency domain realization of a multichannel blind deconvolution algorithm based on the natural gradient,” Proc. ICA2003, pp.543–548, April 2003.

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.7 JULY 2005

1654

[33] A.D. Back and A.C. Tsoi, “Blind deconvolution of signals using acomplex recurrent network,” Proc. Neural Networks for Signal Process., pp.565–574, 1994.
[34] R.H. Lambert and A.J. Bell, “Blind separation of multiple speakers
in a multipath environment,” Proc. ICASSP’97, pp.423–426, April
1997.
[35] T.W. Lee, A.J. Bell, and R. Orglmeister, “Blind source separation ofreal world signals,” Proc. ICNN, pp.2129–2135, June 1997.
[36] J.J. Shynk, “Frequency-domain and multirate adaptive filtering,”
IEEE Signal Process. Mag., vol.9, no.1, pp.14–37, Jan. 1992.
[37] H. Sawada, S. Winter, R. Mukai, S. Araki, and S. Makino, “Estimating the number of sources for frequency-domain blind source
separation,” Proc. ICA 2004 (LNCS 3195), pp.610–617, Sept. 2004,[38] A. Bell and T. Sejnowski, “An information-maximization approachto blind separation and blind deconvolution,” Neural Comput., vol.7,
no.6, pp.1129–1159, 1995.
[39] S. Amari, “Natural gradient works eﬃciently in learning,” NeuralComput., vol.10, no.2, pp.251–276, 1998.
[40] J.F. Cardoso, “Blind beamforming for non-Gaussian signals,” IEE
Proceedings-F, pp.362–370, Dec. 1993.
[41] K. Matsuoka, M. Ohya, and M. Kawamoto, “A neural net forblind separation of nonstationary signals,” Neural Netw., vol.8, no.3,
pp.411–419, 1995.
[42] R.O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag., vol.34, pp.276–280, March
1986.
[43] H. Sawada, R. Mukai, and S. Makino, “Direction of arrival estimation for multiple source signals using independent component analysis,” Proc. International Symposium on Signal Processing and its
Applications, pp.411–414, July 2003.
[44] B.D. Van Veen and K.M. Buckley, “Beamforming: A versatile approach to spatial filtering,” IEEE ASSP Mag., vol.5, pp.4–24, April
1988.
[45] J.-F. Cardoso, “Source separation using higher order moments,”
Proc. ICASSP’89, vol.4, pp.2109–2112, May 1989.
[46] V.C. Soon, L. Tong, Y.F. Huang, and R. Liu, “A robust method for
wideband signal separation,” Proc. ISCAS’93, vol.1, pp.703–706,
May 1993.
[47] R.O. Duda, P.E. Hart, and D.G. Stork, Pattern Classification, 2nd ed.,Wiley Interscience, 2000.
[48] J. Anemüller and T. Gramss, “On-line blind separation of moving
sound sources,” Proc. ICA’99, pp.331–334, 1999.
[49] A. Koutras, E. Dermatas, and G. Kokkinakis, “Blind speech separation of moving speakers in real reverberant environment,” Proc.
ICASSP 2000, pp.1133–1136, 2000.
[50] R. Mukai, H. Sawada, S. Araki, and S. Makino, “Blind source separation for moving speech signals using blockwise ICA and residualcrosstalk subtraction,” IEICE Trans. Fundamentals, vol.E87-A, no.8,
pp.1941–1948, Aug. 2004.

Shoji Makino
received the B.E., M.E., andPh.D. degrees from Tohoku University, Japan, in
1979, 1981, and 1993, respectively. He joined
NTT in 1981. He is now an Executive Manager at the NTT Communication Science Laboratories. He is also a Guest Professor at the
Hokkaido University. His research interests include adaptive filtering technologies and realization of acoustic echo cancellation, blind source
separation of convolutive mixtures of speech.
He received the TELECOM System Technology
Award of the TAF in 2004, the Best Paper Award of the IWAENC in 2003,
the Paper Award of the IEICE in 2005 and 2002, the Paper Award of the
ASJ in 2005 and 2002, the Achievement Award of the IEICE in 1997, andthe Outstanding Technological Development Award of the ASJ in 1995. He
is the author or co-author of more than 200 articles in journals and conference proceedings and has been responsible for more than 150 patents. He
is a member of the Conference Board of the IEEE SP Society and an Associate Editor of the IEEE Transactions on Speech and Audio Processing.
He is also an Associate Editor of the EURASIP Journal on Applied Signal Processing. He is a member of the Technical Committee on Audio andElectroacoustics of the IEEE SP Society as well as the Technical Committee on Blind Signal Processing of the IEEE Circuits and Systems Society.
He is also a member of the International ICA Steering Committee and theOrganizing Chair of the ICA2003 in Nara. He is the General Chair of the
IWAENC2003 in Kyoto. He was a Vice Chair of the Technical Committee
on Engineering Acoustics of the IEICE and the ASJ. He is an IEEE Fellow,
a council member of the ASJ, and a member of the EURASIP.

Hiroshi Sawada
received the B.E., M.E.
and Ph.D. degrees in information science from
Kyoto University, Kyoto, Japan, in 1991, 1993
and 2001, respectively. In 1993, he joined NTT
Communication Science Laboratories. From1993 to 2000, he was engaged in research on thecomputer aided design of digital systems, logic
synthesis, and computer architecture. Since
2000, he has been engaged in research on array signal processing, blind source separation
for convolutive mixtures, and speech enhancement. He received the best paper award of the IEEE Circuit and System
Society in 2000. He is a senior member of the IEEE and a member of the
ASJ.

Ryo Mukai
received the B.S. and the M.S.
degrees in information science from the University of Tokyo, Japan, in 1990 and 1992, respectively. He joined NTT in 1992. From 1992 to
2000, he was engaged in research and development of processor architecture for network service systems and distributed network systems.
Since 2000, he has been with NTT Communication Science Laboratories, where he is engaged
in research of blind source separation. His current research interests include digital signal processing and its applications. He received the Sato Paper Award of the ASJ
in 2005. He is a member of the Technical Committee on Blind Signal Processing of the IEEE CAS Society. He is a senior member of the IEEE, amember of ACM, the ASJ, and the IPSJ.

MAKINO et al.: BLIND SOURCE SEPARATION OF CONVOLUTIVE MIXTURES OF SPEECH IN FREQUENCY DOMAIN

1655

Shoko Araki
received the B.E. and the
M.E. degrees in mathematical engineering andinformation physics from the University of Tokyo, Japan, in 1998 and 2000, respectively. In
2000, she joined NTT Communication Science
Laboratories, Kyoto. Her research interests include array signal processing, blind source separation applied to speech signals, and auditory
scene analysis. She received the TELECOM
System Technology Award from the Telecommunications Advancement Foundation in 2004,
the Best Paper Award of the IWAENC in 2003 and the 19th Awaya Prize
from Acoustical Society of Japan (ASJ) in 2001. She is a member of theIEEE and the ASJ.