Algorithmic Complexity based Blind Source Separation
for Convolutive Speech MixturesSé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

Abstract
It has been recently shown by Pajunen and Hyvärinen that using
algorithmic complexity based cost functions to perform blindsource separation (BSS) yields very good results, at least in the
instantaneous BSS case. From the theoretical point of view,
such cost functions present numerous advantages over the ones
derived using the theory of ICA. In this paper, we suggest amethod of using algorithmic complexity to perform blind source
separation for convolutive speech mixtures. After deriving theappropriate cost function, we show how linear prediction can beused to obtain an acceptable approximation for the algorithmic
complexity of a speech signal. The well known properties ofspeech (stationary for time intervals of approximately 30 ms.,
possibility of accurate modeling as an AR process, etc.) aretaken into account to derive this approximation. Finally, after examining the different problems which arise when actually
implementing gradient descent on the final cost function, we
discuss the results of computer simulation which are encouraging in terms of SNR performance and propose some directions
for future work.

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. n source signals s(t) =(s1 (t), . . . , sn (t))T are mixed and the corresponding mixtures
x(t) = (x1 (t), . . . , xm (t))T observed at m sensors, following
xj (t) = n
i=1
l hji (l)si (t − l), where hji (l) represents theimpulse response from source i to sensor j. The goal is to find a
separating system consisting of FIR filters wij (l) of length L to
L−1
produce separated signals yi (t) = m
j=1
l=0 wij (l)xj (t −l) that are as close as possible to the source signals si (t).
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 KU (xN ) of astring xN = x1 x2 . . . xN with respect to a universal computer
U is defined as ([4, 5])(1)
KU (xN ) = min |pU (xN )|,

 

 

pU (xN )the minimum length over all programs pU (xN ) on a universal
computer U that print xN and halt. The universal computer U
is henceforth considered fixed and its mention thus omitted (see[4] for more extensive explanations). The BSS problem canthen be defined as [5]:Th.P1.8

Having observed mixtures x(t) = f (s(t)), find a mixing mapping f̂ such that the total complexity of the mapping f̂ and theseparated signals y(t) is minimized.
Applying this definition to convolutive BSS leads to the minimization of the cost function
1(2)
JK W(l), y(t) = K W−1 (l) + K[y(t)],Twhere the separating system W(l) = [wij (l)] is assumed to
be invertible and y(t) = (y1 (t), . . . , yn (t))T is the vector of
separated signals, its components yi (t) being of length T for alli.
When using the complexity minimization criterion, no assumptions about the distribution of the source signals need tobe made. Moreover, both time correlations and higher order
statistics are taken into account to perform BSS, as opposed tostandard ICA-based methods which use only one of these two
criteria (see [5] for more details).









2. Proposed Method
We consider here the estimation of only one source signal sα (t),with α ∈ {1, . . . , n}. In the remainder of this section, we denote by s(t) the source signal to be estimated (the index α is
dropped for more clarity), y(t) is the length T estimate for s(t),and w(l) = (w1 (l), . . . , wm (l)) denotes the separating FIR filL−1
ters. We hence have y(t) = m
j=1
l=0 wj (l)xj (t − l).
In the case of the above one-unit algorithm, the unmixing
system w(l) 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
cost function (2) then becomes
1

JK
(3)
[y(t)] = K[y(t)].
T
In order to actually perform BSS, we first need to derive an approximative expression K̂[y(t)] for the algorithmic complexity
K[y(t)] of a time sequence y(t). This is the object of sec. 2.1(a similar discussion can also be found in [6]). Finally, we use
the cost function1JK̂[y(t)] = K̂[y(t)](4)
T
to implement BSS. Sec. 2.2 is devoted to the calculation of thegradient of (4) using the approximation K̂[y(t)] for K[y(t)]obtained in sec. 2.1. In sec. 2.3, we discuss how to impose anormalization contraint on w(l) in order to avoid solutions of
the form w(l) = 0.

 

2.1. Algorithmic Complexity of a Time Sequence
Consider a time sequence y(t) = (y(0), . . . , y(T − 1)). It canbe shown (see [4]) that if its samples y(0), . . . , y(T − 1) are
i.i.d. according to the probability density function p(y) of a

IV - 3127

random variable Y , then1
K[y(t)] = H[y(t)],(5)
lim
T →∞ T
where the notation H[y(t)] is used for the entropy H(Y ) of therandom variable Y (for convenience, a similar notation will be
used in all expressions involving expectations in the remainderof this paper).
However, if y(t) 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 ŷ(t) fory(t) following

 α y(t − τ )τ0ŷ(t) =

τ

(6)

τ =1

and define the sequence of the residuals δy(t) of y(t) as
t = τ 0 , . . . , T − 1,(7)δy(t) y(t) − ŷ(t),the idea being to approximate the complexity of the sequencey(t) by that of the sequence δy(t). Note that in order to perform linear prediction, the process y(t) must be wide-sense stationary. Moreover, if y(t) can be modeled as an AR process oforder τ0 , then the sequence δy(t) is white and equal to the innovations process of y(t). If its samples are also i.i.d. according
to the probability density function p(δy) of a random variable
δY , then the equality1K[δy(t)] = H[δy(t)](8)
lim
T →∞ T − τ0holds, where H[δy(t)] denotes the entropy of the random variable δY . However, it is difficult to find natural sequences y(t)that can exactly be modeled as AR processes, and even in thatcase, there is no guarantee that the samples of the innovations
process δy(t) will be mutually independent. Equality (8) is thus
only approximative. Nevertheless, for natural signals y(t), it isreasonable to assume that (8) is closer to being verified than(5) since linear prediction, as used in (6) and (7), helps remove
dependencies. The complexity of the sequence y(t) is thus approximated byT
K[y(t)] ≈K[δy(t)] ≈ T H[δy(t)],(9)
T − τ0where the second approximation was obtained using (8). Notethat the complexity of the prediction coefficients (α1 , . . . , ατ0 )has been neglected in (9).
Let us now consider the problem of the estimation of theentropy of the residual H[δy(t)] which appears in (8) and (9).
It is useful to normalize the residual δy to unit variance in orderto remove the effects of scaling. Denoting by σδ the variance ofthe residual, we haveδy(t)+ log σδ(10)H[δy(t)] = Hσδ
If a good approximation (denoted G(·)) of the negative log
density of the pdf of the residual normalized to unit variance− log p( σδy ) is known, we can easily approximate H[δy(t)]δasδyH[δy(t)] ≈ E G+ log σδ .
(11)σδIn most cases, p( σδyδ ) is unknown and we set G(·) equal to thenegative log density of a generic super- or subgaussian randomvariable (depending on the nature of the residual). We finallyobtain the following approximation for the per sample complexity KT [y(t)] T1 K[y(t)] of the time sequence y(t):δy+ log σδ .
(12)K̂T [y(t)] = E G
σδ













    

2.2. Gradient Descent
Now that an approximation for the algorithmic complexity ofa sequence y(t) has been derived, let us return to the problemof blind source separation of convolutive mixtures, where we
L−1extract one component y(t) = m
j=1
l=0 wj (l)xj (t − l) =L−1
TT
l=0 w (l)x(t − l), with w(l) = (w 1 (l), . . . , wm (l)) andTx(t) = (x1 (t), . . . , xm (t)) . We now have to find the FIR
filters w(l) of length L such that the cost function given in (4),[y(t)] = K̂T [y(t)],
(13)
JK̂

 



where K̂T [y(t)] is evaluated using (12), is minimized. We thus∂
calculate the derivative ∂w(q)JK̂
[y(t)] in order to perform gra
dient descent on JK̂ [y(t)]. After some basic algebraic manipulations, we obtain:
L−1
[y(t)]∂JK̂
∂=w T (l) · x(t − l) =
K̂T
∂w(q)
∂w(q)
l=0



E δx(t − q) g
+

1
σδ

1−

1
E δy (t) gσδ δy (t) σδ δy (t) 
σδ

∂
σδ
∂w(q)

 E y(t − τ ) g δy (t)  · ∂ α , (14)
σ
∂w(q)

 w (l)x(t)−α x(t−τ), δy(t) =withx(t δx(t)− l) −α x(t − l − τ ) , σ the variance of δy(t),τ0−1σδ τ =1

τ

δ

τ0
τ =1

τ0
τ =1

L−1
l=0

τ

τ

T

δ

α1 . . . ατ0 the linear prediction coefficients obtained using y(t)and g(·) the derivative of G(·). To actually use this expression∂σδ
in a gradient descent algorithm, we need to evaluate ∂w(q)∂and ∂w(q) ατ for all τ ∈ {1, · · · , τ0 }. The expressions forα1 . . . ατ0 and σδ being quite complex (especially for large τ0 ),this is not an easy task.
We hence restrict ourselves to the case τ0 = 1 (first order∂linear prediction). In this case, the expressions for ∂w(q)σδ∂and ∂w(q) α1 (which are given in the appendix) are reasonably
simple.
2.3. Normalization Constraint
To avoid solutions of the kind w(l) = 0, we need to impose
a normalization constraint on w(l). We propose two different
possibilities. A first option consists in adding a term of the form[y(t)] = K̂T [y(t)], which yields the cost func−µ log σy to JK̂tion
[y(t)] = K̂T [y(t)] − µ log σy .
(15)
JK̂,CLFor µ > 0, this increases the cost of solutions with small valuesof σy2 (which should be avoided). Note also that the value ofthe cost function is left almost unchanged if σy2 ≈ 1. Adapting
the value of µ so as to have σy2 ≈ 1 using a simple control loop
proved to work well in practice: a minimum of the cost function
satisfying σy2 ≈ 1 was always attained (see sec. 4).
Computer simulations also showed that, at least in the caseof speech mixtures, performing constraint free gradient descent
on JK̂
[y(t)] = K̂T [y(t)] actually seemed to work as well: although both σy2 and the amplitude of w(l) decrease during minimization, the latter is completed before either σy2 or w(l) become unacceptably small. This is our second method.
In the sequel, we respectively refer to our two methods ofdealing with the normalization problem as the “control loop”
and “unconstrained” methods.

IV - 3128

3. Convolutive Speech Mixtures

1
Normalized Residual (histogram)
Laplacian pdf L(δ y / σδ )0.9

In this section, we show how to apply the method outlined insec. 2 to the particular case of BSS for convolutive speech mixtures. We thus assume that sequences y(t) and x(t) represent
speech signals sampled at a given frequency fs .

0.8

0.7

p(δ y / σδ )

0.6

0.5

0.4

3.1. Stationarity and AR model

0.3

The properties of speech have already been extensively analyzed: it has among others been shown that it can be considered
stationary during time intervals of Tstat ≈ 30 ms. [8] and thatit can accurately be modeled as an AR process [9], the order ofwhich depends on the particular application in view.
In order to compute the approximation of the per sample
complexity K̂T [y(t)] of y(t), we need to use linear prediction,with the underlying assumption that y(t) is wide-sense stationary. Therefore, since Tstat ≈ 30 ms., we split y(t) into M
segments of length L stat = fs · Tstat following y (t),M −1y(t) =

(i)

0.2

0.1

0
−5

where

y(t)
0









After splitting x(t) into M segments x(i) (t), proceeding aswith y(t) in (16) and (17), equation (14) – where x(t) is replaced by x(i) (t) – can be used to evaluate each one of theterms in the above sum. Linear prediction is also performed
separately for each of the M segments y(i) (t) yielding M sequences δy(i) (t) which are used instead of δy(t) in (14). Notehowever that the equality
y (i) (t) = w (l)x (t − l)

L−1

T

(i)

(19)

l=0

is not quite exact for values of t close to i·L stat or (i+1)·L stat .
Since speech can be modeled as an AR process, we can reasonably assume that the sequences δy(i) (t) are close to being
white, and hence that their samples are close to being uncorrelated. This is a necessary condition for having i.i.d. sequences
δy (i) (t). We can thus hope that (8), where δy(t) is replaced by
δy (i) (t), is approximately verified for each sequence δy(i) (t).
3.2. Approximation of the pdf of the Residual
In this section, we consider the problem of the evaluation of
the negative log density of the pdf of the residual normalized
to unit variance − log p( σδy ) . Evaluating the probability den-





δ

sity p( σδy ) using histograms (see Fig. 1) shows that, at least
δ





 

in the case of the parameters specified in sec. 4,√it is reasonable to assume that − log p( σδy ) ≈ 12 log 2 + 2 σδy , i.e.
δδ
the negative√ 1 log-density of a Laplacian random variable with
λ = 2. This is not a surprising result since it is a well known
1 The generic probability density function of a zero mean Laplaciane−λ|x| .
random variable X is given by p(x) = λ
2

−2

−1

.

0
δ y / σδ

1

2

3

4

5

Table 1: Experimental conditions
Source signal length 6 s.
Direction of sources 120◦ and 50◦ (two sources)
Inter-sensor spacing d = 4 cm. (two sensors)
Reverberation timeTr = 130 ms.
Sampling frequency fs = 8 kHz.

(16)

if i L stat ≤ t < (i + 1) L stat and t < Totherwise.
(17)
If the M segments y(i) (t) are mutually independent, it is rea−1(i)(t)] (see [5]).
sonable to assume that K[y(t)] ≈ M
i=0 K[yThe gradient of K̂T [y(t)] can then be approximately evaluated
as
M −1
∂
∂
(18)
K̂T [y(t)] ≈K̂T y (i) (t) .
∂w(q)
∂w(q)
i=0

−3

Figure 1: Probability density function of the residual (normalusing histograms.
ized to unit variance) p(δy/σδ ). as √estimated√
The Laplacian pdf L(δy/σδ ) = 22 · e− 2 |δy/σδ | is a good
approximation for p(δy/σδ ).

i=0

y (i) (t) =

−4





fact that speech is Laplacian distributed. Note however that− log p( σδyδ ) 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 wellapproximated by 12 log 2 + 2 σδyδ .

 

4. Computer Simulations
The method outlined in secs. 2 and 3 was tested by computer simulation. One source signal was estimated from two
two-signal mixtures, this having been done for twelve different source signal combinations. The experimental procedure
was the following: the mixing system’s impulse responses{hij (l), i, j = 1 . . . 2} were first measured in a real room using
the experimental conditions summarized in Table 1. Twelve different combinations (s1 (t), s2 (t)) of speech signals producedby two male and two female speakers were then mixed following xj (t) = 2i=1 l hji (l)si (t − l), and the proposed oneunit algorithm subsequently applied to extract one speech signal
from each of the twelve mixtures. We used a first order linear
predictor (τ0 = 1), a stationarity time Tstat = 125 ms. and separating filters w(l) = (w1 (l), w2 (l)) of length L = 256 taps.
Although speech signals can only be considered stationary fortime intervals of approximately 30 ms., we set Tstat = 125 ms.
(corresponding to L stat = 1000 taps at fs = 8 kHz) in order to be able to learn filters w(l) of acceptable length (thisbeing because we ideally must have L  L stat ). Using
Tstat = 125 ms. proved to be acceptable for our complexity evaluation purposes. Moreover, w(l) was initially set towinit (l) = δ( L2 ), −δ( L2 ) , we used a step size of 1.0 × 10−3 ,and two different methods (“control loop” and “unconstrained”)
were used to deal with the normalization problem (see sec. 2.3).

 





The results are displayed in Table 2. The SDR (between theoriginal signal s(t) and the extracted signal y(t)) 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

IV - 3129

Ext.
Sig.

SNR
imp.
(dB)

Final
SDR
(dB)

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

31.73
11.60
25.59
31.60
13.58
20.60
10.96
13.22
4.71
5.77
17.61
24.13
17.59

−3.29
−16.88
−3.03
−5.55
−2.12
−6.05
−16.41−17.19−15.11−16.36−18.21−12.92−11.09Average

Ext.
Sig.

SNR
imp.
(dB)

Final
SDR
(dB)

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

23.91
22.62
17.82
8.44
3.78
16.09
20.09
20.86
26.43
21.82
1.79
8.89
16.05

−5.59
−4.18
−4.23
−15.36−4.16−15.15−8.18−7.57−11.63−10.61−2.51−17.98−8.92

Table 2: SNR improvement and final SDR obtained with theproposed algorithm (“control loop” and “unconstrained” ver.sions) for twelve different speech signal combinations. Tr =130 ms.
severely distorted. In the “control loop” version, the distortion
was due to filtering (both high-pass filtering and low-pass filtering were observed), whereas in the “unconstrained” version it
was due to whitening (see Fig. 2). We observed that the distortion introduced by high- or low-pass filtering was much more
unpleasant than the one introduced by whitening when it cameto listening. Therefore, we believe that the “unconstrained” version is more suitable than the “control loop” version for convolutive speech mixture BSS. Table 2 also shows that the algorithm performance depends on the nature of the input signals: for instance, in the “control loop” version, the results weremuch better for combination M1-M2 than for combination F1-F2 .
The fact that the SDR be negative might seem strange atfirst sight. Note however that we are trying to recover the original signal s(t), which is substantially more difficult than recovering a reverberated version of s(t) (observed at one of the sensors), as in most traditional convolutive BSS algorithms basedon ICA. Bearing this in mind, an SDR of −4 dB to −5 dB– comparable to the distortion due to the reverberation timeTr = 130 ms. – is acceptable. Nonetheless, when the SDR isbelow −10 dB, the extracted signal, especially in the “controlloop” version of the algorithm, can become very unpleasant tohear. Still, we believe that algorithmic complexity is a concept
powerful enough to allow to recover the original signal s(t) andwill attempt, in future work, to improve the SDR performanceby using approximations of the cost function (2) which include
the term K[W−1 (t)] (neglected in this paper).

5. Conclusions
We 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 isthat, to achieve separation, the whole signal structure is taken
into account instead 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 approximations of the cost function (2)which include the term K[W−1 (t)]. This will be the subject of
future research.

Power Spectrum Magnitude [dB]

Speaker
Combin.

Control Loop (M2−M1)

UNCONSTRAINED

0
−20
−40
−60
−80
−100

Original (M2)
Output
0

500

1000

1500

2000
Frequency [Hz]

2500

3000

3500

4000

Unconstrained (F2−M1).

Power Spectrum Magnitude [dB]

CONTROL LOOP

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

Original (F2)
Output
0

500

1000

1500

2000
Frequency [Hz]

2500

3000

3500

4000

Figure 2: Power Spectrum Magnitude of the original and extracted signals for the “control loop” (Speaker Combination:.M2-M1, top) and “unconstrained” (Speaker Combination: F2M1, bottom) versions of the proposed algorithm. Tr = 130 ms.

A. Appendix

If we define γ(p)E[ y(t + p) y(t) ] and Ap22(1)
E[ x(t + p) xT (t) ], we have σδ2 = γ (0)−γ
and α1 =
γ(0)
γ(1)
∂
∂
. The derivatives ∂w(q)σδ and ∂w(q)α1 respectively read
γ(0)

∂
1
σδ =
∂w(q)
σδ

 
γ(0)
(A

 A w − γ(1) ·

L−1

l−p

l

l=0

γ 2 (0)

L−1
l=0

and

 A w  (20)L−1
l−p+1 + Al−p−1 ) wl − γ(1)

l−p

l

l=0



L−1

∂
1
α1 = 2γ(0)(Al−p+1 + Al−p−1 ) wl
∂w(q)
γ (0)
l=0

 A w . (21)L−1

− 2 γ (1)

l−p

l

l=0

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,” NeuralComputation, vol. 7, no. 6, pp. 1129–1159, 1995.
[3] K. Matsuoka, M. Ohya, and M. Kawamoto, “A neural net forblind separation 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.

IV - 3130