Stereo Source Separation and Source Counting
with MAP Estimation with Dirichlet Prior
Considering Spatial Aliasing Problem
Shoko Araki, Tomohiro Nakatani, Hiroshi Sawada, and Shoji Makino
NTT Communication Science Laboratories, NTT Corporation
2-4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japan{shoko,nak,sawada,maki}@cslab.kecl.ntt.co.jp

Abstract. In this paper, we propose a novel sparse source separation
method that can estimate the number of sources and time-frequency
masks simultaneously, even when the spatial aliasing problem exists. Recently, many sparse source separation approaches with time-frequency
masks have been proposed. However, most of these approaches require information on the number of sources in advance. In our proposed method,
we model the phase diﬀerence of arrival (PDOA) between microphones
with a Gaussian mixture model (GMM) with a Dirichlet prior. Thenwe estimate the model parameters by using the maximum a posteriori(MAP) estimation based on the EM algorithm. In order to avoid one
cluster being modeled by two or more Gaussians, we utilize a sparse
distribution modeled by the Dirichlet distributions as the prior of the
GMM mixture weight. Moreover, to handle wide microphone spacingcases where the spatial aliasing problem occurs, the indeterminacy ofmodulus 2πk in the phase is also included in our model. Experimental
results show good performance of our proposed method.
Keywords: Dirichlet distribution, prior, number of sources, blind source
separation, sparse, spatial aliasing problem.

1

Introduction

Blind source separation (BSS) is an approach for estimating source signals thatuses only the mixed signal information observed at each microphone. The BSS
technique for speech dealt with in this paper has many applications, including
the hands-free teleconference systems and preprocessing for an automatic speech
recognizer.
Let us formulate the task. Suppose that Ns ≥ 2 speech sources s1 , . . . , sNs areconvolutively mixed and observed at Nm microphones,
Ns 
(1)
xj (t) = i=1
l hji (l) si (t − l), j = 1, . . . , Nm ,where hji (l) represents the impulse response from source i to microphone j.
Our goal is to obtain estimates yi of each source signal si from the microphoneobservations xj without information about the number of sources Ns , the speechsources si or the mixing process hji .
T. Adali et al. (Eds.): ICA 2009, LNCS 5441, pp. 742–750, 2009.
c Springer-Verlag Berlin Heidelberg 2009

Stereo Source Separation and Source Counting

743

Two approaches have been widely studied and employed to solve the BSS
problem: one is based on independent component analysis (ICA) (e.g., [1]) andthe other relies on the sparseness of source signals (e.g., [2]). In this paper, we focus on the latter approach, more speciﬁcally, the time-frequency mask approach[2,3]. With the time-frequency mask approach, we classify the phase diﬀerenceof arrival (PDOA) between microphone observations, and separate each signal
by collecting the observation signal at time-frequency points in each cluster.
In previous work [2,3], to automatically ﬁnd clusters, the number of sources
Ns is assumed to be known. However, in real situations we usually cannot obtain
information on the number of sources Ns in advance. Especially for an underdetermined case (Ns > Nm ), the source counting is diﬃcult, and few papers
have dealt with this problem. Moreover, when the microphone spacing is large,
the spatial aliasing problem occurs. This problem makes it diﬃcult to classify
the PDOA because the phase has the indeterminacy of modulus 2πk in highfrequencies. [4] considered the spatial aliasing problem in a time-frequency mask
approach, however, the number of sources Ns should be known.
In this paper, we propose a novel sparse source separation method that canestimate the number of sources and time-frequency masks simultaneously, evenwhen spatial aliasing occurs. We model the PDOA distribution with a Gaussian mixture model (GMM) with a Dirichlet prior [5], and estimate the model
parameters by using the EM algorithm. In order to avoid one cluster being modeled by two or more Gaussians, thus making it possible to estimate the numberof sources correctly, we propose utilizing a sparse distribution modeled by the
Dirichlet distribution as the prior of the GMM mixture weight. The authorsof [6,7] also derived the EM algorithm, however, they still needed to know thenumber of sources Ns in advance. On the other hand, our proposed algorithm
does not require information on the source number, thanks to the weight prior.
Because the indeterminacy of 2πk in phase is modeled in our GMM, we can also
overcome the diﬃculty in the PDOA clustering even in a spatial aliasing case.
The experimental results with a wide microphone spacing (20 cm) show thatour proposed method can estimate the number of sources and can separate signals
by time-frequency masks obtained by the posterior probability for each cluster.

2

Mixing and Separation Processes

This paper employs a time-frequency domain approach. With an F -point shorttime Fourier transform (STFT), (1) is converted into:Nshji (f )si (n, f ),(2)
xj (n, f ) = i=1
where hji (f ) is the frequency response from source i to microphone j, si (n, f )is the STFT of a source si . f ∈ {0, F1 fs , . . . , FF−1 fs } is a frequency (fs is the
sampling frequency) and n(= 0, · · · , N − 1) is a time-frame index.
In this paper, we assume the sparseness of the sources [2]:(3)xj (n, f ) ≈ hji (f )si (n, f ),where si (n, f ) is a dominant source at the time-frequency slot (n, f ). This isapproximately true for speech signals in the time-frequency domain [2,3].

744

2.1

S. Araki et al.

Separation Method

In this paper, we assume that hji (f ) in (2) is modeled by an anechoic model(e.g., eq. (13) of [3]), that is, the PDOA between microphones is given as:
x1 (n, f )d cos ϕ(n, f )
,
(4)
arg
= 2πf τ (n, f ) = 2πfx2 (n, f )vwhere τ (n, f ) = d cos ϕ(n, f )/v is the time diﬀerence of arrival (TDOA), ϕ(n, f )is the dominant source direction at the time-frequency (n, f ), and d and v denotethe microphone spacing and the sound speed.
First, by assuming the source sparseness, we calculate the PDOA at eachtime-frequency slot by the left-side of (4). Then, by considering the frequency
dependence in the PDOA, we classify the PDOA values in some way. For example, if there is no spatial aliasing problem and we know the number ofsources, the k-means clustering algorithm can be applied to TDOA τ (n, f ) =1
2πf arg [x1 (n, f )/x2 (n, f )]. Our method, which considers the aliasing problem
and the unknown source number, is introduced in the following section.
Finally, we estimate the separated signals yi (n, f ) with time-frequency masks
Mi (n, f ), which extract time-frequency points of members in the i-th cluster:(5)yi (n, f ) = x1 (n, f )Mi (n, f ).

3

Proposed Method

3.1

Problems in PDOA Clustering

The ﬁrst problem for the PDOA clustering is the spatial aliasing problem. As
can be seen in (4), when the frequency f or the microphone spacing d are large,
arg [x1 (n, f )/x2 (n, f )] = 2πf d cos ϕ(n, f )/v exceeds ±π. However, since the arg
operation has the indeﬁniteness of modulus 2πk, (4) should be:
x1 (n, f )2πf τ (n, f ) = arg+ 2πk = o(n, f ) + 2πk,(6)x2 (n, f )where o(n, f ) = arg [x1 (n, f )/x2 (n, f )], −π ≤ o(n, f ) < π, and k is an integer.
Note that we can observe just o(n, f ), and k is unknown when the source direction ϕ is unknown. This is the spatial aliasing problem. Figure 1 gives anexample of the observed PDOA o(n, f ) for a wide microphone spacing of 20 cm.
In the next subsection, k in (6) is considered as a hidden variable.
The second problem occurs when we apply a GMM ﬁtting method for an
unknown number of mixtures. Figure 2 shows an example. Here we have twoclusters in the histogram (Fig. 2(a)). Figure 2(b) shows the ﬁtting result of
GMM of eight Gaussians, to Fig. 2 (a). From Fig. 2(b), we can see that multiple
Gaussians are ﬁt to each cluster. However, we expect just one Gaussian for eachpeak, in order to estimate the number of sources by counting the number ofdominant Gaussians. In this paper, in order to avoid the case where one cluster
is modeled by two or more Gaussians, we propose utilizing a sparse distribution
for the prior of the GMM mixture weight parameter in the next subsection.

Stereo Source Separation and Source Counting

745

Fig. 1. Example PDOA for a microphone spacing of 20 cm and a sampling rate of 16
kHz. The PDOA for a source at 70◦ and a source at 150◦ are drawn individually forillustrative purposes.

Fig. 2. Example GMM ﬁtting result with and without prior. (a) Histogram of two
Gaussians, (b) estimated Gaussians of GMM without prior (φ = 1.0), (c) estimatedGaussians of GMM with prior (φ = 0.9).

3.2

Probabilistic Model

To begin with, let us consider that we observe one source from one direction.
Hereafter, a notation onf = o(n, f ) is utilized. Because the spatial aliasing issuein (6) can be considered as a phase wrapping problem, we can model the PDOA
with a Gaussian distribution by considering the unwrapped data onf + 2πk. In
other words, the phase wrapping process can be modeled by summing the Gaussians at intervals of 2π. That is, we assume that the PDOA follows a wrapped
Gaussian distribution [8],
Kf


Kf


1
√
p(onf ; μ, σ) =p(onf , k; μ, σ) =
exp
2πσ 2
k=−Kf
k=−Kf




−(onf + 2πk − 2πf μ)2
,
2σ 2
(7)

where −π ≤ onf < π, μ gives us the expectation value of the TDOA τ ofthe source, σ 2 is the variance of the PDOA, and k is an integer to handle thespatial aliasing (6). The value Kf is a frequency dependent integer, and it canbe determined if we know the microphone spacing d and the frequency f . If we
do not know d, we can set a suﬃciently large value for Kf for all frequencies.
This model is inspired by a wrapped Gaussian model [8].

746

S. Araki et al.

In our observed mixture, we assume that there are a suﬃcient number ofsource signals from diﬀerent directions, where some are dominant and others aremuch less dominant. Each source is modeled by (7). We also assume that the
PDOA for an observed mixture follows a Gaussian mixture model (GMM):


Kf
M


−(onf + 2πk − 2πf μm )2αm
p(onf ; μm , σm ) =
exp
. (8)
2
2
2σm
2πσm
m=1
k=−Kf

We prepare a suﬃcient number M of Gaussians for our GMM model and estimate
2
and weight αm for each Gaussian m.
the mean μm , variance σm
In order to solve the second problem mentioned in Section 3.1, that is, inorder to model the observed PDOA data by allocating one Gaussian to eachsource, we assume the sparseness of the source directions, where each direction
is dominated by at most one source. For this purpose, as the prior of the mixture
weight, we employ the Dirichlet distribution:M1
αφ−1 ,(9)p(α) =
B(φ) m mM
where α = {α1 , · · · , αm , · · · , αM }, m αm = 1, 0 ≤ αm ≤ 1, and B(φ) is thebeta distribution (regularization term). When we set small hyper parameter φ(φ < 1), the prior takes a larger value as the number of mixture weights whose
values are close to zero increases, which is desirable for representing the sparseness of the source direction [5]. In addition, the Dirichlet distribution is known
to be a conjugate prior of the mixture weight [5], and it can be incorporated into
the GMM ﬁtting approach in a computationally eﬃcient manner.
Figure 2 (c) shows a GMM ﬁtting result with prior ((9) with φ = 0.9) for the
distribution in Fig. 2 (a). In spite of utilizing eight Gaussians, we can see thatjust two Gaussians are dominant in Fig. 2 (c). That is, using the prior, morecorrect GMM ﬁtting can be performed.
3.3

Cost Function Based on GMM

Let θ = {αm , μm , σm } be a model parameter set. The observations are o ={o11 , o12 , · · · , onf , · · · , oN F } and power values a = {a11 , a12 , · · · , anf , · · · , aN F },where anf = a(n, f ) = |x1 (n, f )|2 . In the following, Gaussian indices m and k in
the PDOA model (8) are assumed not to be observed, and therefore dealt withas hidden variables.
The cost function of the maximum a posteriori (MAP) estimation is deﬁned
based on a log of a joint probability density function (pdf) as
L(θ) = log p(o, θ) = log p(o|θ) + log p(α) + const.
=

N 
F

n

=

f

N 
F

n

f

f (anf ) log p(onf |θ) + log p(α) + const.

(10)
(11)

⎛
⎞
Kf
M


f (anf ) log⎝p(m, k, onf |θ)⎠ +log p(α)+const., (12)m k=−Kf

Stereo Source Separation and Source Counting

747

where



−(onf + 2πk − 2πf μm )2
αm
exp
.
(13)
p(m, k, onf |θ) = 
2
2
2σm
2πσm
We disregarded the priors of the model parameters except for α in (10), and

anf
(14)
f (anf ) = canf /nfgives a power weight (anf = |x1 (n, f )| in this paper) and controls the importance of the observation relative to the prior term (2nd term of (12)), where cis a control parameter.
M
In (12), the mixture weight α follows the Dirichlet distribution (9) and m αm
= 1, 0 ≤ αm ≤ 1 holds. For the sparse representation of the GMM, φ < 1 ispreferred for the Dirichlet distribution (9). Note that φ = 1 is equivalent to thecase without a prior for the mixture weight.
2

3.4

EM Algorithm

Here we derive an algorithm for estimating parameter θ by the EM algorithm.
The auxiliary function Q is given as(15)
Q(θ|θt ) = E log p(onf ; θ)|onf ; θt


t
=
p(m, k|onf , θ )f (anf ) log p(m, k, onf |θ) + log p(α), (16)
n

f

m

k

t

where θ is the estimate of the parameters after the t-th iteration, andp(m, k, onf |θt )p(m, k|onf , θt ) =  .
tm
k p(m, k, onf |θ )
t

(17)

t

)
)
By setting ∂Q(θ|θ= 0 and ∂Q(θ|θ= 0, we obtain
2
∂μm
∂σm
  
p(m, k|onf , θt )f (anf )(onf + 2πk)n
 f k 
(18)
μt+1
m =
t
n
f
k 2πf p(m, k|onf , θ )f (anf )
  
t
2
n
f
k p(m, k|onf , θ )f (anf )(onf + 2πk − 2πf μm )
2 t+1
  
(σm )
. (19)=tnf
k p(m, k|onf , θ )f (anf )MMoreover, by using the Lagrange multiplier method, m αm = 1 and (14), themixture weight is obtained as follows:
⎫
⎧
⎬
⎨



1
p(m, k|onf , θt )f (anf ) + (φ − 1) . (20)
αt+1
m =
⎭
c + M (φ − 1) ⎩ nfk

Since αm > 0, c > M (1 − φ) must hold from (20).
In the E-step we calculate (17), then in the M-step the parameters θ arecalculated by using (18), (19) and (20). Sometimes αm < 0 occurs. In such acase, we can factor out the corresponding Gaussian (by setting αm = , whereis a very small value) and recalculate the parameters.

748

3.5

S. Araki et al.

Source Counting

Thanks to the Dirichlet prior (9), most of the mixture weight parameter αmbecomes very close to zero and some have dominant values. Sometimes, someweight parameters αm do not come to zero suﬃciently because of very large
variance σm . So, we can determine the number of sources Ns by counting thenumber of Gaussians whose parameters meet conditions αm ≥ and σm ≤ th,where is a suﬃciently small threshod value and th is an appropriate thresholdvalue. = 0.2 and th = π/3 degrees are used in this paper.
3.6

Source Separation

The time-frequency mask Mm (n, f ) for the m-th separated source (see (5)) is
obtained by marginalizing the estimated pdf (17) with respect to k,
Mm (n, f ) = p(m|onf , θ) =Kf
p(m, k|onf , θ).

(21)

k=−Kf

The separated signal is obtained byym (n, f ) = x1 (n, f )Mm (n, f ) = x1 (n, f )p(m|onf , θ).

4Experiments

4.1

Experimental Setup(22)

We performed experiments with measured impulse responses hji in a room whose
reverberation time was 130 ms (see Fig. 9’s setup A of [4]). We utilized two
microphones whose spacing was 20 cm. The numbers of sources Ns were two andthree. Mixtures were made by convolving the measured room impulse responses
and 5-second English speech signals sampled at 16 kHz. The frame size F for
STFT was 1024 (64 ms), and the frame shift was 256 (16 ms).
In the EM algorithm, we utilized M = 8 Gaussians. From the microphonespacing and sampling rate, the aliasing problem occured above 850 Hz. In our
implementation, Kf = K = 5 was utilized for all frequencies f . For the comparison with an aliasing-unconsidered case, Kf = K = 0 for all frequencies f was
also tested. As the hyper parameter for (9), we utilized φ = 0.9 for our proposed
method and φ = 1.0 for a conventional EM algorithm that corresponds to thecase without any prior for the mixture weights. The number of iterations was
10, and the control parameter c for (14) was 5.
We evaluated the signal-to-interference ratio (SIR) as a separation performance measure, and the signal-to-distortion ratio (SDR) as a sound quality
measure. Their deﬁnitions can be found in [3]. We calculated SIR and SDR values for the separated sources that are counted as the sources by the method in
Section 3.5. We conducted 20 trials with diﬀerent speech source combinations
and location combinations, and then averaged the results.

Stereo Source Separation and Source Counting

749

Fig. 3. Example spectra of (A)(B) sources, (C)(D) observations, (E)(F) separated signals. Ns = 2, φ = 0.9 and K = 5.
Table 1. Experimental results. Input SIR was 0.0 [dB] (Ns = 2), and −3.1 [dB]
(Ns = 3).
Ns
2

3

4.2

Accuracy of N̂s estimation [%] Performance [dB]
φ K N̂s :1 2
3 4 5 6 7 8 Output SIR SDR
0.9 5
100
11.5
12.4
0.9 0
35 60 5
11.4
7.5
1.0 5
0 5 40 35 20
8.8
9.7
0.9 5
15 75 10
7.5
8.2
0.9 0
5
25 50 20
2.0
8.7
1.0 5
0 90 10
6.9
7.5

Results

Figure 3 shows the example spectra of sources, observations, and separated
sources for a two source case. The source directions were 70◦ and 150◦ , whose example PDOA is shown in Fig. 1. By comparing the source spectra Fig. 3 (A)(B)and the separation spectra Fig. 3 (E)(F), it can be seen that the spatial aliasingproblem does not occur in most frequencies. However, it is also seen that at
the frequencies where the PDOA of two sources lap over each other, say around
1500, 3000, 4500, 6000, 7500 Hz (see Fig. 1), the signals are not separated well.
Such phenomena can be seen in the separated spectra Fig. 3 (E)(F).
Table 1 reports the experimental results. In the table, φ = 0.9 means theresults with sparse prior (9) and φ = 1.0 indicates the results without a prior.
K = 5 and K = 0 mean the spatial aliasing is considered and unconsidered,respectively. The percentage values are shown where the method estimates thenumber of sources as N̂s within 20 trials. The average separation performance
results, SIR and SDR in dB, are also reported.
From Table 1, we can see that with the prior (φ = 0.9) by considering thealiasing (K = 5), the number of sources is almost perfectly estimated. On theother hand, without the prior, the number of sources is overestimated, and theaccuracy rate was quite low.

750

S. Araki et al.

As for the separation performance, we obtained better performance by using
the prior (φ = 0.9) than without the prior (φ = 1.0) when K = 5. When we
did not consider the spatial aliasing, K = 0, the separation performance was ofcourse poor, especially when Ns = 3.

5

Conclusion

We proposed a speech source separation method that can estimate both thenumber of sources and separation masks. We model the PDOA with a GMM,
where the phase indeﬁniteness in spatial aliasing cases is considered. We employ
the Dirichlet distribution as the prior of the GMM mixture weight to model
each cluster by a single Gaussian. Our experimental results show that the proposed method can estimate the number of sources correctly. We also conﬁrmedthat the proposed method gives good separation performance in a room withreverberation time of 130 ms.

References
1. Hyvärinen, A., Karhunen, J., Oja, E.: Independent Component Analysis. John Wiley& Sons, Chichester (2001)
2. Yilmaz, O., Rickard, S.: Blind separation of speech mixtures via time-frequency
masking. IEEE Trans. Signal Processing 52(7), 1830–1847 (2004)
3. Araki, S., Sawada, H., Mukai, R., Makino, S.: Underdetermined blind sparse source
separation for arbitrarily arranged multiple sensors. Signal Processing 77(8), 1833–
1847 (2007)
4. Sawada, H., Araki, S., Mukai, R., Makino, S.: Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source
separation. IEEE Trans. Audio, Speech and Language Processing 15(5), 1592–1604(2007)
5. Bishop, C.M.: Pattern recognition and machine learning. Springer, Heidelberg(2008)
6. Mandel, M., Ellis, D., Jebara, T.: An EM algorithm for localizing multiple soundsources in reverberant environments. In: Proc. Neural Info. Proc. Sys. (2006)
7. O’Grady, P., Pearlmutter, B.: Soft-LOST: EM on a mixture of oriented lines. In:
Puntonet, C.G., Prieto, A.G. (eds.) ICA 2004. LNCS, vol. 3195, pp. 430–436.
Springer, Heidelberg (2004)
8. Smaragdis, P., Boufounos, P.: Learning source trajectories using wrapped-phasehidden markov models. In: Proc. of WASPAA 2005, pp. 114–117 (2005)