IEICE TRANS. FUNDAMENTALS, VOL.E96–A, NO.10 OCTOBER 2013

1928

PAPER

Special Section on Sparsity-aware Signal Processing

Bayesian Nonparametric Approach to Blind Separation of Infinitely
Many Sparse Sources
Hirokazu KAMEOKA†,††a) , Member, Misa SATO†††∗ , Takuma ONO† , Nobutaka ONO†††† ,and Shigeki SAGAYAMA† , Nonmembers

SUMMARY
This paper deals with the problem of underdeterminedblind source separation (BSS) where the number of sources is unknown.
We propose a BSS approach that simultaneously estimates the number ofsources, separates the sources based on the sparseness of speech, estimates
the direction of arrival of each source, and performs permutation alignment.
We confirmed experimentally that reasonably good separation was obtainedwith the present method without specifying the number of sources.
key words: underdetermined blind signal separation, speech, sparseness, permutation alignment, Bayesian nonparametrics, direction of arrival, Dirichlet process, stick-breaking construction, variational inference

1.

Introduction

Blind Source Separation (BSS) is a technique for separating out individual source signals from microphone inputs
when the transfer characteristics between sources and microphones are unknown. BSS is potentially useful for thedevelopment of such applications as hands-free teleconference systems and automatic 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 diﬃcult to pre-specify
the exact number of all possible sources present in real environments. Many conventional BSS algorithms are designed
to use the number of sources as the input when performing
separation, and most of these 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 thenumber of sources to be inferred 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
Manuscript received January 22, 2013.
Manuscript revised May 23, 2013.
†
The authors are with the Graduate School of Information Science and Technology, The University of Tokyo, Tokyo, 113-8656
Japan.
††
The author is with NTT Communication Science Laboratories, NTT Corporation, Atsugi-shi, 243-0198 Japan.
†††
The author is with School of Engineering, The University ofTokyo, Tokyo, 113-8656 Japan.
††††
The author is with National Institute of Informatics, Tokyo,
101-8430 Japan.
∗
Presently, with Graduate School of Engineering, The University of Tokyo.
a) E-mail: kameoka@hil.t.u-tokyo.ac.jp
DOI: 10.1587/transfun.E96.A.1928

to those assumptions. For example, if the observed signals
outnumber the sources, we can employ independent component analysis (ICA) [1], which estimates the separation
matrix (the inverse of the mixing matrix) such that the independence of the source estimates is maximized. However,to achieve a BSS algorithm that works without assuming thenumber of sources, we must always consider an underdetermined case where there are fewer observations than sources.
In an underdetermined case, there are an infinite number ofsolutions for source signals even if the mixing process isknown. The independence assumption is too weak to allow
us to determine a unique solution and so directly applying
ICA will not work in this case. In underdetermined situations, we typically need a stronger assumption than independence.
One successful approach for underdetermined BSS involves utilizing the sparse nature of speech [2]–[7], which
relies on the fact that the time-frequency components ofspeech are near zero across most of the time-frequency grid
points. Since the time-frequency components of speechrarely overlap even when multiple speakers are speaking simultaneously, the main focus of this approach is how to design a time-frequency mask with which we can extract only
the components 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 suﬃ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 diﬀerent frequency bins that are considered to originate from thesame source in order to construct a separated signal. Somemethods such as [7] are designed to perform frequency-binwise source separation followed by permutation alignment.
However, permutation aligment and source separation problems should be solved in a cooperative manner since theclues used for permutation alignment can also be helpful

c 2013 The Institute of Electronics, Information and Communication Engineers
Copyright KAMEOKA et al.: BAYESIAN NONPARAMETRIC APPROACH TO BLIND SEPARATION OF INFINITELY MANY SPARSE SOURCES

1929

for source separation. Thus, we consider it important to develop a method that can simultaneously perform frequencybin-wise source separation and permutation alignment.
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 basedon the sparseness of speech, and (3) permutation alignment.
2.

Mixing Model

We first consider a situation where N source signals are captured by M microphones. Here, let ym (ωk , tl ) ∈ C be theshort-time Fourier transform (STFT) component observedat the m-th microphone, and si (ωk , tl ) ∈ C be the STFT
component of the i-th source. ωk is the angular frequency
corresponding to the k-th frequency bin and tl is the timeof the l-th frame, respectively. If we assume that the lengthof the impulse response from a source to a microphone issuﬃciently shorter than the frame length of the STFT, theobserved signal can be approximated fairly well by an instantaneous mixture in the frequency domain:y(ωk , tl ) =Nai (ωk )si (ωk , tl ) + n(ωk , tl ),

(1)

i=1

where y(ωk , tl ) = (y1 (ωk , tl ), . . . , y M (ωk , tl ))T . ai (ωk ) =(ai,1 (ωk ), . . . , ai,M (ωk ))T is the frequency array response forsource i, which is assumed to be time-invariant. n(ωk , tl ) 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 at each time-frequency point,
as with [2]–[7]. By using zk,l ∈ {1, . . . , N} to denote the (unknown) active source index at time-frequency point (ωk , tl ),(1) can be rewritten asy(ωk , tl ) = azk,l (ωk )s(ωk , tl ) + n(ωk , tl ).

(2)

Notice that the superscript i is dropped from si (ωk , tl ) in(2) as it is no longer necessary since we are assumingsi (ωk , tl ) = 0 for i  zk,l . Namely, s(ωk , tl ) signifies the component of an active source at time-frequency point (ωk , tl ).
For convenience of notation, we hereafter use subscripts k
and l to indicate ωk and tl .
3.

conditioned on a1:N,k = {a1,k , . . . , aN,k }, sk,l and zk,l , whereNC (x; μ, Σ) ∝ exp(−(x − μ)H Σ−1 (x − μ)). Using this model,
Izumi et al. formulated an Expectation-Maximization (EM)algorithm for the maximum likelihood estimation of the unknown parameters by treating zk,l as a latent variable [6], [7].
To flexibly incorporate hierarchical prior distributions overthe unknown parameters as well as the number of sources,
we employ a Bayesian approach.
3.2 Mixture of Infinitely Many Sparse Sources
We do not usually obtain any information about whichsource is active at each time-frequency point. Thus, we
regard zk,l as a latent variable and consider describing its
generative process. If the number of sources is N, it would
be natural to assume that the probability of choosing an index zk,l from the set of source indices, {1, . . . , N}, can be
described as a categorical distribution
zk,l |π ∼ Categorical(zk,l ; π),

3.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 component nk,l follows a complex normal distribution with mean
0 and covariance Σk(n) . Then, from (2), yk,l is also normally
distributed such that
(3)

(4)

where
Categorical(x; y) = y x (with y = (y1 , . . . , yN ) andNy= 1). The i-th element of π defines how likely theii=1
source index i is to be chosen. Since we also have no information about π, we consider describing its generative process using a “symmetric” distribution. For the convenienceof the following analysis, we assume that π has been generated from a symmetric Dirichlet distribution
π ∼ Dirichlet(π; α0 /N, . . . , α0 /N),

Generative Modelyk,l |a1:N,k , sk,l , zk,l ∼ NC (yk,l ; azk,l ,k sk,l , Σk(n) ),Fig. 1 Spectrograms of speech signals uttered by three speakers, each ofwhich is indicated by a diﬀerent color. Typically, the time-frequency components of speech rarely overlap even when multiple speakers are speaking
simultaneously.

(5)

N zi −1
where Dirichlet(y; z1 , . . . , zN ) ∝ i=1
yi . The shape of the
Dirichlet distribution is governed by a concentration hyperparameter α0 .
Thus far, we have considered the case of a finite number N of sources. It can be shown that with a size-biased
permutation of (5) followed by taking the limit N → ∞, theabove generative processes (3), (4) and (5) approachπ ∼ GEM(π; α0 ),zk,l |π ∼ Categorical(zk,l ; π),

(6)
(7)

IEICE TRANS. FUNDAMENTALS, VOL.E96–A, NO.10 OCTOBER 2013

1930of the source. Specifically, with M = 2 microphones, thecomplex array response for a source at direction θ such that0 ≤ θ < 2π is defined as a function of ω depending on θ


1
(11)
h(θ, ω) = jωB cos θ/C ,eFig. 2 Illustration of the stick-breaking construction. The constructionof π = (π1 , π2 , . . . , π∞ ) can be understood metaphorically as follows. Starting with a stick of length 1, we break it at v1 , assigning π1 as the length of
the stick we just broke oﬀ. We now recursively break the other portion to
obtain π2 , π3 and so forth.

yk,l |a1:∞,k sk,l , zk,l ∼ NC (yk,l ; azk,l ,k sk,l , Σk(n) ).

(8)GEM(π; α0 ) is called the stick-breaking process [8], where
the letters stand for Griﬃths, Engen and McCloskey. Fori = 1, 2, . . . ∞, the i-th mixture weight πi is generated via
vi ∼ Beta(vi ; 1, α0 ),πi = vii−1
(1 − v j ),

(9)
(10)

j=1

where Beta(y; z1 , z2 ) ∝ yz1 −1 (1 − y)z2 −1 . This is known asa constructive definition of the Dirichlet process [9], which
can be thought of as an infinite dimensional Dirichlet distribution. The construction of π = (π1 , π2 , . . . , π∞ ) can beunderstood metaphorically as follows. Starting with a stickof length 1, we break it at v1 , assigning π1 as the lengthof the piece of stick that we have just broken oﬀ. Now
recursively break the other portion to obtain π2 , π3 and soforth (see Fig. 2). π ∼ GEM(π; α0 ) thus produces exponentially decaying weights in expectation. This means that thesource with a larger index will be less likely to be active andthus simple models with fewer sources are favored, given
observed signals.
3.3 Mixture of Direction-of-Arrivals (DOAs)
Here we describe the generative process of the frequency
response ai,k of the mixing system.
So far we have treated ai,k as an independent parameter across k. If the index i indicates an identical source
across ω1 , . . . , ωK , ai,k will have a certain structure that canbe described using the property of acoustic wave propagation. We thus expect that the incorporation of an appropriate
constraint into ai,k would help solve both the permutationalignment problem and the frequency-wise source separation problem simultaneously through parameter 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 diﬀerence between
the microphones depends only on the direction of arrival(DOA) of the source. Since the time delay between two
microphones corresponds to the phase diﬀerence of the frequency response of the microphone array, the complex array response can be expressed explicitly by using the DOAs

where j is the imaginary unit, B [m] is the distance between
the two microphones, and C [m/s] is the speed of sound. If
the DOA θi of source i is known, the frequency response
ai,k should be equal to h(θi , ωk ). However, due to such factors as the plane wave assumption and the narrowband instantaneous mixture approximation, the actual frequency response ai,k may diverge from the “ideal frequency response”
h(θi , ωk ) to some extent. One way to simplify the process of
this kind of deviation is to assume a probability distribution
on ai,k with the expected value of h(θi , ωk ). Here, we assume for convenience that ai,k is generated from a complex
normal distribution with mean h(θi , ωk ). Note that we do
not usually obtain any information about the direction fromwhich each source is emanating. Thus, we regard the DOA
of each source as a latent variable and further consider describing its generative process. As explained in detail below,
the entire generative process of ai,k can then be described asa “mixture of DOAs”.
Let us firstly introduce a discrete set of D possible directions, ϑ1 , . . . , ϑD , which are all assumed to be constants.
For instance, consider defining ϑd as ϑd = (d − 1)π/D,(d = 1, . . . , D), which means dividing π into D equal angles.
We then assume that each source is emanating from one ofthese directions. First, we consider the generative processof the DOA θi of source i. For each source i, an index ciof direction is drawn according to a categorical distributionρ = (ρ1 , . . . , ρD )ci |ρ ∼ Categorical(ci ; ρ).

(12)

By using ci , θi is then given as
θi = ϑci .

(13)

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

(14)

As mentioned above, the frequency response ai,k is assumed
to be generated from a complex normal distribution withmean h(ϑci , ωk ), given ci ,ai,k |ci ∼ NC (ai,k ; h(ϑci , ωk ), Σ(a)
k ),

(15)

where Σk(a) denotes the covariance of the complex normal
distribution, which is assumed to be a constant. An illustration of the generative process of ai,k described above isshown in Fig. 3.
Overall, the entire generative model is described inplate notation in Fig. 4.

KAMEOKA et al.: BAYESIAN NONPARAMETRIC APPROACH TO BLIND SEPARATION OF INFINITELY MANY SPARSE SOURCES

1931

argmin KL(q(Θ) p(Θ|Y)),

(17)

q∈Q

where KL(· ·) denotes the Kullback-Leibler (KL) divergence between its two arguments, i.e.,

q(Θ)
dΘ.
(18)
KL(q(Θ) p(Θ|Y)) =q(Θ) logp(Θ|Y)
Fig. 3 Illustration of the generative process of the complex array
response ai,k based on the DOA mixture model.

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 divergence is minimized when its twoarguments are exactly equal. Of course, solving this optimization problem is just as intractable as directly computing the posterior. Although it may appear that no progresshas been made, restricting q(Θ) to belong to a family of distributions of simpler form than p(Θ|Y) allows us to obtain
principled approximate solutions.
For our model, we define the set of approximate distributions Q as those that factor as follows:
Q = q : q(A)q(S )q(Z)q(V)q(C)q(ρ) .

Fig. 4

4.

Plate notation of present overall generative model.

Approximate Posterior Inference

4.1 Variational Bayesian Approach
In this section, we describe an approximate posterior inference algorithm for our generative model based on variational inference [10]. The random variables of interest in
our model are
A = {ai,k }i,k
: complex array response for source i,: component of active source at (ωk , tl ),S = {sk,l }k,l
: index of active source at (ωk , tl ),Z = {zk,l }k,l: stick breaking proportion,
V = {vi }i: index of direction for source i,
C = {ci }iρ = (ρ1 , . . . , ρD ) : mixture weight for each DOA.
We denote the entire set of the above parameters as Θ. Our
goal is to compute the posterior p(Θ|Y) where Y = {yk,l }k,lis a set consisting of the time-frequency components of observed multichannel signals. By using the conditional distributions defined in 3.2 and 3.3, we can write the joint distribution p(Y, Θ) asp(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:(19)

This approximation is often called a naive mean-field approximation. To define q(A), q(V) and q(C), we needto construct distributions on the infinite sets {v1 , v2 , . . .},{a1,k , a2,k , . . .} and {c1 , c2 , . . .}. For this approach to be
tractable, we truncate the variational distribution at somevalue N ∗ by setting q(vN ∗ = 1) = 1. The mixture proportions πi for i > N ∗ will thus be zero, and we can ignore ai,kand ci for i > N ∗ . In practice, we set N ∗ at a suﬃ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 N ∗ becomeslarger, the approximations to the true posterior become moreaccurate.
4.2 Coordinate Ascent
We now present an algorithm for solving the optimization
problem described in (17) and (19). Unfortunately, the optimization problem is non-convex, and it is intractable to findthe global optimum. However, we can use a simple coordinate ascent algorithm to find a local optimum. Notice that(18) can be written as
KL(q(Θ) p(Θ|Y))q(Θ)
dΘ + log p(Y).
=q(Θ) logp(Y, Θ)(20)

As the log evidence log p(Y) is fixed with respect to q(Θ),
minimizing the first term, which is known as the (negative) variational free energy, minimizes the KL divergenceof p(Θ|Y) from q(Θ). The algorithm optimizes one factor
in the mean-field approximation of the posterior at a time
while fixing all other factors. It can be shown using the calculus of variations that the “optimal” distribution for eachof the factors can be expressed as:q̂(X) ∝ exp EΘ\X [log p(Y, Θ)],(21)

IEICE TRANS. FUNDAMENTALS, VOL.E96–A, NO.10 OCTOBER 2013

1932

where X indicates one of the factors and EΘ\X [log p(Y, Θ)] is
the expectation of the joint probability of the data and latentvariables, taken over all variables except X. The mean-field
update equations for the variational distributions are givenin the following form:NC (ai,k ; mi,k , Γi,k ),(22)

NC (sk,l ; μk,l , σ2k,l ),(23)

Categorical(zk,l ; φk,l ),(24)

Beta(vi ; γi,0 , γi,1 ),(25)

Categorical(ci ; ψi ),(26)

q̂(ρ) ← Dirichlet(ρ; ζ1 , . . . , ζD ),(27)

q̂(A) ←i,kq̂(S ) ←k,lq̂(Z) ←

k,lq̂(V) ←iq̂(C) ←

i

where
Γ−1
i,k =




[φk,l ]i |μk,l |2 + σ2k,l Σk(n) −1 + Σk(a) −1 ,lmi,k = Γi,k Σk(n) −1



(28)

[φk,l ]i μ∗k,l yk,ll+ Σk(a) −1[ψi ]d h(ϑd , ωk ) ,

(29)

d




1
(n) −1
H
=
[φ
]
tr
m
m
+
Γ
,
Σ
i
i,k
i,k
k,l
i,k
k
σ2k,l
i


[φk,l ]i mHi,k Σk(n) −1 yk,l ,μk,l = σ2k,l

(30)
(31)

i


ϕi,k,l = exp 2Re μk,l yHk,l Σk(n) −1 mi,k− (|μk,l |2 + σ2k,l )tr mi,k mHi,k + Γi,k Σk(n) −1+ Ψ(γi,0 ) − Ψ(γi,0 + γi,1 )+Ψ(γ j,1 ) − Ψ(γ j,0 + γ j,1 ) ,

i−1

j=1

ϕi,k,l,[φk,l ]i = N ∗i =1 ϕi ,k,l
γi,0 = 1 +[φk,l ]i ,

(32)
(33)

k,l
N 


N∗
β0 
+
ψi,d .
ζd =
D i=1

(36)

[·]i denotes the i-th element of a vector, tr[·] denotes the traceof a matrix, and Ψ(·) denotes the digamma function.
Finally, the STFT components of the i-th separated signal can be obtained by φi,k,l μk,l . Since q(Θ) is an approximation to the true posterior p(Θ|Y), φi,k,l μk,l corresponds to an
approximation of the minimum mean squareerror estimator
of the i-th source signal, i.e., ŝi,k,l = E 1[zk,l = i]sk,l |YE[zk,l = i|Y]E[sk,l |Y] = φi,k,l μk,l , where 1[·] denotes the indicator function that takes the value 1 if its argument is trueand 0 otherwise.
5.

Related Work

It has been brought to our attention that several papers related to this work have been presented independently by different authors [11], [12], after the publication of our conference paper [17].
Taghia et al. proposed in [11] a variational Bayesian
approach for the approximate inference of the parameters of
Sawada’s model [7]. This method requires a model selection procedure to estimate the number of sources. The variational lower bound (an approximation of the log marginal
likelihood) is used as the criterion for the model selection.
By contrast, our method avoids selecting the number ofsources explicitly and instead pushes this task down into
the inference algorithm, incorporating automated Bayesianmodel selection in the inference procedure.
Otsuka et al. proposed in [12] a combined source separation and source localization method based on a Bayesian
approach. Similar to the idea described in 3.3, they described the prior distribution of the spatial covariance matrixof each source using a mixture model of Wishart-distributed
spatial covariance matrices, each of which is associated witha latent DOA variable. By contrast, the present model
describes the prior distribution of the frequency array response of each source using a mixture model of Gaussiandistributed frequency array responses, each of which is associated with a latent DOA variable. Although these models
appear slightly diﬀerent, the concepts are very similar. It
would be interesting to compare the diﬀerent ways in which
the use of these models aﬀect the source separation ability.
This should be investigated in the future.
6.

Experiment∗γi,1 = α0 +[φk,l ]i ,(34)

j=i+1 k,l
2Re h(ϑd , ωk )H Σk(a) −1 mi,ki,d = expk− tr h(ϑd , ωk )h(ϑd , ωk )H Σk(a) −1 + Ψ(ζd ) ,i,d[ψi ]d = D,
d =1 i,d(35)

We evaluated the performance of the proposed method interms of source separation ability.
We used stereo speech signals with a sampling rate of
16 kHz as test signals, which we obtained by mixing threespeech signals [13] (one male and two female speakers) using a measured room impulse response [14] (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. The geometry setting for the test mixtures

KAMEOKA et al.: BAYESIAN NONPARAMETRIC APPROACH TO BLIND SEPARATION OF INFINITELY MANY SPARSE SOURCES

1933

Fig. 5

Geometry setting for text mixtures.

is shown in Fig. 5. 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
16 ms.
Our preliminary experiments revealed that ζd was
likely to be trapped in local optima due to the spatial aliasingthat occurs at high frequencies. To avoid this, we adopted
the following procedure: We first ran the variational inference algorithm using only the low-frequency region of theobserved signals, after which we gradually increased the frequency range up to the Nyquist frequency during the iteration. Σk(n) and Σk(a) were set respectively at I and 10−1.5 × I.
D was set at 180. All the variational parameters were initialized randomly. After convergence, each separated signal
was obtained by multiplying μk,l by φi,k,l . Figure 6 showssome examples of the spectrograms of separated source signals obtained with the present method. Owing to the property of the stick-breaking process, sources that are dominant at many time-frequency points are likely to be assigned
small indices. Thus, we considered the separated signals of
indices i = 1, 2, 3 to be the source estimates, which we used
for the evaluation. We chose Sawada’s method described in[7] as a comparison. The source code was kindly provided
to us by its author. In this method, the number of sourcesmust be specified manually. The following results report
the performance in terms of the Signal-to-Distortion Ratio(SDR) [15]. The SDR is expressed in decibels (dB), and ahigher SDR indicates superior quality.
The present method was tested with various settingsof the truncation level N ∗ . As for Sawada’s method, it
was tested with various settings of the assumed number
N of sources. Figure 7 shows the average SDRs obtainedwith Sawada’s and the present methods with various N andN ∗ settings. As expected, the performance of the present
method improves with increasing N ∗ , while that of Sawada’s
method deteriorates significantly when the assumed number of sources departs from the actual number. It is worthemphasizing that Sawada’s method requires prior knowledge of the exact number of sources to achieve good separation, while the present approach does not. With Sawada’s
method, each source signal tends to be further separated into
multiple components when the specified number of sources
is greater than the actual number. Unlike with Sawada’s

Fig. 6

Separated source signals obtained with present method.

Fig. 7 SDRs obtained with conventional and present methods with
diﬀerent N and N ∗ settings.

method, such a phenomenon did not occur with the present
method. This is because the present method searches for an
explanation for an observed signal using a model with as
low a complexity as possible.
Figure 8 shows an example of the estimate of the interchannel phase diﬀerence of each source. The slope ofarg([mi,k ]2 /[mi,k ]1 ) along the frequency direction k corresponds to the DOA estimate of source i. Hence, if the DOAsof the sources were estimated correctly and the permutation
problem was successfully solved, arg([mi,k ]2 /[mi,k ]1 ) withthe same index i should lie on a straight line along k. Figure 8 thus reveals that the DOA estimation and permutation alignment both worked successfully with the present
method.
It is important to note that the DOA mixture model
described in 3.3 assumes a completely anechoic environment. Moreover, the assumption of the sparseness of speech

IEICE TRANS. FUNDAMENTALS, VOL.E96–A, NO.10 OCTOBER 2013

1934

Fig. 8 An example of the estimate of the inter-channel phase diﬀerence of each source plotted in adiﬀerent color. The slope of arg([mi,k ]2 /[mi,k ]1 ) along the frequency direction k corresponds to the DOA
estimate of source i. For each i, arg([mi,k ]2 /[mi,k ]1 ) lies on a straight line along k. This means that boththe DOA estimation and permulation alignment worked successfully with the present method.

described in Sect. 2 does not hold in reverberant environments. Indeed, the performance of the present method waspoor when using real data recorded in reverberant environments. This implies the need for a combined generativemodel of the present model and a reverberation process, e.g.,
by adopting the idea described in [16].
7.

Conclusion

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
an estimation of the number of sources, source separation
based on the sparseness of speech, an estimation of the direction of arrival of each source, and permutation alignment.
We confirmed experimentally that reasonably good separations were obtained with the present method without specifying the number of sources.
Blind separation of moving sound sources is an important task to be solved for real applications. One possible
extension to the present method involves representing the
DOA indicator variable ci as a time sequence, i.e., {ci,l }1≤l≤L ,and modeling its generative process using a hidden Markov
model. This is currently under development.
It should be noted that this paper is an extended journalversion of our conference papers [17], [18].
Acknowledgments
We thank Dr. Hiroshi Sawada (NTT) for providing us withthe source code of his method. We also thank Mr. Yosuke
Izumi (previously with The University of Tokyo) for providing us with Fig. 1.
References[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,” Adv.
Neural Info. Process. Syst., pp.953–960, 2006.
[4] S. Araki, H. Sawada, R. Mukai, and S. Makino, “Underdeterminedblind 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 two-stage blindsource separation combining SIMO-ICA and binary masking,” Proc.
9th International Workshop on Acoustic Echo and Noise Control(IWAENC 2005), pp.229–232, 2005.
[6] Y. Izumi, N. Ono, and S. Sagayama, “Sparseness-based 2ch BSS
using the EM algorithm in reverberant environment,” Proc. 2007
IEEE Workshop on Applications of Signal Processing to Audio andAcoustics (WASPAA 2007), pp.147–150, 2007.
[7] H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive blind source separation via frequency bin-wise clustering andpermutation alignment,” IEEE Trans. Audio Speech Language Process., vol.19, no.3, pp.516–527, 2010.
[8] J. Sethuraman, “A constructive definition of Dirichlet priors,” Statistica Sinica, vol.4, pp.639–650, 1994.
[9] T.S. Ferguson, “A Bayesian analysis of some nonparametric problems,” Annals of Statistics, vol.1, no.2, pp.209–230, 1973.
[10] D.M. Blei and M.I. Jordan, “Variational inference for Dirichlet process mixtures,” Journal of Bayesian Analysis, vol.1, no.1, pp.121–
144, 2006.
[11] J. Taghia, N. Mohammadiha, and A. Leijon, “A variational Bayes
approach to the underdetermined blind source separation with automatic determination of the number of sources,” Proc. 2012 IEEE
International Conferences on Acoustics, Speech and Signal Processing (ICASSP2012), pp.253–256, March 2012.
[12] T. Otsuka, K. Ishiguro, H. Sawada, and H.G. Okuno, “Bayesian unification of sound source localization and separation with permutation resolution,” Proc. Twenty-Sixth AAAI Conference on Artificial
Intelligence (AAAI-12), pp.2038–2045, July 2012.
[13] 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,” Transactions on Speech Communication,
pp.357–363, 1990.
[14] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada,“Acoustical sound database in real environments for sound scene
understanding and hands-free speech recognition,” Proc. 2nd International Conference on Language Resources & Evaluation (LREC
2000), pp.965–968, 2000.
[15] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio Speech
Language Process., pp.1462–1469, 2006.

KAMEOKA et al.: BAYESIAN NONPARAMETRIC APPROACH TO BLIND SEPARATION OF INFINITELY MANY SPARSE SOURCES

1935[16] T. Yoshioka, T. Nakatani, M. Miyoshi, and H.G. Okuno, “Blind Separation and Dereverberation of Speech Mixtures by Joint Optimization,” IEEE Trans. Audio Speech Language Process., vol.19, no.1,
pp.69–84, 2011.
[17] H. Kameoka, M. Sato, T. Ono, N. Ono, and S. Sagayama, “Bayesiannonparametric approach to underdetermined sparse BSS,” Proc. ASJ
meeting, 1-1-19, pp.713–716, March 2012 (in Japanese).
[18] H. Kameoka, M. Sato, T. Ono, N. Ono, and S. Sagayama, “Blindseparation of infinitely many sparse sources,” Proc. 13th International Workshop on Acoustic Signal Enhancement (IWAENC 2012),
Sept. 2012.

−



−1
(ai,k − h(ϑci , ωk ))H Σ(a)(ai,k − h(ϑci , ωk ))
k

i,k

A

=


i,k

1[zk,l = i] 2Re[s∗k,l aHi,k Σk(n) −1 yk,l ]l− |sk,l |2 aHi,k Σk(n) −1 ai,k+H (a) −1−11[ci = d] 2Re[aHi,k Σ(a)
h(ϑ
,
ω
)]
−
a
Σ
a
d
k
i,k ,
i,k k
k

d

(A· 8)x

Appendix A: Probability Density Functions
The circular complex normal distribution with mean μ andcovariance P is defined as1
NC (x; μ, P) =exp (x − μ)H P−1 (x − μ)(A· 1)|πP|= exp − xH P−1 x + 2Re μH P−1 x− μH P−1 μ − log |πP| , (A· 2)

where = denotes equality up to a term independent of x and1[·] denotes the indicator function that takes the value 1 if its
argument is true and 0 otherwise. The expectation of (A· 8)
taken over all the variables except A is given asE[zk,l = i] 2Re[E[s∗k,l ]aHi,k Σk(n) −1 yk,l ]i,kl− E[|sk,l |2 ]aHi,k Σk(n) −1 ai,k−1
+ E[ci = d] 2Re[aHi,k Σk(a) −1 h(ϑd , ωk )] − aHi,k Σ(a)
a
,
i,k
k
d

(A· 9)and soE[x] = μ,(A· 3)

E[xxH ] = P + μμH ,(A· 4)

where E[·] denotes expectation.
The Dirichlet distribution of order I with parameters
α1 , . . . , αI > 0 is defined asII
Γ(αi ) xαi −1 ,(A· 5)
Dirichlet(x; α) = i=1
I
Γ( i=1αi ) i=1 i
Ifor all x1 , . . . , xI > 0 satisfying i=1
xi = 1. The density iszero outside this open (I − 1)-dimensional simplex. Γ(·) isthe gamma function. The mean and entropy of the Dirichletdistributed variable are given by
E[xi ] = Iαii =1 αiE[zk,l = i] = q(zk,l = i) = [φk,l ]i
E[ci = d] = q(ci = d) = [ψi ]d
E[sk,l ] = μk,l(A· 10)(A· 11)(A· 12)

E[|sk,l |2 ] = |μk,l |2 + σ2k,l .

(A· 13)

We can thus confirm that (A· 9) is equal up to a constantterm to a quadratic form − i,k (ai,k − mi,k )H Γ−1i,k (ai,k − mi,k )where mi,k and Γi,k are given by (29) and (28), respectively.
Hence, q(A) is shown to be the product of multivariate complex normal distributions with mean mi,k and covariance Γi,k
q(A) =NC (ai,k ; mi,k , Γi,k ).
(A· 14)i,k,E[log xi ] = Ψ(αi ) − Ψ

where(A· 6)

αi .

B.2 Update Equation for q(S )(A· 7)

i

Note that when I = 2 the Dirichlet distribution reduces to a
Beta distribution.

The term in log p(Y, Θ) depending on S is given as−(yk,l − azk,l ,k sk,l )H Σk(n) −1 (yk,l − azk,l ,k sk,l )k,lS

Appendix B: Derivation of Update Equations
Here we show how we obtained the variational update equations. Recall that each variational update equation can beobtained via the formula given by (21).
B.1 Update Equation for q(A)The term in log p(Y, Θ) depending on A is given as−(yk,l − azk,l ,k sk,l )H Σk(n) −1 (yk,l − azk,l ,k sk,l )
k,l

=−


i

k,l



1[zk,l = i] |sk,l |2 tr ai,k aHi,k Σk(n) −1− 2Re s∗k,l aHi,k Σk(n) −1 yk,l .

(A· 15)The expectation of (A· 15) taken over all the variables except
S is given as−
E[zk,l = i] |sk,l |2 tr E[ai,k aHi,k ]Σk(n) −1
k,l

where

i



− 2Re s∗k,l E[aHi,k ]Σk(n) −1 yk,l ,(A· 16)

IEICE TRANS. FUNDAMENTALS, VOL.E96–A, NO.10 OCTOBER 2013

1936

E[zk,l = i] = [φk,l ]i ,(A· 17)E[ai,k ] = mi,k ,(A· 18)

B.4 Update Equation for q(V)E[ai,k aHi,k ] = mi,k mHi,k + Γi,k .

(A· 19)The term in log p(Y, Θ) depending on V is given as

H

H



We can thus confirm that (A· 16)
 is equal up to a constantterm to a quadratic function − k,l |sk,l − μk,l |2 /σ2k,l where
μk,l and σ2k,l are given by (31) and (30), respectively. Hence,q(S ) is shown to be the product of complex normal distributions with mean μk,l and variance σ2k,lq(S ) =NC (sk,l ; μk,l , σ2k,l ).
(A· 20)k,l=log vzk,lThe term in log p(Y, Θ) depending on Z is given as− (yk,l − azk,l ,k sk,l )H Σk(n) −1 (yk,l − azk,l ,k sk,l )

i

+

=

i

j=i+1 k,l
(1 − v j )i−1
(A· 21)log(1 − v j ) .

j=1

The expectation of (A· 21) taken over all the variables except
Z is given as
1[zk,l = i] 2Re[E[sk,l ]yHk,l Σk(n) −1 E[ai,k ]k,li− E[|sk,l |2 ]tr[E[ai,k aHi,k ]Σk(n) −1 ]+ E[log vi ] +

i−1

j=1

q(Z) =

E[zk,l = j] + α0 − 1 log(1 − vi ) .

We can thus confirm that the exponential of (A· 27) is equal
up to a constant factor toq(V) =
Beta(vi ; γi,0 , γi,1 ),(A· 28)iwhere γi,0 and γi,1 are given by (33) and (34), respectively.
B.5 Update Equation for q(C)The term in log p(Y, Θ) depending on C is given as− (ai,k − h(ϑci , ωk ))H Σk(a) −1 (ai,k − h(ϑci , ωk ))+ log ρci
C

=


i

1[ci = d]d− h(ϑd , ωk ) Σk
(A· 23)(A· 24)i2Re[aHi,k Σk(a) −1 h(ϑd , ωk )]k
H (a) −1

where

Hence,

N∗ 


i,k


E[log(1 − v j )] , (A· 22)

E[log vi ] = Ψ(γi,0 ) − Ψ(γi,0 + γi,1 ),E[log(1 − v j )] = Ψ(γ j,1 ) − Ψ(γ j,0 + γ j,1 ).


1[zk,l = j] + α0 − 1 log(1 − vi ).

j=i+1 k,l

j=1

− |sk,l |2 tr[ai,k aHi,k Σk(n) −1 ] + log vi +(A· 26)

k,l

+

i

k,l

(α0 − 1) log(1 − vi )
1[zk,l = i] log vi



i

z
k,l −1
1[zk,l = i] 2Re[sk,l yHk,l Σk(n) −1 ai,k ]

i=1

k,l
∞ 


k,l



j=1

The expectation of (A· 26) taken over all the variables except
V is given as 
E[zk,l = i] log vi(A· 27)

B.3 Update Equation for q(Z)Z∞(1 − v j ) +

 

k,l

+ log vzk,lzk,l −1h(ϑd , ωk ) + log ρd .

(A· 29)The expectation of (A· 29) taken over all the variables except
C is given as
1[ci = d]2Re[mHi,k Σk(a) −1 h(ϑd , ωk )]idk
H (a) −1

− h(ϑd , ωk ) Σkq(zk,l )h(ϑd , ωk ) + E[log ρd ] , (A· 30)k,lq(zk,l ) ∝ exp 2Re[μk,l yHk,l Σk(n) −1 mzk,l ,k− (|μk,l |

2

where

+ σ2k,l )tr[(mzk,l ,k mHzk,l ,k + Γzk,l ,k )Σk(n) −1 ]+ Ψ(γzk,l ,0 ) − Ψ(γzk,l ,0 + γzk,l ,1 )+zk,l −1j=1Ψ(γ j,1 ) − Ψ(γ j,0 + γ j,1 ) .

(A· 25)E[log ρd ] = Ψ(ζd ) − Ψ




ζd .

(A· 31)d

We can thus confirm that the exponential of (A· 30) is equal
up to a constant factor toq(C) =
Categorical(ci ; ψi ),(A· 32)i

KAMEOKA et al.: BAYESIAN NONPARAMETRIC APPROACH TO BLIND SEPARATION OF INFINITELY MANY SPARSE SOURCES

1937

where [ψi ]d is given by (35).
B.6 Update Equation for q(ρ)
The term in log p(Y, Θ) depending on ρ is given as
log ρci + (β0 /D − 1)log ρd(A· 33)
i

=


i

d

1[ci = d] log ρd + (β0 /D − 1)dlog ρdd

The expectation of (A· 33) taken over all the variables except
ρ is given as [ψi ]d + β0 /D − 1 log ρd .
(A· 34)di

We can thus confirm that the exponential of (A· 34) is equal
up to a constant factor toq(ρ) = Dirichlet(ρ; ζ1 , . . . , ζD ),(A· 35)where ζd is given by (36).

Hirokazu Kameoka
received B.E., M.S.
and Ph.D. degrees all from the University of Tokyo, Japan, in 2002, 2004 and 2007, respectively. He is currently a research scientist at
NTT Communication Science Laboratories anda Visiting Associate Professor at the Universityof Tokyo. His research interests include computational auditory scene analysis, statistical signal
processing, speech andmusic processing, andmachine learning. He is a member of IEEE, theInformation Processing Society of Japan (IPSJ)
and the Acoustical Society of Japan (ASJ). He received 13 awards over the
past 10 years, including the Yamashita Memorial Research Award in 2005
from IPSJ, the Itakura Prize Innovative Young Researcher Award in 2007
and the Awaya Prize Young Researcher Award in 2008 from ASJ, IEEE
Signal Processing Society 2008 SPS Young Author Best Paper Award in
2009, and IEICE ISS Young Researcher’s Award in Speech Field in 2011.
Misa Sato
received her B.E. degree fromthe University of Tokyo, Japan, in 2012. She iscurrently an MSc student at the Graduate Schoolof Engineering, the University of Tokyo, Japan.
Her research interests include natural language
processing and machine learning.

Takuma Ono
received the B.E. and M.S.
degrees from the University of Tokyo, Japan, in2010 and 2012, respectively. Since then, he has
been with Zenkyoren (National Mutual Insurance Federation of Agricultural Cooperatives).

Nobutaka Ono
received the B.E., M.S., andPh.D. degrees in Mathematical Engineering andInformation Physics from the University of Tokyo, Japan, in 1996, 1998, 2001, respectively.
He joined the Graduate School of Information
Science and Technology, the University of Tokyo, Japan, in 2001 as a Research Associate andbecame a Lecturer in 2005. He moved to the
National Institute of Informatics, Japan, as an
Associate Professor in 2011. His research interests include acoustic signal processing, specifically, microphone array processing, source localization and separation,
music signal processing, audio coding and watermarking, and optimization algorithms for them. He is the author or co-author of more than 90
articles in international journals and conference proceedings. He was a
Tutorial speaker at ISMIR2010 and will organize a special session in EUSIPCO2013. Dr. Ono has been an Associate Editor of the IEEE Transactions on Audio, Speech and Language Processing since 2012 and an Associate Editor of Acoustic Science and Technology since 2012. He is a chairof SiSEC (Signal Separation Evaluation Campaign) committee in 2013. He
is a senior member of the IEEE Signal Processing Society and a memberof the Acoustical Society of Japan (ASJ), the Information Processing Society of Japan (IPSJ), the Institute of Electrical Engineers of Japan (IEEJ),
and the Society of Instrument and Control Engineers (SICE). He received
the Sato Paper Award and the Awaya Award from ASJ in 2000 and 2007,
respectively, and received the Igarashi Award at the Sensor Symposium on
Sensors, Micromachines, and Applied Systems from IEEJ in 2004.

Shigeki Sagayama
received the B.E., M.S.,
and Ph.D. degrees from the University of Tokyo, Tokyo, Japan, in 1972, 1974, and 1998, respectively, all in mathematical engineering andinformation physics. He joined Nippon Telegraph and Telephone Public Corporation (currently, NTT) in 1974 and started his career inspeech analysis, synthesis, and recognition at
NTT Labs in Musashino, Japan. From 1990,
he was Head of the Speech Processing Department, ATR Interpreting Telephony Laboratories,
Kyoto, Japan where he was in charge of an automatic speech translation
project. In 1993, he was responsible for speech recognition, synthesis, anddialog systems at NTT Human Interface Laboratories, Yokosuka, Japan.
In 1998, he became a Professor of the Graduate School of Information
Science, Japan Advanced Institute of Science and Technology (JAIST),
Ishikawa. In 2000, he was appointed Professor at the Graduate Schoolof Information Science and Technology (formerly, Graduate School of Engineering), the University of Tokyo. On his retirement from the Universityof Tokyo in 2013, he became a Project Professor at the National Instituteof Informatics (NII). His major research interests include the processing
and recognition of speech, music, acoustic signals, handwriting, and images. He was the leader of anthropomorphic spoken dialog agent project(Galatea Project) from 2000 to 2003. Prof. Sagayama received the National
Invention Award from the Institute of Invention of Japan in 1991, the ChiefOﬃcial’s Award for Research Achievement from the Science and Technology Agency of Japan in 1996, and other academic awards including Paper
Awards from the Institute of Electronics, Information and Communications
Engineers, Japan (IEICEJ) in 1996 and from the Information Processing
Society of Japan (IPSJ) in 1995. He is a member of the Acoustical Societyof Japan, IEICEJ, and IPSJ.