UNDERDETERMINED SOURCE SEPARATION FOR COLORED SOURCES
Stefan Winter†,1 , Walter Kellermann‡ , Hiroshi Sawada† , Shoji Makino††

NTT Communication Science Laboratories,
NTT Corporation
2-4 Hikaridai Seika-cho Soraku-gun Kyoto
619-0237 Japan‡

Department of Multimedia Communication and Signal
Processing, University Erlangen-Nuremberg
Cauerstr. 7, 91058 Erlangen
Germany

ABSTRACT

2. PROBLEM FORMULATIONThis contribution focuses on the source separation stage as important part of underdetermined blind source separation (BSS). So
far nearly all approaches for underdetermined BSS assume independently, identically distributed (i.i.d.) sources. They completely
ignore the redundancy that is in the temporal structure of colored sources like speech signals. Instead, we propose multivariate
models based on the multivariate Student’s t or multivariate Gaussian distribution and investigate their potential for underdetermined
BSS. We provide a simple yet effective filter based on the sources’autocorrelations for recovering the sources as basis for further advances in underdetermined BSS. The challenge is estimating thefilter coefficients blindly. The experimental results support the ideathat source separation for underdetermined BSS can be reduced tothe separation of their autocorrelations.

With sp (t) ∈ R denoting the p-th source signal (1 ≤ p ≤ P ) andÃ ∈ RQ×P the mixing matrix, we obtain mixed signals xq (t) ∈ R(1 ≤ q ≤ Q, Q < P ) by
x1 (t)s1 (t)




..
..
(1)
 + n(t).
 = Ã .
.
xQ (t)sP (t)

1. INTRODUCTION
Blind source separation (BSS) describes techniques that aim at separating P signals if only Q mixed versions of the original signals
are available. The need for BSS arises for example if the signals of
simultaneous conversations are captured by several microphones.
Most BSS approaches assume that there are at least as many microphones as source signals (Q ≥ P ), which is called (over-) determined BSS.
Instead, we consider underdetermined BSS, where we have lessmicrophones than source signals (Q < P ). Only few approaches
have been proposed so far [2, 5, 17, 18], and the separation quality in terms of interference suppression and signal distortion is still
not as good as with (over-) determined BSS. This is in particular
true if wideband signals like speech signals are involved. The difficulty is that in contrast to (over-) determined BSS the solution ofunderdetermined BSS goes beyond system identification. Even if
the mixing system is fully identified, additional effort is required toseparate the mixtures.
So far, nearly all approaches to solve the latter problem assume
independently, identically distributed (i.i.d.) sources [2, 5, 17, 18].
While this assumption may serve well as a first order approximation, it completely ignores the redundancy that is in the temporal
structure of colored sources like speech signals. Therefore, we propose multivariate models, that take the temporal correlation explicitly into account and investigate their potential. We consider here
only linear, instantaneous mixtures of speech signals in the time domain. This paper concentrates on the separation and assumes thatthe mixing matrix is known.
After formulating the problem analytically in Sec. 2, we provide statistical models for Bayesian inference in Sec. 3. In Sec. 4 we
propose a closed-form and numerical approach for minimum meansquare error (MMSE) source estimation based on the posterior distribution of the sources. Section 5 elaborates more on autocorrelation estimation, which turns out to be essential for high-quality
underdetermined source separation. Section 6 presents experimental results which are discussed in Section 7.

n(t) ∈ RP denotes noise added to the sensors. Let
S p = [ sp (0)···sp (T − 1) ]T ,

1≤n≤P

(2)

denote a frame of length T of the p-th original speech signal with
autocorrelation
rp (τ ) = E{sp (t)sp (t + τ )}.

(3)(·)T denotes the transpose. We summarize the P frames of thedifferent source signals by the vector


S1


S =  ...  ∈ RP ·T(4)SP
For each source S p the frame-dependent autocorrelation matrix
Rp ∈ RT ×T is given by a symmetric Toeplitz matrix. Its first row
is defined by[ rp (0)···rp (T − 1) ](5)

We summarize the P autocorrelation matrices by


R1
0


P ·T ×P ·T
..
R=
∈R
.
0
RP
Extending the scalar elements Ãqp of Ã to diagonal matrices


Ãqp
0


T ×T
..
Aqp = 
,
∈R
.
0
Ãqp
we define an extended mixing matrix
A11 · · · A1P

Q·T ×P ·T..
..
A =  ...
∈R.
.
AQ1 · · · AQP

(6)

(7)

(8)

Similar to the source signals we define
1 The author is on leave from the Department of Multimedia Communi-

cation and Signal Processing, University Erlangen-Nuremberg

X q = [ xq (0)···xq (T − 1) ]T ,

1≤q≤Q

(9)

In our first speech model, we model each frame S p by a zeromean, multivariate Student’s t distributionp(S p ) = ST (S p |αp , Σp )Figure 1: Mixing process for underdetermined BSS with Q = 2,
P =3

(14)

with degrees of freedom αp and scale Σp . Modeling the speech
signal by a multivariate Student’s t distribution has several advantages:• It is a good approximation of real speech signals. This does notjust follow from experimental results but is also based on resultswith spherically invariant random processes (SIRPs) [3] and theunivariate Student’s t distribution [5].
• The multivariate Student’s t distribution can be expressed as amultivariate Gaussian mixture model (MGM) with an infinite
number of mixtures [1] as follows
ST (S p |αp , Σp ) =Z
N (S p |0, Rp ) · IW(Rp |αp , Σp )dRp .

(15)

and summarize the X q by



X1


X =  ...  ∈ RQ·T
XQ

(10)

This results in the compact description of the mixing process for
one frame by
X = AS + N(11)

as illustrated in Fig. 1. N ∈ RQ·T is derived from n(t) in a similar
way as X is derived from x(t).
The final goal in BSS is the estimation of signals Y that resemble the original signals S as closely as possible. Since only
the mixed signals are available, this is at best possible up to arbitrary permutation and scaling. With (over-) determined BSS it issufficient to estimate the mixing matrix A or its inverse. Since themixing matrix is square and therefore invertible (assuming that A
is non-singular), the inverse can be used to separate the signals. In
contrast to (over-) determined BSS, here the estimation of the mixing matrix A is not sufficient. Even if A is available, estimating
the original signals from the mixtures poses a problem on its own,
since A cannot be simply inverted. In the following we concentrateon estimating the original signals and assume that the mixing matrix
is known or can otherwise be estimated [2, 14, 15, 16].
3. MODELS
The most general approach for estimating the unknown source signals is based on Bayesian inference. It yields the posterior distribution p(θ|X) of the desired parameters θ (here: source signals)
by accounting for their prior distribution p(θ) and the likelihoodp(X|θ) of the observed data X (here: mixed signals). From Bayes’rule we obtain [12]p(θ|X) ∝ p(X|θ) · p(θ).

(12)

So far, nearly all approaches to underdetermined BSS assume
independent and identical prior distributions. In contrast, we propose in this section two multivariate priors tailored to speech signals. We also define the likelihood derived from the commonly assumed Gaussian sensor noise model.
In both speech models we make the common assumption thatthe speech signals are mutually independent.
p(S) =Yp(S p )

(13)p

We further assume that the speech signal is approximately stationary within a frame of appropriately chosen length T .

N and IW denote a multivariate Gaussian and inverse Wishart
distribution, respectively. According to this relation, the multivariate Student’s t distribution is the marginal distribution obtained from the joint distribution of the speech signal and its autocorrelation. This leads to a hierarchical Bayesian model witha multivariate Gaussian distribution for the speech signal S p
p(S p |Rp ) = N (S p |0, Rp )(16)

and an inverse Wishart distribution for the autocorrelation matrix Rp
p(Rp |αp , Σp ) = IW(Rp |αp , Σp ).

(17)

The expected value of Rp is [8]E {Rp } =

1
Σ−1
p .
αp − T − 1(18)

The hierarchical model will be used in Sec. 4 to obtain conditional posterior PDFs.
• The inverse Wishart distribution is the conjugate prior of a multivariate Gaussian likelihood with respect to its correlation matrix [7]. This means that the corresponding posterior distribution
is again an inverse Wishart distribution, which is also exploited
in Sec. 4.
The degrees of freedom αp and the scale Σp in (17) can either
be considered deterministic or random. We summarize all parameters in the parameter vector θ and obtain the joint prior distributionaspθ (θ)==

p(S, R, A, σ 2 , α, Σ)(19)2p(S|R)p(R|α, Σ)p(A)p(σ )p(α)p(Σ) (20)with α = [ α1 · · · αP ]T . The relation between the different
variables is illustrated in Fig. 2.
As a special case of the first speech model we propose an analytically feasible, multivariate Gaussian prior as second speechmodelp(S p ) = lim ST (S p |αp , Σp ) = N (S p |0, Rp )
αp →∞

(21)

with zero mean vector and covariance matrix Rp .
For deriving the likelihood we assume additive white Gaussian
noise with identical variance σ 2 across all sensors. Based on the
mixing process (11), this assumption results in the Gaussian likelihoodp(X|A, S, σ 2 ) = N (X|AS, σ 2 ).

(22)

4.1.2 Implementation
Scale

Degrees of
freedom

We now derive the conditional PDFs for our specific problem basedon the multivariate Student’s t speech model in Sec. 3.
For the conditional PDF of the speech signal we obtain

Correlation∝p(S|θ −S , X)

∝
∝
Noise
variancep(S|R) · p(X|A, S, σ 2 )(24)
2N (S|0, R) · N (X|AS, σ )N (S|µS , ΣS )

(25)
(26)

with

−1
1
ΣS = R−1 + 2 AT Aσ

Mixingmatrix(27)

andFigure 2: Graphical modelµS
4. BAYESIAN INFERENCE−1
σ 2 R−1 + AT A
AT X.

(28)

For the conditional PDF of the autocorrelation we obtain

The general goal in Bayesian inference is to estimate the posterior
PDF of unknown parameters based on the available data. Except for
few special cases like the Gaussian speech model (21) with Gaussian likelihood (22), the posterior cannot be derived in closed-form.
Instead, it requires numerical methods based for example on Monte
Carlo techniques [4]. Based on our multivariate models in the previous section, we derive in the following MMSE estimates of thesources in closed-form and by numerical approximation. Further
details are provided in the appendix.
4.1 Multivariate Student’s t prior
In order to estimate the sources based on the multivariate Student’s
t prior (14), we employ a Gibbs sampler, which is a Monte Carlotechnique. In contrast to our multivariate approach, the Gibbs sampler was applied in [5] to an i.i.d. speech model based on a univariate Student’s t distribution.
4.1.1 Basic idea
The Gibbs sampler approximates the posterior PDF by sampling
iteratively from appropriately chosen conditional PDFs. In fact, theconditional PDFs form a first-order Markov chain, whose stationary
PDF is the desired posterior PDF [4, 8].
Let θ = {θ1 , . . . , θK } denote the unknown parameters withprior PDF pθ (θ) and X the available data with likelihoodpX|θ (X|θ). We further define θ −i as summarizing all parameters
without parameter θ i . For the Gibbs sampler we need the conditional PDFs p(θ i |θ −i , X), which can be shown to be [5]p(θ i |θ −i , X) ∼ pθ (θ)pX|θ (X|θ)

=

(23)

Once they are derived, we can sample iteratively from each conditional PDF using the previous samples as summarized in Algorithm
1. As soon as a sufficient number of samples is drawn, they apAlgorithm 1: General Gibbs sampler
Initialize θ (e.g. randomly) ;for k=1..K do
for i do
(k+1)
(k,k+1)θi∼ p(θ i |θ −i
, x);
end
end
proximate the posterior PDF. They can further be used to calculate
point estimates Y (MAP, MMSE, ...) of the original signals S. In
order to avoid transient effects a certain number of initial samples
should be discarded [9, 13].

p(R|θ −R , X)∝∝

p(S|R) · p(R|α, Σ)Y
N (S p |0, Rp ) · IW(Rp |αp , Σp )

(29)

p

∝

Y



−1 
−1
IW Rp |αp + 1, S p S T
p + Σpp

To successfully implement this Gibbs sampler, appropriate
choices for the degrees of freedom αp and the scale Σp are necessary. It is common in Bayesian inference to manually tune thedegrees of freedom. Usually there is also sufficient prior data available to train the scale (e.g. speaker adaptation for speech recognition). Here we face the problem, that we do not have sufficient priordata. We consider the following options:• We could consider the scale as random variable as well andchoose an uninformative prior (e.g. Jeffrey’s prior, uniform
prior, ...). So far experiments with uninformative priors did not
yet lead to reasonable results. Therefore, we don’t consider thisoption here.
• Knowing that the choice of the scale Σp directly influences the
autocorrelation by (18), we can set it to a specific value, determined for example by some preprocessing (e.g. `1 -norm minimization [17]) as detailed in the next section.
4.2 Multivariate Gaussian prior
Plugging the Gaussian prior (21) and the likelihood (22) in (12) we
obtain for the Gaussian posterior similar to (24)p(S|X)∝∝

p(S|R) · p(X|A, S, σ 2 )N (S|µS , ΣS ).

(30)

with ΣS and µS given by (27) and (28), respectively.
Since the MMSE estimate of a normally distributed randomvariable is its mean, the resulting estimate of the sources is
Y

=


−1
σ 2 R−1 + AT A
AT X.
|
{z
}

(31)

=W

The unmixing filter W is a generalized inverse modified by the
autocorrelation matrix. The influence of the autocorrelation matrix
R is determined by the noise variance σ 2 . Without modificationby the autocorrelation matrix (31) could not be determined, since
AT A is a singular matrix.

=

5. AUTOCORRELATION
For the proposed speech models the autocorrelation plays a crucialrole. With the multivariate Gaussian model, the estimate is useddirectly. With the multivariate Student’s t model, the estimate canbe used to set the scale Σp , if it is treated as a deterministic value.
It then determines the expected value of the autocorrelation. Therefore, for high quality separation in terms of interference suppression
and signal distortion [6] we need a good estimate of the autocorrelation. In this section we will discuss possible options for estimating
the autocorrelation. In general, we can distinguish between nonparametric and parametric approaches.

r(τ )

=
≈

E {s(t)s(t + τ )}
1
T

T −τ
X−1

s(t)s(t + τ )

(32)
(33)

Non-parametric estimation of the autocorrelation suffers from an insufficient amount of data compared to the number of parameters (ifwe define all time lags of the autocorrelation as parameters). This
results in high variance and low resolution [11]. In general, it is desirable to reduce the number of parameters by an appropriate parametric model that exploits the inherent structure.
Based on the harmonicity of speech we propose approximating
the biased autocorrelation by M harmonic, decaying cosine functions.
!
MX
rp (τ ) ≈γm,p cos(p · f0,p · τ ) exp(−βp τ )
(34)
m=1

Here, γm,p denotes the gain of each cosine, f0,p is the fundamental
frequency and βp is the decay factor. This reduces the number ofparameters to M + 2 for each source.
5.2.1 Least square estimation
If we have access to the autocorrelation that we want to approximate, we can use a least square approach to estimate the parameters.
We apply this method to the following scenarios.
• To explore the possibilities of the parametric approach, we canuse the autocorrelations of the original signals to estimate theparameters.
• Other parameter estimation methods might need to be initialized. Then we can use initial estimates of the sources (e.g. `1 norm minimization [17]) and use their autocorrelations.
• Linear mixtures of mutually independent signals have the property, that their autocorrelations are a linear mixture of the autocorrelation of the original signals.
=
=

=
=

E {xq (t)xq (t − τ )}
(35)

 

 X

X
E 
Ãqp sp (t) · Ãqp sp (t − τ )


p
p


X

E
Ã2qp sp (t)sp (t − τ )
(36)


p
X 2
Ãqp E {sp (t)sp (t − τ )}
(37)
p

(38)

If we assume, that the fundamental frequency f0,p is different
for each speaker, we can estimate the parameters from the mixtures for each source in parallel.
In each scenario a non-parametric estimation is the basis for theparameter estimation. As a general least square cost function for all
scenarios we propose
X X
J (γ, β, f 0 ) =q(rq(τ )−
(39)
τ

X
p

Ã2q,p

M
X

!

!2

γmp,p cos(m · f0,p · τ ) exp(−βp τ ) 

m=1

6. EXPERIMENTAL RESULTS
We performed experiments with three speech signals (two male, one
female) of 1.62 seconds (8kHz sampling rate). We generated two
instantaneous mixtures with the mixing matrix


1
1
1
Ã = 1.3 −0.9 0.8
(40)

t=0

5.2 Parametric estimation

rxq (τ )Ã2qp rsp (τ )p

5.1 Non-parametric estimation
Non-parametric estimation has the advantage that no assumption
about the structure is necessary. For our purpose with only a limited
number T of samples of the signal in question available, we choose
a biased estimation technique [10]. It guarantees that the estimatedautocorrelation is positive semi-definite. For a signal frame s(t),0 ≤ t ≤ T − 1 the biased autocorrelation estimate is given byX

and assumed that the mixing matrix is available. The signals were
processed in frames of T = 256 samples and shifted by 64 samples.
No noise was added to the mixtures, but the noise variance was set
to σ 2 = 10−4 . This led to a good influence of the autocorrelation
matrix R. For approximations of the autocorrelation we set M = 2.
The sources had equal variance of 0.0037. We compared thefollowing approaches:
Separating the mixtures by `1 -norm minimization in
the frequency domain, which is based on a univariate
speech model [17].
L1G Using the results from L1 to estimate the autocorrelation. This estimate was used to further enhance theseparated signals by the Gibbs sampler (24)-(29) (K =1000 iterations).
OG
Estimating the signals with the Gibbs sampler (24)-(29)(K = 10 iterations), whereby the autocorrelation was
estimated from the original signals (33).
OI
With the autocorrelation estimated from the original
signals (33), the signals were obtained with the modified generalized inverse (31).
OAI
Estimating the signals with the modified generalized
inverse (31), whereby each autocorrelation was approximated by two harmonic, exponentially weighted
cosines based on the original signals (34).
OSI
Estimating the signals with the modified generalized inverse (31), whereby the autocorrelation matrix was a
diagonal matrix filled with the first lag of the autocorrelation of the original signals (33).
MAI Estimating the signals with the modified generalized inverse (31), whereby the autocorrelation was approximated directly from the mixtures (34).
L1

For the evaluation of the separation results we used the signalto-distortion ratio (SDR), signal-to-interference ratio (SIR) andsignal-to-artifact ratio (SAR) as defined in [6]. The averaged resultsof the different approaches are shown in Table 1.
7. CONCLUSION
The experimental results in Table 1 together with subjective evaluation suggest that knowing the autocorrelation together with themixing matrix is sufficient to perform high quality underdeterminedsource separation in noiseless environments (OG and OI comparedto L1G, OSI and MAI). In other words, the problem of underdetermined source separation can be reduced to the estimation of the underlying autocorrelations once the mixing matrix is available. The

Method
L1
L1G
OG
OI
OAI
OSI
MAI

SDR
9.28
5.29
12.3
18.21
10.42
11.11
3.45

SIR
14.84
17.65
35.61
34.33
20.75
18.08
13.63

SAR
11.79
7.25
12.56
18.51
11.40
12.89
5.83

∝

p

=

Table 1: Experimental results

A. CONDITIONAL PDFS FOR GIBBS SAMPLER• Conditional PDF of speech signalp(S|θ −S , X) ∝ p(S|R) · p(X|A, S, σ 2 )

∝

∝

(41)

2

N (S|R) · N (X|AS, σ )
1 T −1
exp − S R S ·
2


1
I
exp − (X − AS)T 2 (X − AS)
2
σ






 1 T
1 T
−1
S
exp 
−
R
+
A
A
S−
 2
σ2

|
{z
}
:=S T Σ−1
S S


−S T



1 T
1 T 
1 T

A
X
−
X
AS
+
X
X

2
2
2
σ
σ
σ

| {z }
:=Σ−1
S µS

=

N (S|µS , ΣS )with ΣS and µS given by (27) and (28), respectively.
• Conditional PDF of autocorrelationp(R|θ −R , X) ∝ p(S|R) · p(R|α, Σ)YY∝N (S p |0, Rp ) ·IW(Rp |αp , Σp )
∝

p

p

Y

N (S p |0, Rp ) · IW(Rp |αp , Σp )

p

∝

Y
p

− 21

|Rp |



1 T −1
exp − S p Rp S p ·
2


 1 

1 
T
−1 −1
exp − tr R−1p S p S p − tr Rp Σp
2
2


−1 
Y
−1
IW Rp |αp + 1, S p S T
+
Σ
p
p
p

proposed parameterization of the autocorrelation did not yet leadto satisfying results (OAI, MAI). Therefore, either better autocorrelation models or improved non-parametric estimation methods arenecessary.
The experimental results obtained by the MMSE approach OI
resemble the results of the Gibbs sampler OG. This suggests thatthe exact prior distribution of the sources plays only a minor role
as long as the autocorrelation is taken into account. Therefore, thecomparatively easy and fast solution given by (31) might provide abasis for algorithms that yield better autocorrelation estimates andeventually underdetermined source separation.
Possible options for improving the autocorrelation estimate include linear prediction based models and independent component
analysis (ICA) techniques. Linear prediction has already provenin speech coding, that it allows very compact yet high quality descriptions of speech signals. ICA exploits the statistical independence between the sources and might be helpful for estimating the
involved parameters.

∝



1
1 
−1
|Rp |− 2 (αp +T +1) exp − tr R−1
p Σp2
Y− 21 (αp +T +2)
|Rp |
·

(42)

References[1] J.M. Bernardo and A.F.M. Smith. Bayesian theory. Wiley
series in probability and statistics. Wiley, 2000.
[2] A. Blin, S. Araki, and S. Makino. Underdetermined blindseparation of convolutive mixtures of speech using timefrequency mask and mixing matrix estimation. IEICE Trans.
Fundamentals, E88-A(7):1693–1700, 2005.
[3] H. Brehm and W. Stammler. Description and generation of
spherically invariant speech-model signals. Signal Processing,
12:119–141, 1987.
[4] A. Doucet and X. Wang. Monte Carlo methods for signal processing: a review in the statistical signal processing context.
IEEE Signal Processing Magazine, 22(6):152–170, Nov 2005.
[5] C. Févotte and S. J. Godsill. A bayesian approach for blindseparation of sparse sources. Technical report, Cambridge
University Engineering Dept., January 2005.
[6] C. Févotte, R. Gribonval, and E. Vincent. BSS EVAL toolbox
user guide – Revision 2.0. Technical Report 1706, IRISA,
April 2005.
[7] C. Fraley and A.E. Raftery. Bayesian regularization for normal mixture estimation and model-based clustering. Technical Report 486, University of Washington, Statistics, August
2005.
[8] A. Gelman, J.B. Carlin, H.S. Stern, and D.B. Rubin. Bayesiandata analysis. Chapman & Hall, 1995.
[9] D.J.C. MacKay. Information Theory, Inference, and Learning
Algorithms. Cambridge University Press, 7 edition, Aug 2004.
[10] D.G. Manolakis, V.K. Ingle, and S.M. Kogon. Statistical andadaptive signal processing. McGraw-Hill, 2000.
[11] J.M. Mendel. Tutorial on higher-order statistics (spectra) insignal processing and system theory: theoretical results andsome applications. Proceedings of the IEEE, 79(3):278–305,
March 1991.
[12] A. Papoulis and S.U. Pillai. Probability, random variables,and stochastic processes. McGraw-Hill, 4th edition, 2002.
[13] D.B. Rowe.
Multivariate Bayesian statistics: models
for source separation and signal unmixing. Chapman &Hall/CRC, 2003.
[14] F.J. Theis. Mathematics in independent component analysis.
PhD thesis, University of Regensburg, 2002.
[15] L. Vielva, I. Santamaria, C. Pantaleon, J. Ibanez, and D. Erdogmus. Estimation of the mixing matrix for underdeterminedblind source separation using spectral estimation techniques.
In Proc. EUSIPCO 2002, volume 1, pages 557–560, Sep 2002.
[16] S. Winter, H. Sawada, S. Araki, and S. Makino. Overcomplete
BSS for convolutive mixtures based on hierarchical clustering.
In Proc. ICA 2004, pages 652–660, Sept. 2004.
[17] S. Winter, H. Sawada, and S. Makino. On real and complex
valued L1-norm minimization for overcomplete blind source
separation. In 2005 IEEE Workshop on Applications of Signal
Processingto Audio and Acoustics (WASPAA), pages 86–89,
New Paltz, NY, USA, 2005.
[18] M. Zibulevsky and B.A. Pearlmutter. Blind source separationby sparse decomposition. Neural Computations, 13(4):863 –
882, 2001.