International Workshop on Acoustic Echo and Noise Control (IWAENC2003), Sept. 2003, Kyoto, Japan

BLIND SOURCE SEPARATION FOR CONVOLUTIVE MIXTURES BASED ON
COMPLEXITY MINIMIZATION
Sébastien de la Kethulle de Ryhove

Ryo Mukai

Hiroshi Sawada

Shoji Makino

NTT Communication Science Laboratories, NTT Corporation
2–4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619–0237, Japan
delaketh,ryo,sawada,maki  @cslab.kecl.ntt.co.jp
ABSTRACTtion of the cost function
f'gihZjhnj
(2)
43e]lkm7E? 43lkCoqpr Ets e] lu ,j43Cvs : & 43lu is assumed to be inwhere the separating system
vertible and e]Ov<  , ! <      is ther vector of separated
signals, its components <   being of length for all 8 .
When using the complexity minimization criterion, no assumptions about the distribution of the source signals need to be
made. Moreover, both time correlations and higher order statistics are taken into account to perform BSS, as opposed to standard
ICA-based methods which use only one of these two criteria (see[5] for more details).

Using algorithmic complexity to perform blind source separation(BSS) was first proposed by Pajunen. This approach presents
the advantage of taking the whole signal structure into accountto achieve separation, whereas standard ICA-based methods only
use either time-correlations or higher order statistics in order todo so. Another advantage of this approach is that no assumptions
about the probability distribution of the source signals need to be
made. However, although algorithmic complexity based methods
have been shown to outperform standard ICA algorithms in the
instantaneous BSS case, they haven’t been applied to convolutive
BSS to the present date. In this paper, we show that it is also
possible to use algorithmic complexity as a separating criterion to
perform BSS for convolutive mixtures and suggest a method to doso. Testing the proposed method by computer simulation yielded
results which are encouraging in terms of SNR performance.

2. PROPOSED METHOD
We consider here the estimation of only one source signal w  ,with xzy!   . In the remainder of this section, we denote by  p the source signal to be estimated(the index x isr
dropped for more clarity), <I is the length estimate for  ,and {43G|:}43! :~"43  denotes the separating FIR fil"-@ 
ters. We hence have < X * &  * />#? :C&'43B'&X5t3 .
In the case of the above one-unit algorithm, the unmixing system {43 is not invertible and the cost function (2) cannot be used
due to the presence of the first term, which we choose here to ignore (this is also done in [6] in the linear BSS case). The costfunction (2) then becomesfg s <I lu# (3)
pr Ets <Ilu4

1. INTRODUCTION
Blind source separation (BSS) is a method for recovering a setof source signals from the observation of their mixtures without any prior knowledge about the mixing process. We consider
here the convolutive mixture case, i.e.  source signals      are mixed and the corresponding mixtures  !! #"$  observed at % sensors, following,,,
,.'&()+* *0/21(& 43 6573 , where 1(& 43 represents thegoal is to find a
impulse response from source 8 to sensor 9 . The,
separating system consisting, of FIR filters : &'43 oflength ; to-@  ,
"produce separated signals < = * &  *0/>#? : &A43B'&(C5D3
,that are as close as possible to the source signals  .
Many methods based on independent component analysis(ICA) have been proposed to solve the above BSS problem. We
do not review them here, extensive literature on the subject having been published over the past years (see for instance [1, 2, 3]).
We however give an outline of the less well known algorithmic
complexity based approaches.
The Kolmogorov (or algorithmic) complexity EGFCIH6 of astring  H J  IKL with respect to a universal computer M
H
is defined as ([4, 5])
ENFC H OQP
GWZYRTS []\'^ _ F` H  ^ 
(1)
UVX

In order to actually perform BSS, we first need to derive an approximative expression Et
 s <Ilu for the algorithmic complexity
Ets <Ilu of a time sequence <I . This is the object of sec. 2.1(a similar discussion can also be found in [6]). Finally, we use thecost functionfg  s <I lu# (4) s < lu
pr Et
to implement BSS. Sec. 2.2 is devoted to the calculation of the gradient of (4) using the approximation Et
 s <Ilu for Ets <Ilu obtained
in sec. 2.1. In sec. 2.3, we discuss how to impose a normalization
contraint on {43 in order to avoid solutions of the form {43C7 .
2.1. Algorithmic Complexity of a Time SequencerConsider a time sequence <IOv<4! <Ir  5
  . It can beshown (see [4]) that if its samples <I4< 5  p are i.i.d. acp
cording to the probability density function _ < of a random vari
able  , thenthe minimum length over all programs _ FC H  on a universal
computer M that print aH and halt. The universal computer M
is henceforth considered fixed and its mention thus omitted (see[4] for more extensive explanations). The BSS problem can thenbe defined as [5]:Havingobserved mixtures  Ccb(4  , findd a mixing mappingdb such that the total complexity of the mapping b and the separated
signals e` is minimized.
Applying this definition to convolutive BSS leads to the minimiza-

RnP 
(5)
pr Ets < luI7ts <Ilu4 #
where the notation ts < lu is used for the entropy t of therandom variable  (for convenience, a similar notation will be
used in all expressions involving expectations in the remainder ofthis paper).

303

However, if <I  is a natural signal, there will in general bedependencies among its samples, and the above theorem cannot
be used. Therefore, in order to help remove these dependencies,
we use linear prediction (see [7]) to find an estimate <#  for < 
following<# O{43 of length ; such that the cost function given in (4),f g 
(13)s <Ilua E   s <Ilu4
where E  s <Ilu is evaluated using (12), is minimized. We thusf 
calculate the derivative * .W - \ g s <Ilu in order to perform gradif g *,+
ent descent on
s <Ilu . After some basic algebraic manipulations, we obtain:/ f g /s <Ilu
>a ?
/
/
E
{  433   L5 35 4i
 1 / -@{ 0
{ 0 2(6)
x <IL
5 #
and define the sequence of the residuals 2<I of <I as@r
 <I  <I 5 <#
C  ! 5 (7) 
p
the idea being to approximate the complexity of the sequence< by that of the sequence  <I . Note that in order to perform linear
prediction, the process <  must be wide-sense stationary. @ Moreover, if <I  can be modeled as an AR process of order  , thenthe sequence 2<  is white and equal to the innovations processof <I . If its samples are also i.i.d. according to the probability
density function _   <  of a random variable 2 , then the equality(8)
RnP r p @ Ets  <I lu#7ts  <I lu
5 

a$
holds, where ts  <Ilu denotes the entropy of the random variable2 . However, it is difficult to find natural sequences <I that canexactly be modeled as AR processes, and even in that case, there is
no guarantee that the samples of the innovations process 2<I will
be mutually independent. Equality ( ) is thus only approximative.
Nevertheless, for natural signals <I  , it is reasonable to assume
that (8) is closer to being verified than (5) since linear prediction,
as used in (6) and (7), helps remove dependencies. The complexityof the sequence <I is thus approximated byrr@ Ets 2< l
(9)Ets <Iluruts  <I lu 
5
where the second approximation was obtained using (8). Notethat the complexity of the prediction coefficients xLx  
has been neglected in (9).
Let us now consider the problem of the estimation of the entropy of the residual ts 2<Ilu which appears in (8) and (9). It isuseful to normalize the residual 2< to unit variance in order to remove the effects of scaling. Denoting by  the variance of theresidual, we have2< (10)
ts  <Ilu7
o
 

-

"

o
5

p

p

 ;:

 



5

p

 

"
1



 X50768% 2<

 <$ 68%



-

 :

"
1

4

 9&

<IL5?a'68% <&

2<$

 -/ <$

 9&

4>=

/

 

{ 0
/4@/{ 0

x
=



(14)

 X5?# ,  <I ]v* />#-? @  {D`43
with   >   X5 *
 x
 h

~503m5c*
m503O5 # k ,   the variance of 2<I , x
x   x A the linear prediction coefficients obtained using <I
and 6I   the derivative of    . To actually use this expression ina gradient descent algorithm, we need to evaluate * .W - \   and@*,+* W.- \ xfor all y BBC   . The expressions@ for xC  x ,p*,+and ' being quite complex (especially for large  ), this is not aneasy task.
@(first order
We hence restrict ourselves to the case linear prediction). In this case, the expressions for p * .W - \   and*,+* .W - \ x] (which are given in the appendix) are reasonably simple.
*,+

2.3. Normalization Constraint
To avoid solutions of the kind {43C0 , we need to impose a nor
malization constraint on {43 . We propose two different possibili ties. A first option consists in adding a term of the form 58D
f 
to g s <Ilua E  s <Ilu , which yields the cost  functionf F
g   E G'H s <Ilua E  s < lu
(15)
5 D
  

, this increases the cost of solutions with small values of
For DIK  (which should be avoided). Note also that the value of the cost
Kfunction is left almost unchanged if  . Adapting the valuepKof D so as to have  using a simple control loop proved topwork well in practice: a minimumof the cost function satisfyingK 
was always attained (see sec. 4).
pComputersimulations also showed that, at least in the case ofspeech mixtures, performing constraint free gradient descent on
f g 
s < luO E   s <Ilu actually seemed to work as well: althoughKboth   and the amplitude of {43 decrease during minimization,Kthe latter is completed before either   or {43 become unacceptably small. This is our second method.
In the sequel, we respectively refer to our two methods of dealing with the normalization problem as the “control loop” and “unconstrained” methods.

 

If a good approximation (denoted   ) of the negative log
densityof the pdf of the residual normalized to unit varianceh 
5ts  <I lu as_  ! lk is known, we can easily approximate

2<
ts  <Ilu #"$ %
o
 
(11)
 

 !

p

1

 '&(In most cases, _ is unknown and we set   equal to thenegative log density of a generic super- or subgaussian randomvariable (depending on the nature of the residual). We finally obtain the following approximation for the per sample complexity
Es <Ilu   Ets <Ilu of the time sequence <I   :


2<
(12)
E   s < lu#)"$ %o'  &(2.2. Gradient Descent
Now that an approximation for the algorithmic complexity of asequence <  has been derived, let us return to the problem ofblind source separation of convolutive mixtures, where we ex"-@ 
tract one component <I
* &  *0/>#? : & 43B & 5v30-I@ * />a? { C43   5c3 , with {43 : 43! :~"$43  and C   !  "     . We now have to find the FIR filters

3. CONVOLUTIVE SPEECH MIXTURESIn this section, we show how to apply the method outlined in sec. 2
to the particular case of BSS for convolutive speech mixtures. We
thus assume that sequences <I and   represent speech signals
sampled at a given frequency JLK .

304

1
Normalized Residual (histogram)
Laplacian pdf L(δ y / σ )δTable 1. Experimental conditions
Source signal length #$Direction of sources
' and ( )' (two sources)&p %Inter-sensor spacing r*. ,+$-Pi (two sensors)
Reverberation time
XP
p
Sampling
frequency
J K #0/2143 h variance 5_  !  k . Evaluating the probability density _  ! using histograms (see Fig. 1) shows that, at least in the case of 
the parametersspecifiedin sec. 4, it is reasonable to assume that
h 


5


o 5
6
k
!
! , i.e. the negative log-density
K
_

0.9

0.8

0.7

δ

p(δ y / σ )

0.6

0.5

0.4

0.3

0.2

0.1

0
−5

−4

−3

−2

−1

.

0
δ y / σδ

1

2

3

4

%

5

Fig. 1. Probability density function of the residual (normalized.
to unit variance) _  2<    as estimatedusing histograms. The
K   ! 
KLaplacian pdf ;~ 2
is a good approxi<   6 K  ? mation for _  2< ' .
3.1. Stationarity and AR model



7

7

The properties of speech have already been extensively analyzed:
it has among others been shownr that it can be considered stationary during time intervals of
tP [8] and that it canaccurately be modeled as an AR process [9], the order of which
depends on the particular application in view.
In order to compute the approximation of the per sample complexity E  s <Ilu of <I , we need to use linear prediction, with theunderlying assumptionthat <I is wide-sense stationary. Therer r DP
fore, since  , we split <I into  segments of)J K 
followinglength ;
  ,
 ?
W \
(16)
<IX , -I@ <
where
,
W \
<
O%87

7 7 9i:5 . 1 This is not a surof a Laplacian random variable with% speech is Laplacian fact thatprising result since it is a well knownh 
distributed. Note however that 5!  k is function of many_different parameters (mixing system, separating system, order of 
the linear predictor, source signals, iteration number,...) and in
some cases may not be well approximated by K
o 5
6
! .
%
%87 7
7

4. COMPUTER SIMULATIONS
The method outlined in secs. 2 and 3 was tested by computer
simulation. One source signal was estimated from two twosignal mixtures, this having been done for twelve different source
signal combinations. The experimental procedure , was the following: the mixing system’s impulse responses 1 & 43 8 l9  were first measured in a real room using the experimenp conditions%
talsummarized in Table 1. Twelve different combinations   K   of speech signals produced by two maleand two female
speakers
were then mixed following  & 7
,
,
,.K 1
4


3



=

5

3

,
and the proposed one-unit algorithm
*
*
&
/

subsequently applied to extract one speech signal from each
@ of the
twelve mixtures. We
used
a
first
order
linearpredictor(),
r 
(XP and separating filters {43mp
a stationarity time:=43 : K 43  of length ; &p %(# taps. Although speech sig%
nals can only be considered;stationary
for time intervals of apr proximately
XP , we set(XP  (corresponding to&p % to be able to learn fil; taps at J K$/21<3 ) in orderters {43 p of acceptable
lengthr (this being because we ideally must

; = ;
(X
P  proved to be accept). Using
have >
p %able for our complexity evaluation purposes.
Moreover, {43 was
initially set to 
{ ? @A? 43mCB5  > K (5 A > K  D , we used a step size of FE? G , and two different methods (“control loop” and “unppconstrained”)
were used to deal with the normalization problem(see sec. 2.3).
The results are displayed in Table 2. The SDR (between
the original signal  and the extracted signal <  ) was measured using the method proposed in [10]. Both the “control loop”
and “unconstrained” versions of the suggested algorithm achieved
a good SNR improvement, but the extracted signal was often
severely distorted. In the “control loop” version, the distortion was
due to filtering (both high-pass filtering and low-pass filtering wereobserved), whereas in the “unconstrained” version it was due towhitening (see Fig. 2). We observed that the distortion introducedby high- or low-pass filtering was much more unpleasant than the
one introduced by whitening when it came to listening. Therefore,
we believe that the “unconstrained” version is more suitable thanthe “control loop” version for convolutive speech mixture BSS.
Table 2 also shows that the algorithm performance depends on thenature of the input signals: for instance, in the “control loop” ver-
r
if 8(;and c8Io A;potherwise.
,
(17)
W \
,
If the  segments <independent,
itisreason are mutually
W \, -I@ ? Ets <
 lu (see [5]). The
able to assume that Ets < lu*
gradient of E  s <I lu can then be approximately evaluated as

 
/
/
,
 ?
W \
/
/
(18)
E   s <IluE<" !A
, -I@{ 0
{ 0,   W \  , proceeding as with
After splitting   into  segments <I,  in (16) and (17), equation (14) – where    is replaced byW\ – can be used to evaluate each one of the terms in the above,,
sum. Linear predictionis also performed separately
for each of theW \W \ segments <  yielding  sequences 2<  which are usedinstead of  <I  in (14). Note however that the equality
,

,
 ?
W \
>a
W \
<
 C -@ {  43 
L5 3
(19)
/


 ;
is not quite exact for values of  close to 8 ;or 8o '.
p,Since speech can be modeled as an AR
process, we can reaW \sonably assume that the sequences  <  are close to being,
white, and hence that their samples are close to being uncorrelated.
W \This is a necessary condition for having i.i.d. sequences 2, < .
W \,  , isWe can thus hope that (8), where 2<I is replacedby 2<W \ .
approximately verified for each sequence 2<<I

3.2. Approximation of the pdf of the Residual

1 The generic probability density function of a zero mean Laplacian ran-

In this section, we consider the problem of the evaluation of thenegative log density of the pdf of the residual normalized to unit

dom variable X is given by HJIKMLON,P KQ ? P 

305

Y 

.

M1-M2
M2-M1
M1-F1
F1-M1
M1-F2
F2-M1
M2-F1
F1-M2
M2-F2
F2-M2
F1-F2
F2-F1

M1
M2
M1
M1
M1
M1
M2
M2
M2
M2
F1
F1

Average

SNR
imp.
(dB)

 . 
p .#
p (p . ( % . # p . (p . #% . #p  .
p + . %)%( .p.#
p + . p
% . (p 
p

Ext.
Sig.

SNR
imp.
(dB)

5 # .%5 p  .)5<( .()(

M1
M1
M1
F1
M1
F2
F1
F1
F2
F2
F1
F2

% .# p%% . %
p .++%
 . 
# . 


5  .
5

.
5 #% . p&(%5 # .+ 
5 p . p5 p ( .p
5 p # .)p# p
5 p .
5 p . %(p&p % . %5p pControl Loop (M2−M1)

UNCONSTRAINED

Final
SDR
(dB)


 .

p .

% . #% # .+)

% .
%'p . %
p

p . 

# . (

Power Spectrum Magnitude [dB]

Ext.
Sig.

Final
SDR
(dB)


5<( . (5 + .
5 + .p 
5 ( . %#
5p + . #5 ( .p (5 p .p5 . (p5. #
5 p p .#
5 p . ( p5 % . p
p 
5 .

%

.

−40
−60
−80
−100

Original (M2)
Output
0

500

1000

1500

2000
Frequency [Hz]

2500

3000

3500

4000

−10
−20
−30
−40
−50
−60
−70

Original (F2)
Output
0

500

1000

1500

2000
Frequency [Hz]

2500

3000

3500

4000

Fig. 2. Power Spectrum Magnitude of the original and extracted
signals for the “control loop” (Speaker Combination: M2-M1, top)
. “unconstrained” (Speaker Combination: F2-M1, bottom) verandr.
sions of the proposed algorithm.
 CP
p
The derivatives * .W - \   and * .W - \ x  respectively readsion, the results were much better for combination M1-M2 than for
combination F1-F2 .
The fact that the SDR be negative might seem strange at firstsight. Note however that we are trying to recover the original signal  , which is substantially more difficult than recovering areverberated version of (  (observed at one of the sensors), as inmost traditional convolutive BSS algorithms based on ICA. Bearing this in mind, an SDR of 5 + to <5 r ( .  – comparable to theOP – is acceptdistortion due to the reverberation timeable. Nonetheless, when the SDR is below 5 p  , the extractedpsignal, especially in the “control loop” version of the algorithm,
can become very unpleasant to hear. Still, we believe that algorithmic complexity is a concept powerful enough to allow to recover
the original signal   and will attempt, in future work, to improve the SDR performance by using japproximations of the costfunction (2) which include the term Ets ? lu (neglected in thispaper).

*,+

/

 ?
>#
-I@
/



*,+

/ U { / 5 XK p   p?
{ 0
 :
4


 ?
 ?
>#
>#
/ U  o / U  A{ / 5X  -I@
X4 -I@ 
?
? ?
p /
/
/

and
/
/

{ 0

x]]Kp/ U { / =
?



 ?
>#
X4 - @  / ? U  o / ? U ?  X{ /4 :/
 ?
>a/ U { /= 5  -@
?
%
p /

(20)

(21)

B. REFERENCES[1] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component
Analysis, John Wiley & Sons, 2001.
[2] A. J. Bell and T. J. Sejnowski, “An information-maximization approach to blind separation and blind deconvolution,” Neural Computation, vol. 7, no. 6, pp. 1129–1159, 1995.
[3] K. Matsuoka, M. Ohya, and M. Kawamoto, “A neural net for blindseparation of nonstationary signals,” Neural Networks, vol. 8, no. 3,
pp. 411–419, 1995.
[4] T. Cover and J. Thomas, Elements of Information Theory, John
Wiley & Sons, 1991.
[5] P. Pajunen, “Blind source separation using algorithmic information
theory,” in Proc. of Independence and Artificial Neural Networks
Workshop (I&ANN’98), Feb. 1998, pp. 26–31.
[6] A. Hyvärinen, “Complexity pursuit: Separating interesting components from time-series,” Neural Computation, vol. 13, no. 4, pp.
883–898, 2001.
[7] J. Proakis and D. Manolakis, Digital Signal Processing, Prentice
Hall, 1996.
[8] A. Oppenheim and R. Schafer, Discrete-time Signal Processing,
Prentice Hall, 1998.
[9] D. Manolakis, V. Ingle, and S. Kogon, Statistical and Adaptive Signal Processing, McGraw-Hill, 2000.
[10] R. Gribonval, E. Vincent, C. Févotte, and L. Benaroya, “Proposals
for performance measurement in source separation,” in Proc. ICA
2003, April 2003, pp. 763–768.

5. CONCLUSIONSWe proposed a method to use algorithmic complexity as a separating criterion to perform BSS for convolutive mixtures. The main
advantage of this method over standard ICA algorithms is that, toachieve separation, the whole signal structure is taken into accountinstead of only time correlations or higher order statistics. The experimental results obtained thus far are encouraging in terms of
SNR performance although the extracted signal is often severely
distorted. We believe that this problem can be solved by using
approximationsof the cost function (2) which include the termj

Ets ? lu . This will be the subject of future research.

A. APPENDIX

U−20

Unconstrained (F2−M1)

Table 2. SNR improvement and final SDR obtained with the pro.posed algorithm (“control loop” and “unconstrained” versions) forr .
 CP
twelve different speech signal combinations.
p

If we define X _ 
"Gs<I o  @ 
 W\ u and_ W \ @ <IKand x "Gs o?W \_   (u , we have   

0Power Spectrum Magnitude [dB]

CONTROL LOOP
Speaker
Combin.



W @ \

W \.

306