Multiple Sound Source Localization Based on
Stochastic Modeling of Spatial Gradient Spectra
Natsuki Ueno1,2 and Hirokazu Kameoka1
1

NTT Communication Science Laboratories, Atsugi, Japan
2
Tokyo Metropolitan University, Hino, Japannatsuki.ueno@ieee.org, hirokazu.kameoka.uh@hco.ntt.co.jp

source localization, and by assuming time-frequency sparsityof the source signals, it can be applied to each time-frequency
component even in the presence of multiple sources [5].
However, its application to multiple source localization problem, including a clustering scheme of the estimated positions
obtained for each time-frequency component, is still to beestablished. We also refer to DNN-based methods [10], [11],whose approaches are different from those based only on
physical properties of the wave propagation, such as the SSCPDE-based method.
Focusing on the capability of the SSC-PDE-based method,
we extend its framework to multiple source localization problem in the presence of noises based on stochastic modelingof the observation. As well as the conventional SSC-PDEbased method, the proposed method is applicable to a small
microphone array as long as a sound pressure and its spatial
gradient at one observation point can be obtained, which has
an advantage in practical feasibility over other methods usedmainly with a large microphone array. We ﬁrst derive theexpectation-maximization (EM) algorithm [12] for a given
number of sources, and then extend it to the variational
Bayesian inference (VBI) algorithm [13] for an unknown number of sources. By numerical experiments of multiple source
localization under a reverberant environment, the performancecomparison of the proposed EM and VBI algorithms were
investigated with the baseline MUSIC algorithm.

Abstract—We propose source localization methods for multiple
sound sources. The proposed method requires only an observationof a sound pressure and its spatial gradient at one ﬁxed point,
which can be realized by a small microphone array. The key idea
is to utilize the partial differential equation relating the observed
signals and the source position, which was originally proposed forthe direct method for the single source localization problem. We
extend this framework using stochastic modeling and proposeda method for the mutliple source localization in the presenceof noises. Two source localization methods are proposed: oneis the expectation-minimization algorithm for a given numberof sources, and the other is the variational Bayesian inference
for an unknown number of sources. By numerical experiments,
the localization accuracies of the two proposed methods arecompared with the baseline method.
Index Terms—source localization, microphone array,
expectation-maximization algorithm, variational Bayesianinference algorithm.

I. I NTRODUCTION
Sound source localization technique using a microphonearray has a wide variety of applications such as sonar, robot
audition, and hearing aids [1], [2]. Especially, it is one of
the major challenges to localize multiple sources using asmall microphone array. Various source localization methods
have been proposed in the literature, such as Multiple Signal
Classiﬁcation (MUSIC) algorithm [3], Generalized CrossCorrelation methods with Phase Transform (GCC-PHAT) [4],method based on the sound source constraint partial differential equation (SSC-PDE) [5], stochastic methods considering
reverberant envionment [6]–[9], and methods using deep neural network (DNN) [10], [11].
In most source localization methods, a large array with alarge number of microphones is generally preferable since
the source positions are estimated on the basis of the timedifference of arrivals of signals observed by the multiple
microphones. On the other hand, the SSC-PDE-based method,which utilizes the partial differential equation relating theobserved signal and the source position, enables source localization only from an instantanesous observation of a soundpressure and its spatial gradient at one observation position.
Therefore, this method can be realized theoretically by using asmall microphone array with a small number of microphones.
The SSC-PDE is formulated for a direct method of single

II. PARTIAL -D IFFERENTIAL -E QUATION -BASEDS TOCHASTIC M ODELING FOR S OURCE L OCALIZATIONAs preliminaries for the proposed methods, this section describes the observation model for the multiple source localization problem, which relates the observation of a sound pressure
and its spatial gradient to the multiple source positions.
A. Sound Source Constraint Partial Differential Equation
We begin by breiﬂy introducing basic theories of the SSCPDE, which was proposed by Ando et al. [5]. Let r0 ∈ R3
be the position vector of the sound source. By assuming thefree-ﬁeld (spherical-wave) propagation from the point source,
the sound pressure at the observation position r ∈ R3 andtime t ∈ R, denoted by f (r, t), is represented by1∥r − r0 ∥f (r, t) =
g t−
,
(1)
∥r − r0 ∥cThis work was supported by JST CREST Grant Number JPMJCR19A3,
Japan.

ISBN: 978-1-6654-6798-8

31

EUSIPCO 2022

Authorized licensed use limited to: NTT. Downloaded on November 07,2023 at 09:28:18 UTC from IEEE Xplore. Restrictions apply.

with θ := {R, n}, i.e., θ denotes the source position, and
 

ωm
1
−1
R + j c  nx
ωm
1


−1
R + j c  ny  ,Cm (θ) := 
(7)
ωm
1

nz 
−1
R +j c1
ym := [F0,m , Fx,m , Fy,m , Fz,m ]T ,

2
2
2
2
Σm := diag σm
, σm
, σm
, σ0,m.

Fig. 1: Example of array conﬁguration. For example, thegradient of F with respect to x direction is approximated by
Fx,m ≈ (F1,m − F2,m )/2d.

(9)Since C(θ) is invertible from det(C(θ)) = −1, we obtain theprobabilistic distribution of the observed signal ym under the
given θ aswhere g is the source signal, c is the speed of sound, and ∥ · ∥denotes the Euclidean norm. Therefore, the spatial gradient off , denoted by ∇f , is given by




1
R
1
R
∇f (r, t) = 2 g t −n+ġ t −
n,
(2)
R
c
cR
c

p(ym |θ) = NC (ym ; 0, Cm (θ)−1 Σm Cm (θ)−H ),(10)where (·)−H denotes the inverse of the Hermitian conjugate
of a matrix.
C. Stochastic Modeling for Multiple Sound Sources
Various kinds of real-world acoustic signals such as humanspeech and music have sparsity in the time-frequency representation. Therefore, many practical situations satisfy well antime-frequency disjointness of the sources, which means eachtime-frequency component of the observed signal is assumed
to be dominated by at most one sound source.
Let m, l, and k ∈ N denote indices of the frequency,
time, and sound source, respectively, where k = 0 denotes
the noise signal and k ̸= 0 denotes the signal derived fromthe point source. We assume that only the zm,l th sourceis active at the frequency m and time l. Then, under the
given zm,l = k, the probability density function of theobservation of the sound pressure and its spatial gradientobtained by the short-time Fourier transform, denoted byym,l := [F0,m,l , Fx,m,l , Fy,m,l , Fz,m,l ]T ∈ C4 , is given in asimilar manner with (10) by

where R := ∥r − r0 ∥, n := (r − r0 )/∥r − r0 ∥, and ġ denotes
the temporal derivative of g. On the other hand, the temporal
derivative of f , denoted by f˙, is given by1∥r − r0 ∥f˙(r, t) =ġ t −.
(3)∥r − r0 ∥c
From (2) and (3), we obtain the partial differential equation1
1∇f (r, t) =f (r, t) + f˙(r, t) n,(4)Rcwhich includes only the observed signal and its spatial andtemporal gradients. This equation is referred to as the soundsource constraint partial differential equation.
B. Stochastic Modeling for Single Sound Source
The SSC-PDE can be rewritten in the frequency domain as1ω−∇F (r, ω) ++jF (r, ω)n = 0,
(5)
R
c

(k)

p(ym,l |zm,l = k, ϑ) = NC (ym,l ; 0, Λm,l )where F is the temporal Fourier transform of f , ω ∈ Rdenotes the angular frequency, and j denotes the imaginary
unit. Here, suppose the temporal spectra of the sound pressure
and its spatial gradient are observed by a small microphonearray, for example, as shown in Fig. 1. Various other types ofarray, such as commercially available ambisonic microphonearrays [14] and acoustic vector sensors [15], [16], can also
be used for this purpose. Let {ωm }m denote the set of
discrete frequencies, and F0,m and [Fx,m , Fy,m , Fz,m ]T denoterespectively the observation of F (r, ωm ) and ∇F (r, ωm ) at
a given ﬁxed observation position r. Then, we have where
ϵx,m , ϵy,m , ϵz,m ∈ C are the observation errors caused by thesensor noises and approximation of the spatial gradient by thespatial subtraction.
Here, we assume that ϵx,m , ϵy,m , ϵz,m follow independently the (circularly-symmetric) complex Gaussian distribu2
tion NC (0, σm) and that F0,m follows the complex Gaussian2
distribution NC (0, σ0,m). Then, from (5), we have
Cm (θ)ym ∼ N (0, Σm )

(8)

with
(k)
Λm,l :=(Cm (θ(k) )−1 Σm,l Cm (θ(k) )−H (k ̸= 0),2
Wm(k = 0)
νm,l


(k)
(k)2
2
2
2
Σm,l := diag σm
, σm
, σm
, σ0,m,l .

(11)

(k)

(12)
(13)

Here, θ(k) := {R(k) , n(k) } denotes the position of the kth
source, and ϑ denotes the set of model parameters consisting(k)22
2
of {θ(k) }k , {σm}m , {σ0,m,l }m,l,k , {νm,l}m,l . The covarianvce
2
matrix for the noise component is modeled by νm,l
Wm , where
Wm represents the normalized time-invariant covariance and2νm,lrepresents the power of the noise depending on both timeand frequency. For example, we can deﬁne Wm for the diffuse
noise ﬁeld in accordance with [17].
III. P ROPOSED M ETHODS
On the basis of stochastic modeling described in Sec. II-C,
we propose the EM and VBI algorithms for multiple soundsource localization.

(6)

32
Authorized licensed use limited to: NTT. Downloaded on November 07,2023 at 09:28:18 UTC from IEEE Xplore. Restrictions apply.

with V := {vk }∞k=0 . Therefore, the joint distribution p(Θ, Y )and the posterior distribution p(Θ|Y ) are obtained using theprior distribution p(ϑ) as

A. EM Algorithm for Given Number of Sources
Let K ∈ N denotes the number of sources, which isassumed to be given in this algorithm, and zm,l is assumed tofollow the categorical distribution:
(0)
(1)
(K)
zm,l ∼ Categorical(πm,l , πm,l , . . . , πm,l ).
(0)

(1)

k=0

with Y := {ym,l }m,l . Therefore, the source localization
problem can be reduced to the maximization problem of thelog likelihood function:
ϑ,π

(16)

q(Θ) = q(N )q(ρ)q(λ)q(ζ)q(γ)q(Z)q(V ).
N = {n(k) }k ,
ρ = {ρ(k)γ = {γm,l }m,l

(17)

(29)

(k)

(30)
(k)2

(31)

2
(γm,l = 1/νm,l
).

(32)

(ζm,l = 1/σ0,m,l ),k=1

q(ρ) =NR (ρ(k) ; µ̄(k) , η̄ (k) ),k=1
Yq(λ) =
Gamma(λm ; ᾱλ,m , β̄λ,m ),q(ζ) =

m

Y

YK ′
m,l

q(γ) =
Y
q(Z) =

(k)

(k)

(k)

Gamma(γm,l ; ᾱγ,m,l , β̄γ,m,l ),(37)

Categorical(zm,l ; {π̄m,l }∞k=0 ),(38)(k)

q(V ) =

(35)
(36)

m,l
m,l

(34)

Gamma(ζm,l ; ᾱζ,m,l , β̄ζ,m,l ),

k=1

Y

(18)

YK ′
k=0

Beta(vk ; ᾱv(k) , β̄v(k) ),(39)

where vMF(·), NR (·), and Gamma(·) denotes the von Mises–
Fisher distribution [20], real Gaussian distribution, and gamma
distribution, respctively. Note that the classes of the distributions from (33) to (37) are kept the same in the iterations by
deﬁning the initial distributions also in the above classes. Here,K ′ ∈ N represents the truncation of the stick-breaking process,which corresponds to the assumption q(zm,l = k) = 0 (k ≥K ′ + 1). Note that this truncation does not ﬁx the complexityof the model but simply restricts the function space for q(Θ)
to a certain extent. Therefore, we can aproximate q(Θ) wellby setting large K ′ . The update rules for parameters in (33)to (37) are summarized in Appendix B.

(20)

m,l

Moreover, from the stick-breaking process, we have


Y
Yk−1
p(Z|V ) =πzm,lπk = v k(1 − vk ) , (22)j=0

Beta(vk ; 1, βs )

),

YK ′

Then, the probability density function of the observed signals
Y under the given parameters ϑ and Z := {zm,l }m,l is givenbyY(zm,l )
p(Y |ϑ, Z) =NC (ym,l ; 0, Λm,l
).
(21)

k=1

= 1/R

(k)

The distributions q(N ), q(ρ), q(λ), q(ζ), q(γ), q(Z), and q(V )are updated alternately in the following forms:
YK ′q(N ) =vMF(n(k) ; ξ̄ (k) , κ̄(k) ),
(33)

j=0

p(V ) = Beta(v0 ; αb , βb )

(ρ

(k)

2
(λm = 1/σm),(k)Moreover, the probabilities {πk }∞k=1 for the point sources
are modeled as random variables following the stick-breaking
process [19] with a hyper parameter βs ∈ [0, ∞):Yk−1
πk = v k(1 − vk ),
(19)

Y∞

(28)

ζ = {ζm,l }m,l,k

Here,
P∞ πk ∈ [0, ∞) is a nonnegative parameter satisfying
k=0 πk = 1, which represents the probability that the kth
source is active. The probability π0 for the noise component ismodeled as a random variable following the Beta distribution
with hyper parameters αn , βn ∈ [0, ∞):

m,l

}k

λ = {λm }mThe number of sources are often unknown in practicalsituations. For such cases, we extend the observation modelin Sec. III-A to unknown number of sources using Dirichlet
process mixture model in a similar manner with [18].
We consider inﬁnitely many sources and assume zm,l follows the categorical distribution for countably inﬁnite indices:vk ∼ Beta(1, βs ).

(27)Here, the model parameters N , ρ, λ, ζ, and γ are deﬁned as

B. VBI Algorithm for Unknown Number of Sources

π0 = v0 ∼ Beta(αn , βn ).

(25)

is locally minimized under the constraint that q(Θ) can be
decomposed as

This maximization problem is difﬁcult to solve in a closed
form; however, the local solution can be obtained by the EM
algorithm [12], where ϑ and π are updated alternately. The
update rules are summarized in Appendix A.

zm,l ∼ Categorical({πk }∞k=0 ).

p(Θ|Y ) = p(Θ, Y )/p(Y ),where Θ := {ϑ, Z, V } denotes all the model parameters. The
posterior distribution p(Θ|Y ) is difﬁcult to obtain in a closed
form; however, its approximate distribution can be obtained by
the VBI algorithm [13]. In the VBI algorithm, the distribution
q(Θ) can be obtained iteratively so that the Kullback–Leibler(KL) divergence between q(Θ) and p(Θ|Y ), deﬁned asZq(Θ)
KL(q(Θ)|p(Θ|Y )) = q(Θ) logdΘ,(26)p(Θ, Y )(K)

maximize log p(Y |ϑ, π).

(24)

(14)

Here, πm,l , πm,l , . . . , πm,l ∈ [0, ∞) are nonnegative paramePK(k)ters satisfying k=0 πm,l = 1. Then, the likelihood function(k)with respect to ϑ and π := {πm,l }m,l,k is given by
Y XK
(k)
(k)
p(Y |ϑ, π) =πm,l NC (ym,l ; 0, Λm,l )
(15)
m,l

p(Θ, Y ) = p(Y |ϑ, Z)p(Z|V )p(V )p(ϑ),

(23)

33
Authorized licensed use limited to: NTT. Downloaded on November 07,2023 at 09:28:18 UTC from IEEE Xplore. Restrictions apply.

Proposed (VBI)
Proposed (EM1)
Proposed (EM2)MUSIC

Proposed (VBI)
Proposed (EM1)
Proposed (EM2)MUSIC(a) Reverberation time 0.5 s (b) Reverberation time 0.5 s
and signal length 2779 ms.
and signal length 1665 ms.

Fig. 2: Experimental settings. The blue circle denotes thecenter position of the microphone array, and the red crosses
denote the true source positions.
Proposed (VBI)
Proposed (EM1)
Proposed (EM2)MUSIC

Proposed (VBI)
Proposed (EM1)
Proposed (EM2)MUSIC

IV. N UMERICAL E XPERIMENTS
We conducted numerical simulations of multiple source
localization in a reverberant environment. A microphone arraywith seven microphones, whose conﬁguration is given by
Fig. 1 with d = 0.03 m, was located in the center of arectangular room, whose size was 6 m × 10 m × 4 m. The
reﬂection coefﬁcients were set as 0.7308 and 0.4566 so thatthe reverberation times from Sabine’s formula were 0.5 s
and 0.2 s, respectively. Three sound sources were locatedaround the microphone array, whose positions were (1, 0, 0) m,(−0.5, 0.87, 0) m, and (−0.5, −0.87, 0) m from the center ofthe room as shown in Fig. 2. Human speeches from SRVDB [21] were used as the source signals. The sampling
frequency of the microphones were 32 kHz, and the frame
length of the short-time Fourier transform were 64 ms withthe half overlap.
We compared the four algorithms: the proposed VBI algorithm (denoted by VBI), the proposed EM algorithm withthe correct number of sources K = 3 (denoted by EM1), theproposed EM algorithm with an incorrect number of sources
K = 6 (denoted by EM2), and the conventional MUSIC
algorithm as a baseline method. For a simple comparison
with the MUSIC algorithm, only the source directions were
evaluated and they were estimated on the xy-plane. The soundP(k)source was detected if m,l πm,l ym,l in the EM algorithm,
P(k)m,l π̄m,l ym,l in the VBI algorithm, or the angular spectrum
in the MUSIC algorithm was larger than the threshold value.
In the proposed VBI algorithm, the estimated source direction
was given by the posterior mean, i.e., ξ̄ (k) . The accuracy ofthe source localization was evaluated by the F-measure, where
the true positive was deﬁned as the number of the estimatedsources whose direction was within an angle ±τ from sometrue source direction, the false positive was deﬁned as thenumber of the estimated sources whose direction was notwithin an angle ±τ from any true source direction, and the true
negative was deﬁned as the number of the true sources whose
direction was not within an angle ±τ from any estimatedsource direction.
Figure 3 shows the F-measure against the angle τ . Here,for each τ , different threshold values for the source detection(c) Reverberation time 0.2 s (d) Reverberation time 0.2 s
and signal length 2779 ms.
and signal length 1665 ms.

Fig. 3: Estimation accuracy with respect to source directions.

were investigated, and the highest F-measures were plotted.
In most cases, the proposed methods achieved close or higher
accuracies than the baseline MUSIC algorithm. Among the
proposed methods, the VBI algorithm achieved the highest
F-measures in most cases, and the number of sources andtheir directions were estimated accurately even though noassumption on the number of sources was required.
V. C ONCLUSION
We proposed sound source localization methods for multiple
source sources. The proposed methods and the baseline MUSIC algorithm were evaluated and compared by the numericalexperiments, and their results indicated that the proposed VBI
algorithm was able to estimate both the number of sources andtheir directions. Expetimental comparison with other stochastic
methods, further evaluation of source localization including
the source distances, and an online extension of the proposed
algorithm are considered as future works.
A PPENDIX
A. Update Rules in Proposed EM Algorithm
The update rules in the proposed EM algorithms are as
below.
h
i

P
(k)
∗
πm,l λm Re ρ(k) + j ωcm F0,m,l f˜m,lm,lhi , (40)n(k) ← P(k)(k) + j ωm F
˜∗
f
0,m,l
m,l πm,l λm Re ρ
m,l
c
h
i
P
(k)
˜H n(k)
π
λ
Re
F
f
m
0,m,l
m,l
m,l m,lρ(k) (= 1/R(k) ) ←, (41)
P
(k)
2
π
λ
|F
|
m
0,m,l
m,l m,lP PK(k) P PK(k)2λm (= 1/σm ) ← 3 l k=1 πm,l / l k=1 πm,l A, (42)
(k)

(k)2

ζm,l (= 1/σ0,m,l ) ← 1/|F0,m,l |2 ,34
Authorized licensed use limited to: NTT. Downloaded on November 07,2023 at 09:28:18 UTC from IEEE Xplore. Restrictions apply.

(43)

2
H
−1
γm,l (= 1/νm,l) ← 1/ym,lWm
ym,l ,
P∞
(k)
(j)
(k)
πm,l ← exp(ϕm,l )/ j=0 exp(ϕm,l ).

(44)

C=2

(0)

H
−2Re[(ρ(k) + jωm /c)F0,m,l f˜m,l]n(k) +∥f˜m,l ∥2 . (48)

B. Update Rules in Proposed VBI Algorithm
The update rules in the proposed VBI algorithms are as
below. Variables without overlines (e.g., ξ (k) for ξ̄ (k) ) denoteparameters for the initial distribution.
ξ̄ (k) ← a/∥a∥, κ̄(k) ← ∥a∥,(49)

← B/C, η̄ ← C,
P PK ′ (k)ᾱλ,m ← 3 l k=1 π̄m,l + αλ,m ,P PK ′ (k)β̄λ,m ← l k=1 π̄m,l D + βλ,m ,

(50)

(k)

(k)

(k)

(k)

(k)

(51)
(52)
(k)

ᾱζ,m,l ← π̄m,l + αζ,m,l , β̄ζ,m,l ← π̄m,l |Y0,m,l |2 + βζ,m,l ,(53)(0)ᾱγ,m,l ← 4π̄m,l + αγ,m,l ,
(54)
H
−1
β̄γ,m,l ← 4π̄m,l ym,lWm
ym,l + βγ,m,l ,
P∞
(k)
(k)
(j)
π̄m,l ← exp(ϕ̄m,l )/ j=0 exp(ϕ̄m,l ),(P(0)π̄ + αn (k = 0)ᾱv(k) ← Pm,l m,l,(k)(k ̸= 0)m,l π̄m,l + 1(βn(k = 0)β̄v(k) ← P.
Pk(j)(k ̸= 0)m,l (1 −j=0 π̄m,l ) + βs
(0)

(55)
(56)
(57)

(58)

(k)

Here, ϕ̄m,l , a, B, C, and D are deﬁned as(0)

ϕ̄m,l = 4(Ψ(ᾱγ,m,l )−log β̄γ,m,l )+Ψ(ᾱv(k) )−Ψ(ᾱv(k) + β̄v(k) )H−1−log det(Wm )−(ᾱγ,m,l /β̄γ,m,l )ym,lWm
ym,l ,

(59)
(k)

(k)

(k)

ϕ̄m,l = 3(Ψ(ᾱλ,m ) − log β̄λ,m ) + Ψ(ᾱζ,m,l ) − log β̄ζ,m,l
(k)

(k)

− (ᾱλ,m /β̄λ,m )D − (ᾱζ,m,l /β̄ζ,m,l )|F0,m,l |2+ Ψ(ᾱv(k) ) + Ψ(β̄v(k) ) − 2Ψ(ᾱv(k) + β̄v(k) ),a=2

X(k) ᾱλ,m

m,l

π̄m,l

β̄λ,m

h
Re

µ̄(k) +j(63)

[1] I. E. Robert, “Robust sound localization: An application of an auditory
perception system for a humanoid robot,” MIT Department of ElectricalEngineering and Computer Science, 1995.
[2] M. Brandstein and D. Ward, Microphone Arrays: Signal Processing
Techniques and Applications, Springer, New York, 2001.
[3] R. Schmidt, “Multiple emitter location and signal parameter estimation,”
IEEE Trans. Antennas Propag., vol. 34, no. 3, pp. 276–280, 1986.
[4] C. Knapp and G. Carter, “The generalized correlation method for
estimation of time delay,” IEEE Trans. Acoust., Speech, Signal Process.,
vol. 24, no. 4, pp. 320–327, 1976.
[5] S. Ando, N. Ono, and T. Nara, “Direct algebraic method for soundsource localization with ﬁnest resolution both in time and frequency,”
in Proc. Int. Congress Sound Vibration, Cairns, July 2007, pp. 3273–
3280.
[6] M. I. Mandel, R. J. Weiss, and D. P. W. Ellis, “Model-based expectationmaximization source separation and localization,” IEEE Trans. Audio,
Speech, Lang. Process., vol. 18, no. 2, pp. 382–394, 2010.
[7] O. Schwartz and S. Gannot, “Speaker tracking using recursive EM
algorithms,” IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 122,
no. 2, pp. 392–402, 2014.
[8] Z. Zohny, S. Naqvi, and J. Chambers, “Variational EM for clustering interaural phase cues in MESSL for blind source separation of speech,” inProc. IEEE Int. Conf. Acoust., Speech, Signal Process., South Brisbane,
April 2015, pp. 3966–3970.
[9] Y. Soussana and S. Gannot, “Variational inference for DOA estimationin reverberant conditions,” in Proc. IEEE European Signal Process.
Conf., A Coruna, September 2019.
[10] S. Adavanne, A. Politis, J. Nikunen, and T. Virtanen, “Sound event
localization and detection of overlapping sources using convolutional
recurrent neural networks,” IEEE J. Sel. Topics Signal Process., vol. 13,
no. 1, pp. 34–48, 2019.
[11] S. Chakrabarty and E. A. P. Habets, “Multi-speaker DOA estimation
using deep convolutional networks trained with noise signals,” IEEE J.
Sel. Topics Signal Process., vol. 13, no. 1, pp. 8–21, 2019.
[12] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihoodfrom imcomplete data via the em algorithm,” J. Royal Statistical Soc.,
B (Methodological), vol. 39, no. 1, pp. 1–38, 1977.
[13] H. Attias, “Inferring parameters and structure of latent variable models
by variational Bayes,” in Proc. the Fifteenth Conf. Uncertainty in
Artiﬁcial Intelligence, Stockholm, July 1999, pp. 21–30.
[14] M. A. Poletti, “Three-dimensional surround sound systems based on
spherical harmonics,” J. Audio Eng. Soc., vol. 53, no. 11, pp. 1004–
1025, 2005.
[15] A. Nehorai and E. Paldi, “Acoustic vector-sensor array processing,”
IEEE Trans. Signal Process., vol. 42, no. 9, pp. 2481–2491, 1994.
[16] D. Levin, E. A. P. Habets, and S. Gannot, “Maximum likelihood
estimation of direction of arrival using an acoustic vector-sensor,” J.
Acoust. Soc. Am., vol. 131, no. 2, pp. 1240–1248, 2012.
[17] R. K. Cook, R. V. Waterhouse, R. D. Berendt, S. Edelman, and M. C. T.
Thompson Jr., “Measurement of correlation coefﬁcients in reverberant
sound ﬁelds,” J. Acoust. Soc. Am., vol. 27, no. 6, pp. 1072–1077, 1955.
[18] H. Kameoka, M. Sato, T. Ono, N. Ono, and S. Sagayama, “Bayesiannonparametric approach to blind separation of inﬁnitely many sparsesources,” IEICE Trans. Fundamentals, vol. E96-A, no. 10, pp. 1928–
1937, 2013.
[19] J. Sethuraman, “A constructive deﬁnition of Dirichlet priors,” Statistica
Sinica, vol. 4, no. 2, pp. 639–650, 1994.
[20] K. V. Mardia and P. E. Jupp, Directional Statistics, John Wiley & Sons,
Chichester, 2009.
[21] K. Takahashi, “SRV-DB,” [Online], http://www.it.cei.uec.ac.jp/SRV-DB,(Accessed on March 7, 2022).

2
A = (ρ(k)2 + ωm/c2 )|F0,m,l |2(k)|F0,m,l |2 + η (k) ,R EFERENCES(k)

ϕm,l = 3 log λm + log ζm,l − λm A − ζm,l |F0,m,l |2 (k ̸= 0),
(47)

µ̄

β̄λ,m

where Ψ(·) is the digamma function.

H−1
ϕm,l = 4 log γm,l −log det(Wm )−γm,l ym,lWm
ym,l , (46)

(k)

π̄m,l

2
D = (µ̄(k)2 + η̄ (k) + ωm/c2 )|F0,m,l |2H(k)]ξ̄ (k) +∥f˜m,l ∥2 , (64)
− 2Re[(µ̄ + jωm /c)F0,m,l f˜m,l

(k)

(k)

(k) ᾱλ,m

m,l

(45)

Here, f˜m,l := [Fx,m,l , Fy,m,l , Fz,m,l ]T , and ϕm,l and A are
deﬁned as

(k)

X

(60)

i
ωm 
∗
F0,m,l f˜m,lc(k) (k)+κ ξ ,(61)hiX(k) ᾱλ,mH
B=2π̄
Re F0,m,l f˜m,lξ̄ (k) + η (k) µ(k) ,m,l m,l β̄λ,m(62)35Authorized licensed use limited to: NTT. Downloaded on November 07,2023 at 09:28:18 UTC from IEEE Xplore. Restrictions apply.