International Workshop on Acoustic Signal Enhancement 2012, 4-6 September 2012, Aachen

BLIND SEPARATION OF INFINITELY MANY SPARSE SOURCES
Hirokazu Kameoka1,2 , Misa Sato3 , Takuma Ono1 , Nobutaka Ono4 , Shigeki Sagayama1
12Graduate School of Information Science and Technology, The University of TokyoNTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation
3
School of Engineering, The University of Tokyo4
Principles of Informatics Research Revision, National Institute of Informatics
ABSTRACTThis paper deals with the problem of underdetermined blindsource separation (BSS) where the number of sources is unknown. We propose a BSS approach that simultaneously
estimates the number of sources, separates the sources basedon the sparseness of speech, and performs permutation alignment. We conﬁrmed experimentally that reasonably good
separation was obtained with the present method without
specifying the number of sources.
Index Terms— Underdetermined blind source separation, sparseness, Dirichlet process, variational inference
1. INTRODUCTION
Blind Source Separation (BSS) is a technique for separatingout individual source signals from microphone inputs when
the transfer characteristics between sources and microphones
are unknown. BSS is potentially useful for the developmentof such applications as hands-free teleconference systems andautomatic meeting transcription systems. In meeting situations, for instance, it is likely that the number of participants(speakers) will change during the meeting or a loud, unexpected noise such as a door slamming will occur in the room.
Thus, it is often difﬁcult to pre-specify the exact number ofall possible sources present in real environments. Many conventional BSS algorithms are designed to use the number ofsources as the input when performing separation, and most ofthese algorithms do not work well if the assumed and actual
numbers of sources are not the same. This paper presents a
novel BSS algorithm that allows the number of sources to beinferred along with the separation.
To estimate the unknown mixing matrix and source signals solely from observed signals, we must make some assumptions about the sources, and formulate an appropriate
optimization problem based on criteria designed according tothose assumptions. For example, if the observed signals outnumber the sources, we can employ independent component
analysis (ICA) [1], which estimates the separation matrix (theinverse of the mixing matrix) such that the independence ofthe source estimates is maximized. However, to achieve a
BSS algorithm that works without assuming the number ofsources, we shall always consider an underdetermined casewhere there are fewer observations than sources. In an underdetermined case, there are an inﬁnite number of solutions forsource signals even if the mixing process is known. The independence assumption is too weak to allow us to determine
a unique solution and so directly applying ICA will not workin this case. In underdetermined situations, we typically needa stronger assumption than independence.
One successful approach for underdetermined BSS involves utilizing the sparse nature of speech [2–6], which relies on the fact that the time-frequency components of speech
are near zero across most of the time-frequency bins. Since
the time-frequency components of speech rarely overlap eachother even when multiple speakers are speaking simultaneously, the main focus of this approach is how to design atime-frequency mask with which we can extract only thecomponents of target speech from the mixture.
The signals observed at each microphone can be modeled as a convolutive mixture of source signals. To exploit
the sparse nature of speech, we must convert it to a timefrequency representation. If we assume the use of a shorttime Fourier transform (STFT) to obtain a time-frequency
representation with a frame length sufﬁciently longer thanthe length of the impulse response from a source to a microphone, an observed signal can be approximated fairly wellby an instantaneous mixture in the frequency domain. BSS
based on this observation model is called frequency domain
BSS. While frequency domain BSS allows for a fast implementation compared with BSS that uses a time domain convolutive mixture model, it requires us to solve an additional
problem called the permutation alignment problem. That is,
we must group together the separated components of different
frequency bins that are considered to originate from the same
source in order to construct a separated signal.
Motivated by the above, this paper proposes a novel BSS
approach that simultaneously performs (1) an estimation of
the number of sources, (2) source separation based on thesparseness of speech, and (3) permutation alignment.
2. MIXING MODELWe ﬁrst consider a situation where K source signals are captured by M microphones. Here, let ym (ω, t) ∈ C be theshort-time Fourier transform (STFT) component observed at
the m-th microphone, and sk (ω, t) ∈ C be the STFT component of the k-th source. 1 ≤ ω ≤ Ω and 1 ≤ t ≤ T are
the frequency and time indices, respectively. If we assume
that the length of the impulse response from a source to a mi-crophone is sufﬁciently shorter than the frame length of the
STFT, the observed signal can be approximated fairly well by
an instantaneous mixture in the frequency domain:y(ω, t) =K∑

ak (ω)sk (ω, t) + n(ω, t),

(1)

k=1

where y(ω, t) = (y1 (ω, t), . . . , yM (ω, t))T . ak (ω) =(ak,1 (ω), . . . , ak,M (ω))T is the frequency array response
for source k, which is assumed to be time-invariant. n(ω, t)
is assumed to comprise all kinds of components that cannot be expressed by the instantaneous mixture representation(e.g., background noise and reverberant components).
We now exploit the sparseness of speech and assume that
only one source is active in each time-frequency bin. By using zω,t ∈ {1, . . . , K} to denote the (unknown) active source
index at time-frequency bin (ω, t), (1) can be rewritten asy(ω, t) = azω,t (ω)s(ω, t) + n(ω, t).

3. GENERATIVE MODEL3.1. Generative process of observed signals
Here we describe the generative process of an observed signal
on the basis of (2). Let us assume that the noise componentnω,t follows a complex normal distribution with mean 0 and(n)covariance Σω . Then, from (2), y ω,t is also normally distributed such thaty ω,t |a1:K,ω , sω,t , zω,t ∼ NC (y ω,t ; azω,t ,ω sω,t , Σ(n)ω ), (3)conditioned on a1:K,ω = {a1,ω , . . . , aK,ω }, sω,t and zω,t ,where NC (x; µ, Σ) ∝ exp(−(x − µ)H Σ−1 (x − µ)).
3.2. Mixture of inﬁnite sparse sources
We do not usually obtain any information about which sourceis active at each time-frequency bin. Thus, we regard zω,t as alatent variable and consider describing its generative process.
If the number of sources is K, it would be natural to assume
that the probability of choosing an index zω,t from the setof source indices, {1, . . . , K}, can be described as a discrete
distribution(4)

where Discrete(x; y) = yx (with y = (y1 , . . . , yK ) and∑K
k=1 yk = 1). The k-th element of π deﬁnes how likely
the source index k is to be chosen. Since we do not also have
any information about π, we consider describing its generative process using a “symmetric” distribution. For the convenience of the following analysis, we assume that π has been
generated from a symmetric Dirichlet distribution
π ∼ Dirichlet(π; α0 /K, . . . , α0 /K),

π ∼ GEM(π; α0 ),zω,t |π ∼ Discrete(zω,t ; π),y ω,t |a1:∞,ω sω,t , zω,t ∼ NC (y ω,t ; azω,t ,ω sω,t , Σ(n)
ω ).

(5)

(6)
(7)
(8)

GEM(π; α0 ) is called the stick-breaking process [7], which is
known as a constructive deﬁnition of the Dirichlet process [8]and is given byvk ∼ Beta(vk ; 1, α0 ),(2)Notice that the superscript k is dropped from sk (ω, t) in (2) as
it is no longer necessary since we are assuming sk (ω, t) = 0
for k ̸= zω,t . Namely, s(ω, t) signiﬁes the component of anactive source at time-frequency bin (ω, t). For convenienceof notation, we hereafter use subscripts to indicate ω and t.

zω,t |π ∼ Discrete(zω,t ; π),∏K
where Dirichlet(y; z1 , . . . , zK ) ∝ k=1 ykzk −1 . The shapeof the Dirichlet distribution is governed by a concentration
hyperparameter α0 .
Thus far, we have considered the case of a ﬁnite number
K of sources. It can be shown that the limit of the above
generative processes (3), (4) and (5) as K goes to inﬁnity canbe described explicitly as

πk = vkk−1∏(1 − vl ),

(9)
(10)

l=1

where Beta(y; z1 , z2 ) ∝ y z1 −1 (1−y)z2 −1 . π ∼ GEM(π; α0 )produces exponentially decaying weights in expectation. Thismeans that the source with a larger index will be less likely
to be active and thus simple models with fewer sources arefavored, given observations.
3.3. Mixture of Direction-of-Arrivals (DOAs)
We describe here the generative process of the frequency response ak,ω of the mixing system.
So far we have treated ak,ω as an independent parameter
across ω. If the index k indicates an identical source acrossω, ak,ω will have a certain structure that can be describedusing the property of acoustic wave propagation. We thus
expect that incorporating an appropriate constraint into ak,ωwould help solve the permutation alignment problem throughparameter inference. If each source is assumed to be located
far from the microphones so that the signal can be treated approximately as a plane wave, the interchannel time difference
between the microphones depends only on the direction ofarrival (DOA) of the source. Since the time delay between
two microphones corresponds to the phase difference of thefrequency response of the microphone array, the complex array response can be expressed explicitly by using the DOAsof the source. Speciﬁcally, with M = 2 microphones, thecomplex array response for a source at direction θ such that0 ≤ θ < 2π is deﬁned as a function of ω depending on θ[]1ℓθ,ω = ȷωD cos θ/C ,(11)ewhere ȷ is the imaginary unit, D [m] is the distance between
the two microphones, and C [m/s] is the speed of sound. If
the DOA θk of source k is known, the frequency responseak,ω should be equal to ℓθk ,ω . However, due to such factors
as the plane wave assumption and the narrowband instantaneous mixture approximation, the actual frequency responseak,ω may diverge from the “ideal frequency response” ℓθk ,ω

A = {ak,ω }k,ω : complex array response for source k,
S = {sω,t }ω,t
: component of active source at (ω, t),
Z = {zω,t }ω,t
: index of active source at (ω, t),
V = {vk }k: stick breaking proportion,
C = {ck }k: index of direction for source k,ρ = (ρ1 , . . . , ρN ) : mixture weight for each DOA,
which we denote by Θ. Our goal is to compute the posterior p(Θ|Y ) where Y = {y ω,t }ω,t is the set consisting ofobserved multichannel signals. By using the conditional distributions deﬁned in 3.2 and 3.3, we can write the joint distribution p(Y, Θ) as
Fig. 1. Plate notation of the present generative model.
to some extent. One way to simplify the process of this kindof deviation is to assume a probability distribution on ak,ωwith the expected value of ℓθk ,ω . Here, we assume for convenience that ak,ω is generated from a complex normal distribution with mean ℓθk ,ω . Note that we do not usually obtain
any information about which direction each source is emanating from. Thus, we regard the DOA of each source asa latent variable and further consider describing its generative process. As explained in detail below, the entire generative process of ak,ω can then be described as a “mixture of
DOAs”.
Let us now introduce a discrete set of N possible directions, ϑ1 , . . . , ϑN , which are all assumed to be constants.
For instance, consider deﬁning ϑn as ϑn = (n − 1)π/N ,(n = 1, . . . , N ) (dividing π in N equal angles). We assume
that each source is emanating from one of these directions.
First, we consider the generative process of the DOA θk ofsource k. For each source k, an index ck of direction is drawn
according to a discrete distribution ρ = (ρ1 , . . . , ρN )ck |ρ ∼ Discrete(ck ; ρ).

(12)

By using ck , θk is then given asθk = ϑck .

(13)

As with 3.2, we assume that ρ has been generated from asymmetric Dirichlet distributionρ ∼ Dirichlet(ρ; β0 /M, . . . , β0 /M ).

(14)

As mentioned above, the frequency response ak,ω is assumed
to be generated from a complex normal distribution, given ck ,ak,ω |ck ∼ NC (ℓθk ,ω , Σ(a)
ω ),

(15)

(a)

where Σω denotes the covariance of the complex normal distribution. The entire generative model is described in Fig. 1.
4. APPROXIMATE POSTERIOR INFERENCE
In this section, we describe an approximate posterior inference algorithm for our generative model based on variational
inference. The random variables of interest in our model arep(Y, A, S, Z, V, C, ρ)= p(Y |A, S, Z)p(Z|V )p(V )p(A|C)p(C|ρ)p(ρ), (16)but to obtain the exact posterior p(Θ|Y ), we must computep(Y ), which involves many intractable integrals.
We can express this posterior variationally as the solution
to an optimization problem:
argmin KL(q(Θ)∥p(Θ|Y )),

(17)

q∈Q

where KL(·∥·) denotes the Kullback-Leibler (KL) divergence
between its two arguments. Indeed, if we let Q be the family of all distributions over Θ, the solution to the optimization
problem is the exact posterior p(Θ|Y ), since KL divergenceis minimized when its two arguments are exactly equal. Of
course, solving this optimization problem is just as intractable
as directly computing the posterior. Although it may appear
that no progress has been made, having a variational formulation allows us to consider tractable choices of Q in order toobtain principled approximate solutions.
For our model, we deﬁne the set of approximate distributions Q to be those that factor as follows:{}
Q = q : q(A)q(S)q(Z)q(V )q(C)q(ρ) .
(18)
To deﬁne q(A), q(V ) and q(C), we need to construct distributions on the inﬁnite sets {v1 , v2 , . . .}, {a1,ω , a2,ω , . . .}and {c1 , c2 , . . .}. For this approach to be tractable, we truncate the variational distribution at some value K ∗ by settingq(vK ∗ = 1) = 1. The mixture proportions πk for k > K ∗will thus be zero, and we can ignore ak,ω and ck for k > K ∗ .
In practice, we set K ∗ at a sufﬁciently large integer. It is important to emphasize that truncating the variational distribution does not mean that the true posterior itself is truncated.
As the truncation level K ∗ becomes larger, the approximations to the true posterior become more accurate.
We now present an algorithm for solving the optimization
problem described in (17) and (18). Unfortunately, the optimization problem is non-convex, and it is intractable to ﬁnd
the global optimum. However, we can use a simple coordinate
ascent algorithm to ﬁnd a local optimum. The algorithm optimizes one factor in the mean-ﬁeld approximation of the posterior at a time while ﬁxing all other factors. The mean-ﬁeld
update equations for the variational distributions are given in
the following form:∏q̂(A) =NC (ak,ω ; mk,ω , Γk,ω ),
(19)
k,ω

the present methods with various K and K ∗ settings. As expected, the performance of the present method improves withincreasing K ∗ , while that of Sawada’s method deteriorates
signiﬁcantly when the assumed number of sources departs
from the actual number.
6. CONCLUSION

Fig. 2. SDRs obtained with conventional and present methods
with different K and K ∗ settings.
q̂(S) =∏

7. REFERENCES
2
NC (sω,t ; µω,t , σω,t
),

(20)

Discrete(zω,t ; ϕω,t ),(21)

Beta(vk ; γk,0 , γk,1 ),(22)

Discrete(ck ; ψ k ),

(23)

ω,t

q̂(Z) =∏
ω,tq̂(V ) =∏kq̂(C) =∏kq̂(ρ) =∏

This paper aimed at developing a BSS algorithm that works
well even when the number of sources is unknown and proposed a novel BSS approach that simultaneously performs anestimation of the number of sources, source separation basedon the sparseness of speech, and permutation alignment. We
conﬁrmed experimentally that reasonably good separations
were obtained with the present method without specifying thenumber of sources.

Dirichlet(ρ; ζk,1 , . . . , ζk,M ).

(24)

k

5. EXPERIMENT
We evaluated the performance of the proposed method interms of source separation ability.
We used a stereo speech signal with a sampling rate of
16kHz as a test signal, which we obtained by mixing threespeech signals [9] (two female and one male speakers) using a measured room impulse response [10] (in which the
distance between the microphones was 5 cm and the reverberation time was 0 ms). The three sources were spaced 30
degrees apart. To compute the STFT components of the observed signal, the STFT frame length was set at 64 ms and a
Hamming window was used with an overlap length of 16ms.
(a)(n)Σω and Σω were set respectively at I and 10−1.5 × I. M
was set at 180. All the variational parameters were initialized randomly. The variational inference algorithm was run
for 100 iterations. After convergence, each separated signal was obtained by multiplying µω,t by ϕk,ω,t . We chose
Sawada’s method described in [11] as a comparison. In thismethod, the number of sources must be speciﬁed manually.
The following results report the performance in terms of the
Signal-to-Distortion Ratio (SDR) [12]. The SDR is expressedin decibels (dB), and a higher SDR indicates superior quality.
The present method was tested with various settings of thetruncation level K ∗ . As for Sawada’s method, it was tested
with various settings of the assumed number K of sources.
Fig. 2 shows the average SDRs obtained with Sawada’s and[1] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, John Wiley & Sons, 2001.
[2] Ö. Yılmaz and S. Rickard, “Blind separation of speech mixtures via time-frequency masking,” IEEE Trans. Signal Process., vol. 52, no. 7, pp. 1830–1847, 2004.
[3] M. I. Mandel, D. P. W. Ellis, and T. Jebara, “An EM algorithm
for localizing multiple sound sources in reverberant environments,” in Adv. Neural Info. Process. Syst., 2006, pp. 953–960.
[4] S. Araki, H. Sawada, R. Mukai, and S. Makino, “Underdetermined blind sparse source separation for arbitrarily arranged
multiple sensors,” Signal Process., vol. 87, no. 8, pp. 1833–
1847, 2007.
[5] Y. Mori, H. Saruwatari, T. Takatani, S. Ukai, K. Shikano,
T. Hiekata, and T. Morita, “Real-time implementation of twostage blind source separation combining SIMO-ICA and binary masking,” in Proc. IWAENC ’05, 2005, pp. 229–232.
[6] Y. Izumi, N. Ono, and S. Sagayama, “Sparseness-based 2ch
BSS using the EM algorithm in reverberant environment,” inProc. WASPAA ’07, 2007, pp. 147–150.
[7] J. Sethuraman, “A constructive deﬁnition of Dirichlet priors,”
Statistica Sinica, vol. 4, pp. 639–650, 1994.
[8] T. S. Ferguson, “A Bayesian analysis of some nonparametricproblems,” Annals of Statistics, vol. 1, no. 2, pp. 209–230,
1973.
[9] A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri,
H. Kuwabara, and K. Shikano, “ATR Japanese speech
database as a tool of speech recognition and synthesis,” Trans.
Speech Communication, pp. 357–363, 1990.
[10] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada, “Acoustical sound database in real environments forsound scene understanding and hands-free speech recognition,” in Proc. LREC ’00, 2000, pp. 965–968.
[11] H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment,” IEEE Trans. ASLP, vol.
19, no. 3, pp. 516–527, 2010.
[12] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans.
ASLP, pp. 1462–1469, 2006.