2007 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

October 21-24, 2007, New Paltz, NY

A TWO-STAGE FREQUENCY-DOMAIN BLIND SOURCE SEPARATION METHOD FOR
UNDERDETERMINED CONVOLUTIVE MIXTURES
Hiroshi Sawada, Shoko Araki and Shoji Makino
NTT Communication Science Laboratories
NTT Corporation
2-4 Hikari-dai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japan{sawada, shoko, maki}@cslab.kecl.ntt.co.jp
ABSTRACTSTFTThis paper proposes a two-stage method for the blind separationof convolutively mixed sources. We employ time-frequency masking, which can be applied even to an underdetermined case where
the number of sensors is insufficient for the number of sources.
In the first stage of the method, frequency bin-wise mixtures areclassified based on Gaussian mixture model fitting. In the second
stage, the permutation ambiguities of the bin-wise classified signals are aligned by clustering the posterior probability sequences
calculated in the first stage. Experimental results for separating
four speeches with three microphones under reverberant conditions show the superiority of the proposed method over existing
methods based on time-difference-of-arrival estimations or signal
envelope clustering.

Bin-wiseclassification

mixtures are classified by the procedure explained in Sec. 3. Since
the classification is not based on TDOA estimations, but on a frequency bin-wise instantaneous model, it is robust as regards room
reverberations and microphone characteristics as long as the frame
size in the frequency analysis is long enough to cover the main partof the impulse responses. In the second stage, the bin-wise classified signals are grouped together for the same original source.
This task is almost identical to that of the permutation problem[8, 9] which occurs in frequency-domain ICA-based BSS for a determined case (M ≥ N ). We adapt our recent idea [10], which
we developed for a determined case, to an underdetermined case.
Section 4 describes the procedure in detail.

The technique for estimating individual source components fromtheir mixtures at multiple sensors is known as blind source separation (BSS) [1]. With acoustic applications of BSS, such as solving
a cocktail party problem, signals are mixed in a convolutive manner. Let s1 , . . . , sN be source signals and x1 , . . . , xM be sensor(microphone) observations. The convolutive mixture model is formulated asxj (t) =l hjk (l) sk (t − l),j = 1, . . . , M,

2. SYSTEM OVERVIEW(1)

where t represents time and hjk (l) represents the impulse response
from source k to microphone j.
If we have a sufficient number of microphones (determinedcase), i.e., M ≥ N , linear filters that are estimated for example byindependent component analysis (ICA) [2] effectively separate themixtures. However, if the number of microphones is insufficient(underdetermined case), i.e., M < N , such linear filters do notwork well. Instead, time-frequency (T-F) masking [3, 4, 5] or l1 norm minimization [6, 7] is widely used to separate such underdetermined mixtures. Most of the existing techniques for underdetermined convolutive BSS rely on time-difference-of-arrival (TDOA)
estimations for each source at the microphones. They work effectively under low reverberant conditions. However, under severe
reverberant conditions, TDOA estimations become unreliable andsuch techniques do not work well. Also, any microphone directivity complicates the TDOA estimation.
This paper proposes a new method for underdetermined convolutive BSS. The method follows the frequency-domain approach,and consists of two stages. In the first stage, frequency bin-wise

978-1-4244-1619-6/07/$25.00 ©2007 IEEE

Permutationalignment

Figure 1: System flow

1. INTRODUCTIONN
k=1

InverseSTFT

T-F masking

139

Let us begin with a system overview of the proposed method.
Figure 1 shows the system flow, and Fig. 2 shows some examples of spectrograms and classification results for a 3-source 2microphone case.
2.1. Short-time Fourier transform (STFT)First, the sensor observations (1) are converted into frequencydomain time-series signals xj (f, t) by an STFT. The time index
t is then down-sampled with a distance equal to the STFT frame
shift. If the STFT frame size is long enough to cover the main partof the impulse responses hjk (l), the time-domain convolutions in(1) can be well approximated with frequency-domain multiplications:xj (f, t) =N
k=1 hjk (f )sk (f, t) ,j = 1, . . . , M.

(2)

Let us rewrite (2) in a vector notation:N
k=1 hk (f )sk (f, t)x(f, t) =T(3)where hk = [h1k , . . . , hM k ] is the vector of frequency responses
from source sk to all sensors, and x = [x1 , . . . , xM ]T is called an
observation vector.

2007 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics(a) Sources(b) Mixtures

October 21-24, 2007, New Paltz, NY(c) Bin-wise classification(d) Permutation aligned
classification(e) Separated signals

Figure 2: Example spectrograms for sources (a) mixtures (b) and separated signals (e), and classification results (c) and (d).

source in an observation vector x(f, t). The procedure consists of
two stages, which are also described in this section and the next.
The first stage classifies observation vectors x(f, t) in a frequency bin-wise manner. Thus, here we omit the frequency dependence in (3) and (4) for simplicity:

2.2. Time-frequency (T-F) masking and inverse STFT
To separate the mixtures, we employ T-F masking, which is effective for separating sparse sources. A sparse source can be characterized by the fact that the source amplitude is close to zero mostof the time. A frequency-domain speech source is a good example.
The mixing model (3) can be approximated to a simpler form
x(f, t) ≈ hk (f )sk (f, t)

x(t) =(4)yk (f, t) =x (f, t) if P (C |x) ≥ P (C |x) , ∀k = kJkp(x|ai , σi ) =

if P (Ck |x) < P (Ck |x) , ∃k = k
(5)
for k = 1, . . . , N where J is the index of an arbitrarily selectedreference microphone that is used to construct the separated signals.
At the end of the flow, time-domain output signals yk (t) arecalculated with an inverse STFT to the frequency-domain separated signals yk (f, t).
0



2
1
||x − (aHi x) · ai ||
exp −
2 M
2
(πσi )σi(7)where ai is the mean vector with unit-norm ||ai || = 1 and σi2 isthe variance. Since (aHi x) · ai is the orthogonal projection of x
onto the subspace spanned by ai , ||x − (aHi x) · ai || represents theminimum distance between the point x and the subspace, which
implies the probability that x belongs to the i-th class.
Now the density function for data x can be given byp(x|θ) =N
i=1 αi p(x|ai , σi )(8)with the parameter set θ = {a1 , σ1 , α1 , . . . , aN , σN , αN }. The
mixture ratios αi should satisfy 0 < αi < 1 and α1 + · · · + αN =1. Since there is some computational difficulty involved in maximizing the log-likelihood Tt log p(x(t)|θ) of T samples directly,

3. BIN-WISE CLASSIFICATIONThis section and the next explain how to calculate the posterior
probability P (Ck |x) that the k-th source is the most dominant

978-1-4244-1619-6/07/$25.00 ©2007 IEEE

(6)

We changed the source subscript from k to i. This is because there
is permutation ambiguity in the frequency bin-wise classifications.
The effect of permutations is demonstrated in Fig. 2 (c). Suchpermutation ambiguities will be aligned in the second stage, which
is detailed in the next section.
We see in (6) that classification can be performed according to
the information on the vectors h1 , . . . , hN . To remove the effectof source amplitude from x, let us normalize them to have unit
norm
x(t)
hi
si (t)
x(t) ←=·.
||x(t)||||hi || |si (t)|Unknown phase si (t)/|si (t)| ambiguity still remains in x(t). To
model such vectors for each source, we follow the idea describedin [6] and employ the following complex Gaussian density function

for sparse sources sk . The summation in (3) is eliminated because
the sources are assumed to be close to zero most of the time. The
subscript k is the index of the most dominant source for each timefrequency slot (f, t). The index k should be identified or estimatedfor each (f, t) to separate the sources by T-F masking.
Let us classify the observation vectors x(f, t) into N classes
c1 , . . . , cN , each of which corresponds to a source signal sk . A
vector x(f, t) should belong to class Ck if the source sk is the mostdominant in the observation x(f, t). Through the processes of“Bin-wise classification” and “Permutation alignment”, which
are depicted in Fig. 1 and explained in the next two sections, we
calculate the posterior probabilities P (Ck |x) for all source indicesk and for all observation vectors x(f, t). Then, T-F masking forsource separation is performed bykN
i=1 hi si (t) ≈ hi si (t) .

140

2007 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

we typically employ an EM algorithm to estimate the parameters
θ. By introducing the posterior probability
αi p(x(t)|ai , σi ),p(x(t)|θ)

f = 766 Hz
Envelope

P (Ci |x(t), θ) =(9)N
i=1 P (Ci |x(t), θ) log αi p(x(t)|ai , σi )Envelope

T
t

(10)

to be maximized in the EM algorithm.
In the E-step, P (Ci |x(t), θ) is calculated by using (9) and (8)based on the current parameters θ. In the M-step, we maximize the
Q function by adjusting the parameters in θ. Considering the factthat both x and ai have a unit-norm, the optimal ai is given by the
eigenvector corresponding to the maximum eigenvalue ofTHt P (Ci |x(t), θ) · x(t)x (t) .

By performing the operations described in the previous section independently for all frequencies f , we have the posterior probability P (Ci |x(f, t)) according to (9) for i = 1, . . . , N and all timefrequency slots (f, t). However, since the class order c1 , . . . , cN
may be different from one frequency to another (Fig. 2 (c)), we
need to reorder the indices such that the same index corresponds
to the same source over all frequencies (Fig. 2 (d)). In other words,
we determine a permutation Πf : {1, . . . , N } → {1, . . . , N } for
all frequencies f , and then update the posterior probabilities by

2

3Time (sec)

0

1

2

4

5

3

−0.07 −0.07−0.28 −0.28

978-1-4244-1619-6/07/$25.00 ©2007 IEEE

4

5

4

5

g = 906 Hz
1

1

2

3
Time (sec)f ∈F k=1

ρ(vif , ck ) i=Π (k) ,f(11)

where F represents the set of all frequency bins. The centroid ck
is calculated for each source as the average of activity sequenceswith the current permutation Πf :ck (t) ←1v f (t) i=Π (k) , ∀ k, t ,f|F| f ∈F i

(12)where |F| is the number of elements in the set F. The permutations Πf are optimized to maximize the correlation coefficients ρ
between activity sequences and the current centroid:NΠf ← argmaxΠρ(vif , ck ) i=Π(k) , ∀ f ∈ F .

(13)

k=1

These two operations (12) and (13) are iterated until convergence.
Then we employ a local fine optimization to improve the separation further. It maximizes the sum of the correlation coefficientsover a set of selected frequencies R(f ) for a frequency f :0.44

 ρ(v , v ) ρ(v , v ) ρ(v , v )  0.51 −0.20 −0.29
 ρ(v , v ) ρ(v , v ) ρ(v , v ) = −0.25 0.46 −0.23 .
f
g
1
3
f
g
2
3
ρ(v3f , v3g )1Nand those calculated with the posterior probability sequences are
f
g
1
2
f
g
2
2
ρ(v3f , v2g )0

J ({ck }, {Πf }) = ρ(v , v ) ρ(v , v ) ρ(v , v )  0.01 0.07 −0.10
 ρ(v , v ) ρ(v , v ) ρ(v , v ) = −0.09 0.10 −0.15 ,
f
g
1
1
f
g
2
1
ρ(v3f , v1g )

5

2

f

Such a permutation problem has been extensively studied for
frequency-domain ICA-based BSS applied to a determined case(M ≥ N ). The proposed method follows a major approach basedon the correlation coefficients ρ(vif , vjg ) of two sequences vif (t),vjg (t) that express the separated signal activities [8, 9, 10]. The
correlation coefficient of any two sequences is bounded by −1 ≤ρ(vif , vjg ) ≤ 1 , and becomes 1 if the two sequences are identical.
Although signal envelopes vif (t) ← |yi (f, t)| are widely used
[8, 9] to represent the separated signal activities, here we use theposterior probabilities vif (t) ← P (Ci |x(f, t)) instead. Figures 3
and 4 show examples of separated signal envelopes and posterior
probability sequences, respectively. The correlation coefficients
calculated with the signal envelopes are
f
g
1
3
f
g
2
3
f
ρ(v3 , v3g )4

The posterior probability sequences resulted in higher correlation
coefficients for the same source. This is because the posterior
probabilities represent active signals with values uniformly close
to 1, whereas envelopes represent active signals with various values. In this sense, a posterior probability has similar characteristics
to a dominance measure developed in [10] for determined BSS.
To determine permutations Πf , we cluster such activity sequences vif (t) for each source. For computational efficiency, we
first employ a global rough optimization that maximizes

P (Ck |x) ← P (Ci |x) i=Π (k) , k = 1, . . . , N .

f
g
1
2
f
g
2
2
f
ρ(v3 , v2g )3

Figure 4: Posterior probability sequences for three classes

4. PERMUTATION ALIGNMENT

f
g
1
1
f
g
2
1
f
ρ(v3 , v1g )

2

4

0

T
t P (Ci |x(t), θ).

1

f = 766 HzTH
2t P (Ci |x(t), θ) · ||x(t) − (ai x(t)) · ai ||T
M · t P (Ci |x(t), θ)and αi = (1/T )0

1Other parameters are optimized by
σi2 =

2

Figure 3: Envelopes of three separated signals

Posterior probabilityR=4

g = 906 Hz

we construct the so-called Q function
Q(θ) =

October 21-24, 2007, New Paltz, NYNΠf ← argmaxΠ

0.55

g∈R(f ) k=1

141



ρ(vif , vig ) i=Π(k),i =Π (k) . (14)g2007 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics

380 and 450 ms). We compared four strategies for permutationalignment. “TDOA” corresponds to the existing techniques mentioned in the Introduction. “Envelope” and “Posterior” use signal envelopes |yi | and posterior probability sequences P (Ci |x),respectively, for the permutation alignment procedure described inSec. 4. “Optimal” corresponds to optimal permutations calculatedwith knowledge of the original source signals. This is not a practical solution but used to determine the upper bound of the results.
We observed the following tendencies. “TDOA” performed
moderately well in short distance (60 cm) or low reverberant (130ms) cases, but did not perform well under reverberant conditions
with a distance of 120 cm. “Envelope” did not perform so well inmany cases. The proposed “Posterior” performed the best among
the practical solutions (except “Optimal”) in all cases.
We also made recordings in a room using a portable audio
recorder with two microphones, and separated the mixtures of threespeeches. Sound examples can be found on our web site [11].

Room size:

Loudspeakers

4.45 × 3.55 × 2.5 m

150°

245°

Distance:
60, 120cm

Reverberation time:
RT60 = 130 ~ 450 ms

Sampling rate:
8 kHz

Microphones

On edges of 4cm triangle

STFT frame size:

70°1024 points (128 ms)

315°

STFT frame shift:
256 points (32 ms)

Height of microphones and loudspeakers: 120 cm

Figure 5: Experimental setup

14
12
10
8
6

Source to microphones: 120 cm
16
SIR improvement (dB)

SIR improvement (dB)

Source to microphones: 60 cm
16

OptimalPosterior
Envelope
TDOA

14
12

6. CONCLUSIONS

10This paper presented a method for underdetermined convolutive
BSS. The method has a two-stage hierarchical structure that classifies observation vectors x(f, t) into N sources. Frequency binwise samples along the time axis are classified in the first stage,
and then the posterior probability sequences are clustered along
the frequency axis in the second stage. The experimental results
show the effectiveness and robustness of the proposed method evenunder reverberant conditions.

8
6

200
300
400
Reverberation time (ms)

200
300
400
Reverberation time (ms)Figure 6: Experimental results

The set R(f ) preferably consists of frequencies g where a highcorrelation coefficient ρ(vif , vig ) would be attained for vif and vig
corresponding to the same source. We typically select adjacentfrequencies and harmonic frequencies [9]. The fine local optimization (14) is performed for one selected frequency f at a time, andrepeated until no improvement is found for any frequency f .

7. REFERENCES[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] O. Yilmaz and S. Rickard, “Blind separation of speech mixtures via
time-frequency masking,” IEEE Trans. Signal Processing, vol. 52,
no. 7, pp. 1830–1847, July 2004.
[4] M. Mandel, D. Ellis, and T. Jebara, “An EM algorithm for localizing
multiple sound sources in reverberant environments,” in Advances in
Neural Information Processing Systems 19, B. Schölkopf, J. Platt,and T. Hoffman, Eds. Cambridge, MA: MIT Press, 2007.
[5] S. Araki, H. Sawada, R. Mukai, and S. Makino, “Underdeterminedblind sparse source separation for arbitrarily arranged multiple sensors,” Signal Process., vol. 87, no. 8, pp. 1833–1847, 2007.
[6] P. D. O’Grady and B. A. Pearlmutter, “Soft-LOST: EM on a mixtureof oriented lines,” in Proc. ICA 2004 (LNCS 3195). Springer, Sept.
2004, pp. 430–436.
[7] S. Winter, W. Kellermann, H. Sawada, and S. Makino, “MAP based
underdetermined blind source separation of convolutive mixtures byhierarchical clustering and L1-norm minimization,” EURASIP Journal on Advances in Signal Processing, 2007, Article ID 24717.
[8] 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.
[9] 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 Processing, vol. 12, no. 5, pp. 530–538, Sept. 2004.
[10] H. Sawada, S. Araki, and S. Makino, “Measuring dependence ofbin-wise separated signals for permutation alignment in frequencydomain BSS,” in Proc. ISCAS 2007, 2007, pp. 3247–3250.
[11] [Online]. http://www.kecl.ntt.co.jp/icl/signal/sawada/demo/ubssconv/

5. EXPERIMENTAL RESULTS
We performed experiments to separate four speech sources withthree microphones. We measured impulse responses hjk (l) under the conditions shown in Fig. 5. Mixtures at the microphones
were made by convolving the impulse responses and 6-second English speeches. The computation time was around 8 seconds for
6-second speech mixtures. The program was coded in Matlab
and run on an AMD 2.4 GHz Athlon 64 processor. The separation performance was evaluated in terms of signal-to-interference
ratio (SIR) improvement. The improvement was calculated by
OutputSIRi − InputSIRi for each output i, and we took the average over all outputs. These two types of SIRs are defined as
InputSIRi = 10 log10t|2
l hJ i (l)si (t − l)|t|OutputSIRi = 10 log10

October 21-24, 2007, New Paltz, NY

k=i

t|

2
l hJ k (l)sk (t − l)|2t |yii (t)|2k=i yik (t)|

(dB),

(dB),

where J ∈ {1, . . . , M } is the index of one selected reference sensor, and yik (t) is the component of sk that appears at output yi (t),i.e. yi (t) = N
k=1 yik (t).
Figure 6 shows the average SIR improvements for 8 combinations of 4 speeches. We had two source-to-microphones distances(60 and 120 cm) and six reverberation times (130, 200, 270, 320,978-1-4244-1619-6/07/$25.00 ©2007 IEEE

142