Hierarchical clustering applied to overcomplete BSS for convolutive mixtures
Stefan Winter , Hiroshi Sawada, Shoko Araki, Shoji Makino
NTT Communication Science Laboratories, NTT Corporation
2-4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japanwifan,sawada,shoko,maki @cslab.kecl.ntt.co.jp

Abstract
In this paper we address the problem of overcomplete BSS
for convolutive mixtures following a two-step approach. In
the first step the mixing matrix is estimated, which is thenused to separate the signals in the second step. For estimating
the mixing matrix we propose an algorithm based on hierarchical clustering, assuming that the source signals are sufficiently sparse. It has the advantage of working directly on
the complex valued sample data in the frequency-domain. It
also shows better convergence than algorithms based on selforganizing maps. The results are improved by reducing thevariance of direction of arrival. Experiments show accurate
estimations of the mixing matrix and very low musical tonenoise even in reverberant environment.

1. Introduction
High quality separation of speech sources is an important
prerequisite for further processing like speech recognitionin environments with several active speakers. Often the underlying mixing process is unknown, which requires blindsource separation (BSS). In general we can distinguish twocases depending on the number of sources and the numberof sensors : 

overcomplete BSS(under-) complete BSS 

Since undercomplete BSS () can be reduced to complete BSS () [1] we refer to both by complete BSS.
Most approaches assume complete mixtures [2, 3], but in reality often the contrary is true. While the area of overcomplete BSS has obtained more and more attention [4–12], it
still remains a challenging task.
Several of the proposed algorithms are based on histograms and developed for only two sensors [4–6]. Some
could, in principle, be enhanced for higher dimensions .
But since histograms are based on densities, the so-called
curse of dimensionality [13] sets practical limits on the number of usable sensors. Another problem occurs with complex numbers, which cannot be handled straightforwardly
by histograms, but are necessary if BSS is performed in thefrequency-domain. Some methods approach complex numbers by applying real-valued algorithms to the real and imag-

inary part or amplitude and phase [7, 8], which is not always applicable. Some approaches extract features like thedirection-of-arrival (DOA) or work on the amplitude relation
between two sensor outputs [4, 5, 9, 10]. In both cases only
two sensors can contribute, no matter how many sensors areavailable.
Other algorithms like GeoICA [12] or AICA [11] resemble self-organizing maps (SOM) and could more easily beapplied to convolutive mixtures. However, their convergence
depends on initial values.
In this paper we propose the use of hierarchical clustering
embedded into a two-stage framework of overcomplete BSS
to deal with convolutive mixtures in the frequency-domain.
This method can work directly on the complex valued samples. While it does not limit the usable numbers of sensors,it also prevents the convergence problems which can occur
with SOM based algorithms.
After estimating the mixing matrix in the first stage, amaximum a-posteriori (MAP) approach is applied to finallyseparate the mixtures, assuming statistical independence andLaplacian pdfs for the sources [14].
In Sec. 2 we first explain the general framework before
we give details about the hierarchical clustering in Sec. 3 andthe MAP based source separation in Sec. 4. After this, we
present experimental results in Sec. 5 demonstrating the performance for convolutively mixed speech data in a real roomwith reverberation timems.

The author is on leave from the Chair of Multimedia Communications
and Signal Processing, University Erlangen-Nuremberg.

2. General framework !#"$%'&(&)& * .0/1 2+
,


-
*

43567&)&)&86
9 :> = @:. /  ' 4;2< ;?< / ,A,  CBDA,
(1)@
where /#  denotes the impulse response from source " tosensor 3 .

We will consider a convolutive mixing model with speechsourcesandsensors that
. The mixing
yield linearly mixed signals
can be described byInstead of solving the problem in the time-domain, wechoose a narrowband approach in the frequency domain by
applying a short-time discrete Fourier transform (STDFT).
Thus time-domain signals
andare converted into frequencydomain time-series

O # PQG . < # RI(&)&)&SI .UT# VKWE1M # FHG  <  JI)&)&(&I8 9  LKNM
X'ZYUI [0\ G ]?<^,YUI [0JI)&(&)&SI_] 9 ZYUI [0LKNM

Workshop on Statistcal and Perceptual Audio Processing SAPA-2004, 3 Oct 2004, Jeju, Korea       
S
T
D
F
T

 
I S
n
v T
e D
r
s F
e T

BSR


PSfrag replacements

BMMR

Permutation

DOA
improvement

DOA
estimation

Figure 1: Overall unmixing system, YUI[0 G < ,YUI[0RI(&)&)&?I T ,YUY I [0LKNM
 I Y I?&(&)&JIY  !B   Y , Y  T 9 [C8,Y  "$%I)&(&)& I8*

and 
by an  -point STDFT, respectively.
Thereby 
 (  : sampling frequency; : time"!
 
dependence). Let us define 
as a matrix
whose elements are the transformed impulse responses. We
mixing veccall the column vectors #tors and can approximate the mixing process as*ZYUI [07 ,Y 8X'ZYUI [0&


$

(2)

3. Blind Mixing Model Recovery
Several algorithms have been proposed so far for BMMR.
They usually have in common that they assume a certain degree of sparseness of the original signals. This assumption
leads to the conclusion that the samples in the mixed vector
space cluster around the true mixing vectors #.
This becomes clear when we consider the most sparse casewhen only a single source is active. Let us rewrite (2) as*,YUI[0 ,Y 

9
:,YUI [0 4;?<  ZY  ]  ,YUI[0*&(3)Assuming only one active source at ,YUI[0 means that the vector pointing to the resulting mixed sample *,YUI[0 is a scaled
version of the corresponding mixing vector 8,Y  . However,YThis reduces the problem from convolutive to instantaneous mixtures in each frequency bin . For simplicity wewill omit the dependence on frequency and time. Switching
to the frequency domain has the additional advantage thatthe sparseness of the sources is increased [7]. This is very
important, since the hierarchical clustering is based on theassumption of sparse sources, which is approximately true
for speech signals.
The disadvantage of narrowband BSS in the frequency
domain is the permutation problem, which results in wrongalignments of the frequency bins. In our framework we use
a DOA based method to reduce the permutation problem [3].
We also apply the minimum-distortion-principle [2] to solve
the scaling problem.
In complete BSS the mixing matrix  is square and (assuming full rank) invertible. Therefore the BSS problem canbe solved by either inverting an estimate of the mixing matrix
or directly estimating its inverse and solving (2) for .

X

However, this approach does not work in overcomplete
BSS where the mixing matrix is not invertible. Therefore
we follow a two-stage approach as utilized in [7, 12] consisting of blind mixing model recovery (BMMR) and blindsource recovery (BSR). To estimate the mixing matrix in the
BMMR step, we propose the use of hierarchical clustering
as described in detail in Sec. 3. To eventually separate thesignals in the BSR step, we utilize a MAP based approach.
Finally the inverse STDFT is applied to obtain time-domain
signals. The overall system is depicted in Fig. 1.

#



#

depending on the actual sparseness of the source signals andalso the mixing system, the mixed signals will also have components of other mixing vectors. Therefore the mixtures willbe spread around the mixing vectors but still form clusters.

3.1. Hierarchical Clustering
To avoid the problems discussed in Sec. 1, such as the curseof dimensionality or poor convergence, we propose the useof a hierarchical clustering algorithm for finding the clusters around the mixing vectors. We follow an agglomerative(bottom-up) strategy. [13]. This means that the starting pointis the single samples, considering them as clusters whichcontain only one object. From there clusters are combined,
so that the number of clusters decreases while the average
number of objects per cluster increases. In the following weassume phase and amplitude normalized samples
&%

'




' (*),+.-0/21(4)where 354 1 denotes the phase of the first component of  .
The combination of clusters into new clusters is an iterative process and based on the distance between the current
clusters. Starting from the normalized samples, the distanceLevel,+.-/-1032 Hz

0.3
10
9
8
c
7
6
5
4
3
2
1

0.2
0.1

*

imag

' ()

PSfrag replacements
1

3

2

5

4

6

7

8

9

10

11 Object

0
-0.1
-0.2

Estimated hi
Original hi

-0.3
0

0.1

0.2

0.3

0.4

Figure 2: Linking closest clusters

real 

0.5

&



0.6

0.7

0.8

Figure 4: Estimation of mixing vectors   <replacements

(of the original signals. Special attention must be paid to theremaining samples (outliers), which are randomly scatteredin the space between the mixing vectors. Usually they are faraway from other samples and will be combined with otherclusters only at higher levels of the clustering process (i.e.
when only few clusters are left). This led us to the idea to set
the final number of clusters at a high number   

54



Figure 3: Illustration of distances
between each pair of clusters is calculated, resulting in a distance matrix. At each level of the iteration the two clusters
with the least distance are combined and form a new binary
cluster (Fig. 2). This process is called linking and repeateduntil the amount of clusters
  has decreased to a predetermined
number  ,( : total number of samples).

For measuring the distance between clusters, we have to
distinguish between two different problems. First we needa distance measure    1   that is applicable todimensional complex vector spaces. While there are several
possibilities, we currently use the Euclidean distance basedon the normalized samples, which is defined by

  

 I 



U I  '  B  JI( B    (5)where 4I   stands for the inner product and  for complex
conjugation.
 1  

  1



 1When a new cluster is formed, we need to enhance thisdistance measure to relate the new cluster to the other clusters. The method we employ here
( is called the nearestneighbor technique. Let  and  denote two ( clusters asillustrated in Fig. 3. Then the distance  
betweenthese clusters is defined as the minimum distance between
its samples by< <I 

U < I '   "! # %$    ! #  U I   &
 



(

1

1

C

  1  

(6)

As mentioned earlier, most of the samples will clusteraround the mixing vectors # , depending on the sparseness

 &

(7)

By doing so, we avoid linking these outliers with the clusters
around the mixing vectors # and therefore distortions. This
results in more robustness. More important, however, is thefact that we avoid combining desired clusters. Since the outliers are often far away from other clusters it might happen
that desired clusters are closer to each other than to outliers.
Due to the distance( calculation the computational cost of
this approach is 6 87[15], where 7 denotes the numberof samples per frequency bin. An example for the resulting
clusters is shown in Fig. 4. Experimental details are given inSec. 5.

 

3.2. Estimation of mixing matrix

Assuming that the clusters around the mixing vectors # havethe highest densities and therefore the highest number ofsamples we finally choose thelargest clusters. Thereby
the number of sources must be known. To obtain the mixing vectors, we average over all samples of each cluster  : I
! #:9

$ "' (8)where  denotes the cardinality of cluster  . Thereby weassume that the influence of other sources has zero mean.
#

'



'

'



'





3.3. Advantages of hierarchical clustering
Among the most important advantages of the described hierarchical clustering algorithm is the fact that it works directly on the sample data in any vector space of arbitrary
dimensions. The only requirement is the definition of a distance measure for the considered vector space. Therefore, itcan easily be applied to complex valued data that occurs infrequency-domain convolutive BSS.
No initial values for the mixing vectors # are required.
This means, in particular, that if the assumption of clusters
with high densities around the mixing vectors is true, thenthe algorithm converges to those clusters.
Besides choosing a distance measure, there is only thesingle parameter  that determines the number of clusters.
Experiments have shown that the choice for this parameteris quite insensitive as long as it is above a certain limit thatwould combine desired clusters. Its choice is, in general, related to the sparseness of the sources. The sparser the signals
are, the smaller the value of  can be chosen, because thenumber of outliers that must be avoided will be smaller.
While the considered signals must have some degree ofsparseness, they do not have to be statistically independent at
this point to obtain clusters.

180



C

Experiments have shown that as long as there are clusters
around the mixing vectors # , the estimation results are ofhigh quality. Even if the assumption of clear clusters is nottrue for all mixing vectors, the remaining ones are not influenced by poor estimation of others. In order to improve thepoorly estimated mixing vectors, we can utilize DOA information.
While the mixing matrix is different for each frequency
bin, the phase difference 3

between two
arbitrary components of a mixing vector #contains information about the relative physical position of
its corresponding source . Assuming a linear sensor array ina far-field situation with plain wave fronts and no reverberation, the DOA is given by   @ < C ?@ /V G @ <  &)&(& @ T #KWM   + <  Y3C  <Z/

(9)

<V/



where denotes the sound velocity. stands for the distance between the sensor that corresponds to the first component
and the sensor that corresponds to the -th component.
Since is theoretically constant for all frequency bins,
we can consider the DOA of the -th signal( as a random variable (RV) with mean and variance . While even the
DOA of the original mixing matrix has a variance larger than(Fig. 5), the results for the estimated mixing matrix can be
improved if the variance of its DOAs is reduced.
For this purpose we define a new RVwith reduced
variance by@@ /< 








"

3



 

 ?      B ^I 5(  5  & (10)While its mean is still   , its variance    can be adjusted by and yields((
0   !  &
(11)

2

140

DOA

120
100
80
60
40
20
0
0

50

100

150

@ <

We apply the new DOA by adjusting the phase of the mixstays fixed, because
ing vectors # . The first component

200

250
300
frequency bin

350

400

450

500

Figure 5: DOAs of original mixing matrix
we want to change only the relative phase. The remaining
components are updated by@ / @ / #"%$ 9'& $)(* 9'+'+  3  (12)@@where '< ',0 < Z denotes the phase of <  and   V isgiven by Y  <V/  )      7&(13)'

3.4. Reduction of DOA variance"

1st source
2nd source
3rd source

160

'

- 1

)

-

3

3



3

Since we do not need absolute DOA information, this improvement fully complies with the blind approach of BSS.

4. Blind Source Recovery
The unmixed signals cannot be directly obtained, because themixing matrix cannot be inverted in overcomplete BSS. Several approaches have been proposed to solve blind source recovery [14]. Among those we chose the shortest-path algorithm which is based on maximum a-posteriori (MAP) estimation, assuming statistical independence and Laplacian
pdfs for the sources. Given the mixed signals  and the mixing matrix  , the sources are recovered byX

9:  .; -0 / 4;?< ]  &

X

'

'

(14)



This equation can be interpreted as finding the shortest-path
decomposition, based on the mixing vectors # for each sample   separately. It means that each sample is assigned toexactlysignals. While (14) can, in general, be solved forreal numbers by linear programming, we explicitly compute 

decompositions and choose the one that  possible
9
minimizes 1  ;?< ]  . Taking a selection of  mixing vectors 2 &)&(& 23 2 , the decomposition is calculated by
X 6G 2 &(&)& 2) 2!K < O " < I)&)&(&I " T 5411I)&(&)& I7 (15)6
Since each mixing vector has  components and assuming
full rank, the inverse exists.

all

#

'

'

#

1

# 1

#

+



^1 1   
(1& ,+.- - 032 Hz

 Y

Y

imag

0.3
Table 1: Experimental conditions
Direction of sources ,   , 
0.2Distance of sensors
 mmLength of source signals
  seconds
0.1*
Reverberation timems'( ) 0 kHz
Sampling frequency Window typevon Hann-0.1
 points
Filter length-0.2
 points
Shifting interval
PSfrag replacements
Cluster threshold  (const ) 100-0.3
Variance factor
0.8

Estimated hi
Original hi0 Source  Average
 Source
13.9
11.7
13.50

Table 2: Detailed performance
Scenario 1
SIR (dB)SDR (dB)

Source
14.8
13.39

6.83

10.55


Scenario 2 Source SourceSIR (dB)SDR (dB)

10.5
7.47

6.4
2.82

Source9.3
5.99

SIR (dB)SDR (dB)

11.1
9.65

9.9
4.05

Source8.9
5.95


Scenario 3 Source Source5.43

 Average
9.95
6.55

We performed experiments with the proposed algorithm usspeech signals andsensors. The signals
ing
were taken from the Acoustical Society of Japan (ASJ) continuous speech corpus. The convolution was done with room
impulse responses that were recorded at our laboratory. Further experimental conditions are given in Table 1.
As performance measure, we used the signal-tointerference ratio (SIR) and the signal-to-distortion ratio(SDR) based on time-domain signals. For this purpose, let
us denote the output signals in the time-domain by. If the separation was ideal, they would equal filtered versions of the input signals. We calculated the
SIR at output by  #"C'&)&(& *   
1  #  SIR S6 (16)
1     
where    is the portion of 1 #  that comes from source
signal  #  and   #  1   B    denotes the interference from other sources.
To calculate the SDR, #  is further decomposed into
 #       >= # @  .  The reference signal    isgiven by    # ' 1 ;?< / #AZ   #B5A, , where 3 denotes thesensor that is closest to the location of source " . The scalar is chosen to minimize the power of #  . Then the SDR at
"

(





(







)

)

0.3

0.4

real 

0.5

& 

0.6

0.7

0.8

Figure 6: Effect of too few clusters"

output is obtained by Average
8.73





0.2

10.26

5. Experimental Results

Q 

0.1

SDR

(

2(! #"$% 1& 1*  '#   )( +-. ,
)

(

(17)

As an upper limit for the performance of the whole system, scenario 1 in Table 2 shows the separation results when
the original mixing matrix is used. This means that the permutation problem does not occur and the BSR part is given
the best possible input.
Scenario 2 gives the results if we use the estimated mixing matrix without reduction of DOA variance. Scenario 3shows the results if the estimated mixing matrix is used together with reduction of DOA variance.
Figure 4 gives an example of the clustering forwith a total of 71 samples in each/0 Hz and 
frequency bin. To visualize the normalized samples, the realpart of the first mixture
versus the imaginary part of the(second mixture is plotted. The largest clusters (black)around the normalized original mixing vectors # (dashed)
can be clearly seen and result in precise estimations (solid).
In contrast, Fig. 6 shows the clustering which resulted
from the same conditions as Fig. 4 but a reduced number ofclusters   . Reducing the number of final clusters led
to the combination of two desired clusters and eventually towrong estimates for two mixing vectors. However, this figure also shows, that the estimation of the remaining mixing
vector is not affected by the estimation of other vectors.
The influence of different values of the variance factor
is depicted in Figs. 7 and 8. The solid lines in Fig. 7 givethe average SIR if the mixing matrix is estimated. Each line
shows the result for three mixed sources that were chosenfrom four different speech signals. In the same way, thedashed lines give the average SDR. The corresponding lines
are marked by the same symbol. For both SIR and SDR avalue ofand, therefore, reduced DOA variance yields
improved results.
For Fig. 8 the original mixing matrix was used instead ofthe estimated mixing matrix. The depicted graphs are averaged over all different mixture combinations. As expected,

6(1

1

P<

Y



F 



P



 

11

[3] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Arobust and precise method for solving the permutation
problem of frequency-domain blind source separation,”in Proc. ICA 2003, 2003, pp. 505–510.

SIR
SDR

10
9
8

SIR, SDR (dB)76[4] O. Yilmaz and S. Rickard, “Blind separation of
speech mixtures via time-frequency masking,” IEEE
Transactions on Signal Processing, Jul 2004. [Online].
Available: http://eleceng.ucd.ie/ srickard/bss.html

5
4
3
2
1

replacements

0

0

0.1

0.2

0.3

0.4

0.5
0.6
Variance factor e

0.7

0.8

0.9

1

Figure 7: Performance comparison with different mixtures
14
SDR
SIR

SIR, SDR (dB)

12

8

4
0

[6] L. Vielva, I. Santamaria, C. Pantaleon, J. Ibanez, and
D. Erdogmus, “Estimation of the mixing matrix forunderdetermined blind source separation using spectral estimation techniques,” in Proc. EUSIPCO 2002,
vol. 1, Sep 2002, pp. 557–560.
[7] P. Bofill and M. Zibulevsky, “Blind separation of more
sources than mixtures using sparsity of their short-time
fourier transform,” in Proc. ICA 2000, Jun 2000, pp.
87–92.

10

6

replacements[5] S. Rickard and O. Yilmaz, “On the approximate Wdisjoint orthogonality of speech,” in Proc. ICASSP
2002, 2002, pp. 529–532.

0.1

0.2

0.3

0.4
0.5
0.6
Variance factor εFigure 8: Influence of variance factor
matrix

0.7

0.8

0.9

1

 on original mixing

 -, where the DOA
the best performance is obtained forand therefore the mixing matrix remains unchanged.
Subjective evaluation of the separated sources showedvery low musical tone noise.

6. Conclusion
We proposed the application of hierarchical clustering embedded into a two-stage framework of overcomplete BSS
for convolutive speech mixtures. This method can work
directly on the complex mixture samples. It also prevents
the convergence problems which can occur with SOM based
methods like GeoICA. Experimental results confirmed thatthe assumption of sparseness and, therefore, clusters aroundthe mixing vectors is sufficiently fulfilled for convolutively
mixed speech signals in the frequency domain.

7. References[1] S. Winter, H. Sawada, and S. Makino, “Geometricalinterpretation of the PCA subspace method for overdetermined blind source separation,” in Proc. ICA 2003,
2003, pp. 775–780. [Online]. Available: http://www.
kecl.ntt.co.jp/icl/signal/ica2003/cdrom/index.htm[2] K. Matsuoka, “Independent component analysis andits applications to sound signal separation,” in Proc.
IWAENC 2003, Kyoto, Sep 2003, pp. 15–18.

[8] P. Bofill, “Underdetermined blind separation of delayedsound sources in the frequency domain,” Neurocomputing, vol. 55, no. 3-4, pp. 627–641, Oct 2003.
[9] S. Araki, S. Makino, A. Blin, R. Mukai, and H. Sawada,“Blind separation of more speech than sensors with lessdistortion by combining spareseness and ica,” in Proc.
IWAENC 2003, 2003, pp. 271–274.
[10] A. Blin, S. Araki, and S. Makino, “Blind source separation when speech signals outnumber sensors using
a sparseness - mixing matrix estimation (SMME),” inProc. IWAENC 2003, 2003, pp. 211–214.
[11] K. Waheed and F. M. Salem, “Algebraic overcomplete
independent component analysis,” in Proc. ICA 2003,
2003, pp. 1077–1082.
[12] F. Theis, “Mathematics in independent component
analysis,” Ph.D. dissertation, University of Regensburg, 2002. [Online]. Available: http://homepages.
uni-regensburg.de/ thf11669/phdthesis.html[13] T. Hastie, R. Tibshirani, and J. Friedman, The elementsof statistical learning: data mining, inference, and prediction, ser. Springer Series in Statistics. SpringerVerlag, 2002.
[14] L. Vielva, D. Erdogmus, and J. C. Principe, “Underdetermined blind source separation using a probabilistic
source sparsity model,” in Proc. ICA 2001, 2001, pp.
675–679.
[15] F. Murtagh, “Comments on ‘Parallel algorithms for hierarchical clustering and cluster validity’,” IEEE Trans.
on Pattern Analysis and Machine Intelligence, vol. 14,
no. 10, pp. 1056–1057, Oct 1992.