Xiaofei Li

Yutong Ban

Laurent Girin

Alameda-Pineda Xavier

Radu Horaud

X Alameda-Pineda

Online Localization and Tracking of Multiple Moving Speakers in Reverberant Environments

Keywords: Inter-channel acoustic features, reverberant environments, sound-source localization, multiple target tracking, speaker tracking, Bayesian variational inference, expectationmaximization

We address the problem of online localization and tracking of multiple moving speakers in reverberant environments. The paper has the following contributions. We use the direct-path relative transfer function (DP-RTF), an interchannel feature that encodes acoustic information robust against reverberation, and we propose an online algorithm well suited for estimating DP-RTFs associated with moving audio sources. Another crucial ingredient of the proposed method is its ability to properly assign DP-RTFs to audio-source directions. Towards this goal, we adopt a maximum-likelihood formulation and we propose to use exponentiated gradient (EG) to efficiently update source-direction estimates starting from their currently available values. The problem of multiple speaker tracking is computationally intractable because the number of possible associations between observed source directions and physical speakers grows exponentially with time. We adopt a Bayesian framework and we propose a variational approximation of the posterior filtering distribution associated with multiple speaker tracking, as well as an efficient variational expectation maximization (VEM) solver. The proposed online localization and tracking method is thoroughly evaluated using two datasets that contain recordings performed in real environments.

I. INTRODUCTION

The localization and tracking of multiple speakers in real world environments are very challenging tasks, in particular in the presence of reverberation and ambient noise and of natural conversations, e.g. short sentences, speech pauses and frequent speech turns among speakers. Methods based on time differences of arrival (TDOAs) between microphones, such as generalized cross-correlation [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF], are typically used for single-speaker localization, e.g. [START_REF] Chen | Time delay estimation in room acoustic environments: an overview[END_REF]. In the case of multiple speakers, beamforming-based methods, e.g. steered-response power (SRP) [START_REF] Dibiase | Robust localization in reverberant rooms[END_REF], and subspace methods, e.g. multiple signal classification (MUSIC) [START_REF] Ishi | Evaluation of a music-based real-time sound localization of multiple sound sources in real noisy environments[END_REF], are widely used. The W-disjoint orthogonality (WDO) principle [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] assumes that the audio signal is dominated by a single audio source in small regions of the time-frequency (TF) domain. This assumption is particularly valid in the case of speech signals. Applying the short-time Fourier transform (STFT), or any other TF representation, inter-channel localization features, such as the interaural phase differences (IPDs) [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], can be extracted. In [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], multiple-speaker localization is based on the histogram of inter-channel features, which is suitable only in the case where there is no wrapping of phase measures. In [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF], a Gaussian mixture model (GMM) is used as a generative model of the inter-channel features of multiple speakers, with each GMM representing one speaker, and each GMM component representing one candidate inter-channel time delay. An expectation maximization (EM) algorithm iteratively estimates the component weights and assigns the features to their corresponding candidate time delays. This method overcomes the phase ambiguity problem by jointly considering all frequencies in the likelihood maximization procedure. After maximizing the likelihood, the azimuth of each speaker is given by the component that has the highest weight in the corresponding GMM. The complex-valued version of IPD, i.e. the pair-wise relative phase ratio (PRP), is used in [START_REF] Dorfan | Tree-based recursive expectationmaximization algorithm for localization of acoustic sources[END_REF]. Instead of setting one GMM for each speaker, a single complex Gaussian mixture model (CGMM) is used for all speakers with each component representing one candidate speaker location. After maximizing the likelihood of the PRP features, with an EM algorithm, the weight of each component represents the probability that there is an active speaker at the corresponding candidate location. Therefore, for an unknown number of speakers, counting and localization of active speakers can be jointly carried out by selecting components with large weights.

The inter-channel features and associated localization methods mentioned above assume a direct-path propagation model: hence, they perform poorly in reverberant environments. To overcome this limitation, several TDOA estimators based on system identification were proposed in [START_REF] Huang | Adaptive multichannel time delay estimation based on blind system identification for acoustic source localization[END_REF]- [START_REF] Kowalczyk | Blind system identification using sparse learning for TDOA estimation of room reflections[END_REF]. In [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] it is proposed to use the DP-RTF as a TF-domain interchannel localization feature robust against reverberation. The estimation of the DP-RTF is based on the identification of the room impulse response (RIR) in the STFT-domain, i.e. the convolutive transfer function (CTF) [START_REF] Avargel | System identification in the short-time Fourier transform domain with crossband filtering[END_REF], [START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF]. Overall, the method of [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] combines the merits of robust TDOA estimators [START_REF] Huang | Adaptive multichannel time delay estimation based on blind system identification for acoustic source localization[END_REF]- [START_REF] Kowalczyk | Blind system identification using sparse learning for TDOA estimation of room reflections[END_REF] and of the WDO assumption mentioned above.

To localize moving speakers, one-stage methods such as SRP and MUSIC can be directly used using frame-wise spatial spectrum estimators. In contrast, methods based on inter-channel features require to assign frame-wise features to speakers in an adaptive/recursive way, e.g. the smoothed histogram method of [START_REF] Pavlidi | Real-time multiple sound source localization and counting using a circular microphone array[END_REF]. Similar to [START_REF] Dorfan | Tree-based recursive expectationmaximization algorithm for localization of acoustic sources[END_REF], [START_REF] Schwartz | Speaker tracking using recursive EM algorithms[END_REF] uses one CGMM for each predefined speaker; the model is plugged into a recursive EM (REM) algorithm in order to update the Input: STFT of microphone signals:

{x i t,f } I,F -1 i=1,f =0
Eqs. (3), [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] build the vectors: {x m t,f , y m t,f } M,F -1 m=1,f =0

Alg. 1 computes CTFs and inverse covariances:

{ã t,f } F -1 f =0 , { Pt,f } F -1 f =0
DP-RTF estimation and consistency test:

C t = {{ĉ i t,f } i∈I f } F -1 f =0
Eqs. [START_REF] Roman | Binaural tracking of multiple moving sources[END_REF], [START_REF] Evers | Bearing-only acoustic tracking of moving speakers for robot audition[END_REF] Predicted features:

{c i,d f } F -1,I,D f =0,i=2,d=1 {ã t-1,f } F -1 f =0 { Pt-1,f } F -1 f =0 {w t-1,d } D d=1 {µ t-1,n } N n=1 {Γ t-1,n } N n=1
Online DP-RTF Estimation (Section II)

Multiple Speaker Localization (Section III)

Multiple Speaker Tracking (Section IV)

Candidate directions

Fig. 1: Flowchart of the proposed multiple-speaker localization and tracking methodology.

mixture's weights.

Speaker tracking methods are generally based on Bayesian inference which combines localization with dynamic models in order to estimate the posterior probability distribution of audio-source directions, e.g. [START_REF] Roman | Binaural tracking of multiple moving sources[END_REF]- [START_REF] Ban | Exploiting the complementarity of audio and visual data in multi-speaker tracking[END_REF]. Kalman filtering and particle filtering were used in [START_REF] Liang | Robust tracking of moving sound source using multiple model kalman filter[END_REF] and in [START_REF] Vermaak | Nonlinear filtering for speaker tracking in noisy and reverberant environments[END_REF], respectively, for tracking a single audio source. In order to address the problem of multiple speakers, possibly with unknown and time-varying number of speakers, additional discrete latent variables are needed, i.e. observation-to-speaker assignments, as well as speaker-birth and -death processes, e.g. [START_REF] Ba | An online variational bayesian model for multi-person tracking from cluttered scenes[END_REF], [START_REF] Gebru | Audio-visual speaker diarization based on spatiotemporal Bayesian fusion[END_REF]. Sampling-based methods were widely used, e.g. extended particle filtering [START_REF] Fallon | Acoustic source localization and tracking of a time-varying number of speakers[END_REF]- [START_REF] Cevher | Acoustic multitarget tracking using direction-of-arrival batches[END_REF], or sequential Monte Carlo implementation of the probability hypothesis density (PHD) filter [START_REF] Vo | Tracking multiple speakers using random sets[END_REF], [START_REF] Ma | Tracking an unknown time-varying number of speakers using tdoa measurements: A random finite set approach[END_REF]. However, the computational burden of samplingbased methods can be prohibitive in practice. Under some assumptions, the multiple-target tracking GMM-PHD filter of [START_REF] Vo | The gaussian mixture probability hypothesis density filter[END_REF] has an analytical solution and is computationally efficient: it was adopted for multiple-speaker tracking in [START_REF] Evers | Bearing-only acoustic tracking of moving speakers for robot audition[END_REF].

In this paper we propose a method for the simultaneous localization and tracking of multiple moving speakers (please refer to Figure 1 for a method overview). The paper has the following original contributions:

• Since we deal with moving speakers or, more generally, with moving audio sources, DP-RTF features are computed using the online CTF estimation framework presented in [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF], based on recursive least squares (RLS), rather than using the batch CTF estimation of [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] which assumes static audio sources. The online RLS algorithm has a faster convergence rate than the least mean squares (LMS) algorithms described in [START_REF] Huang | Adaptive multichannel time delay estimation based on blind system identification for acoustic source localization[END_REF], [START_REF] Doclo | Robust adaptive time delay estimation for speaker localization in noisy and reverberant acoustic environments[END_REF]. This is important when dealing with multiple moving sources, where the adaptive estimator is required to quickly switch between multiple sources and to deal with moving sources.

• A crucial ingredient of multiple speaker localization is to properly assign acoustic features, i.e. DP-RTFs, to audiosource directions. We adopt the maximum-likelihood formulation of [START_REF] Dorfan | Tree-based recursive expectationmaximization algorithm for localization of acoustic sources[END_REF]. We propose to use EG [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF] to update the source directions from their current estimated values.

The EG-based recursive estimator proposed below is better suited for moving sources/speakers than the batch estimator proposed in [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF]. • The problem of multiple speaker tracking is computationally intractable because the number of possible associations between acoustic features and sources/speakers grows exponentially with time. In this paper we adopt a Bayesian variational approximation of the posterior filtering distribution which leads to an efficient VEM algorithm. In order to deal with a varying number of speakers, we propose a birth process which allows to initialize new speakers at any time.

This paper is an extended version of [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF] which has proposed an online DP-RTF method that has been combined with REM to estimate the source directions. In this paper, while we keep the DP-RTF method of [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF] we propose to use EG. The advantages of using EG instead of REM are described in detail in Section III. Moreover, the multiple speaker tracking method is completely novel.

The paper is organized as follows (please refer to Figure 1). Section II presents the online DP-RTF estimation method. Section III describes the EG-based speaker localization method and Section IV describes the variational approximation of the tracker and the associated VEM algorithm. Section V presents an empirical evaluation of the method based on experiments performed with real audio recordings. Section VI concludes the paper. Supplemental materials are available on our website. 1

II. RECURSIVE MULTICHANNEL DP-RTF ESTIMATION

A. Recursive Least Squares

For the sake of clarity, we first consider the noise-free single-speaker case. In the time domain x i (τ ) = a i (τ ) s(τ ) is the i-th microphone signal, i = 1, . . . , I, where τ is the time index, s(τ ) is the source signal, a i (τ ) is the RIR from the source to the i-th microphone, and denotes the convolution. Applying the STFT and using the CTF approximation, for each frequency index f = 0, . . . , F -1 we have:

x i t,f = a i t,f s t,f = Q-1 q=0 a i q,f s t-q,f , (1)
where x i t,f and s t,f are the STFT coefficients of the corresponding signals, and the CTF a i t,f is a sub-band representation of a i (τ ). Here, the convolution is executed with respect to the frame index t. The number of CTF coefficients Q is related to the reverberation time of the RIR. The first CTF coefficient a i 0,f mainly consists of the direct-path information, 1 https://team.inria.fr/perception/research/multi-speaker-tracking/ thence the DP-RTF is defined as the ratio between the first CTF coefficients of two channels: a i 0,f /a r 0,f , where channel r is the reference channel.

Based on the cross-relation method [START_REF] Xu | A least-squares approach to blind channel identification[END_REF], using the CTF model of one microphone pair (i, j), we have:

x i t,f a j t,f = x j t,f a i t,f
. This can be written in vector form as:

x i t,f a j f = x j t,f a i f , (2)
with a i f = (a i 0,f , . . . , a i Q-1,f ) , where denotes matrix/vector transpose, and x i t,f = (x i t,f , . . . , x i t-Q+1,f ) . The CTF vector involving all channels is defined as a f = (a 1 f , . . . , a I f ) . There is a total of I(I -1)/2 distinct microphone pairs, indexed by (i, j) with i = 1, . . . , I -1 and j = i + 1, . . . , I. For each pair, we construct a cross-relation equation in terms of a f . For this aim, we define:

x ij t,f = [0, . . . , 0 (i-1)Q , x j t,f , 0, . . . , 0 (j-i-1)Q , -x i t,f , 0, . . . , 0 (I-j)Q ] . (3)
Then, for each pair (i, j), we have:

x ij t,f a f = 0. (4)
Let's assume, for simplicity, that the reference channel is r = 1. To avoid the trivial solution a f = 0 of (4), we constrain the first CTF coefficient of the reference channel to be equal to 1. This is done by dividing both sides of (4) by a 1 0,f and by moving the first entry of x ij t,f , denoted by -y ij t,f , to the right side of (4), which rewrites as:

xij t,f ãf = y ij t,f , (5)
where xij t,f is x ij t,f with the first entry removed, and ãf is the relative CTF vector:

ãf = ã1 f a 1 0,f , a 2 f a 1 0,f , . . . , a I f a 1 0,f , (6)
where ã1 f = (a 1 1,f , . . . , a 1 Q-1,f ) denotes a 1 f with the first entry removed. For i = 2, . . . , I, the DP-RTFs appear in [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF] as the first entries of

a i f a 1 0,f
. Therefore, the DP-RTF estimation amounts to solving [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. Equation ( 5) is defined for one microphone pair and for one frame. In batch mode, the terms xij t,f and y ij t,f of this equation can be concatenated accross microphone pairs and frames to construct a least square formulation. For online estimation, we would like to update the ãf using the current frame t. For notational convenience, let m = 1, . . . , M denote the index of a microphone pair, where M = I(I -1)/2. Then let the superscript ij be replaced with m . The fitting error of ( 5) is

e m t,f = y m t,f -xm t,f ãf . (7)
At the current frame t, for the microphone pair m, RLS aims to minimize the error

J m t,f = t-1 t =1 M m =1 λ t-t |e m t ,f | 2 + m m =1 |e m t,f | 2 , ( 8
)
Algorithm 1 RLS at frame t

Input: xm t,f , y m t,f , m = 1, . . . , M Initialization: ã0 t,f ← ãM t-1,f , P 0 t,f ← λ -1 P M t-1,f for m = 1 to M do e m t,f = y m t,f -xm t,f ãm-1 t,f g = P m-1 t,f xm * t,f /(1 + xm t,f P m-1 t,f x m * t,f ) P m t,f = P m-1 t,f -gx m t,f P m-1 t,f ãm t,f = ãm-1 t,f + e m t,f g end for Output: ãM t,f , P M t,f
which sums up the fitting error of all the microphone pairs for the past frames and the microphone pairs up to m for the current frame. The forgetting factor λ ∈ (0, 1] gives a lower weight to older frames, whereas all microphone pairs have the same weight at each frame. To minimize J m t,f , we set its complex derivative with respect to ã * f to zero, where * denotes complex conjugate. We obtain an estimate of ãf at frame t for microphone pair m as:

ãm t,f = R m -1 t,f r m t,f , (9)
with

R m t,f = t-1 t =1 M m =1 λ t-t xm * t ,f xm t ,f + m m =1 xm * t,f xm t,f , r m t,f = t-1 t =1 M m =1 λ t-t xm * t ,f y m t ,f + m m =1 xm * t,f y m t,f .
It can be seen that the covariance matrix R m t,f is computed based on the rank-one modification, thence its inverse, denoted by P m t,f , can be computed using the Sherman-Morrison formula, without the need of matrix inverse. The recursion procedure is summarized in Algorithm 1, where g is the gain vector. The current frame t is initialized with the previous frame t -1. At the first frame, we initialize ã0 1,f as zero, and P 0 1,f as the identity. At each frame, all microphone pairs are related to the same CTF vector that corresponds to the current speaker direction, hence all microphone pairs should be simultaneously used to estimate the CTF vector of the current frame. In batch mode, this can be easily implemented by concatenating the microphone pairs. However, in RLS, to satisfy the rank-one modification of the covariance matrix, we need to process the microphone pairs one by one as shown in [START_REF] Huang | Adaptive multichannel time delay estimation based on blind system identification for acoustic source localization[END_REF] and Algorithm 1. At the end of the iterations over all microphone pairs, ãM t,f is the "final" CTF estimation for the current frame, and is used for speaker localization. The DP-RTF estimates, denoted as ci t,f , i = 2, . . . , I, are obtained from ãM t,f . Note that implicitly we have c1 t,f = 1.

B. Multiple Moving Speakers

So far, the proposed online DP-RTF estimation method has been presented in the noise-free single-speaker case. The noisy multiple-speaker case was considered in [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF], but only for static speakers, i.e. batch mode, and in the two-channel case.

We summarize the principles of this method and then explain in details the present online/multi-channel extension.

1) Estimation of the CTF vector: It is reasonable to assume that the CTF vector doesn't vary over a few consecutive frames and that only one speaker is active within a small region in the TF domain, due to the sparse representation of speech in this domain. Consequently, the CTF vector can be estimated over the current frame and a few past frames. An estimated CTF value, at each TF bin, is then assumed to be associated with only one speaker. The CTF vector computation in the case of multiple speakers can be carried out using the RLS algorithm, presented in Section II-A, by adjusting the forgetting factor λ to yield a short memory.

The forgetting factor λ is set to λ = P -1 P +1 , where P is the number of frames being used. To efficiently estimate the CTF vector ãM t,f of length IQ -1, we need ρ × (IQ -1) equations, where the parameter ρ should be chosen in such a way to achieve a good tradeoff between the validity of the above assumptions and robust estimation of ãM t,f . To guarantee that ρ × (IQ -1) equations are available, we need P = ρ(IQ-1) I(I-1)/2 ≈ ρ 2Q I-1 frames. One may observe that the number of frames needed to estimate ãM t,f decreases as the number of microphones increases.

2) Noise reduction: When noise is present, especially if the noise sources are temporally/spatially correlated, the CTF estimate can be contaminated. In addition, even in a lownoise case, many TF bins are dominated by noise due to the sparsity of speech spectra. To classify the speech frames and noise frames, and to remove the noise, we use the inter-frame spectral subtraction algorithm proposed in [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF], [START_REF] Li | Estimation of relative transfer function in the presence of stationary noise based on segmental power spectral density matrix subtraction[END_REF].

The cross-and auto-power spectral density (PSD) between the convolution vector of the microphone signals, i.e. x i t,f , and the current frame of the reference channel, i.e. x 1 t,f , is computed by averaging the cross-and auto-periodograms over frames. In the present work, we use recursive averaging:

φ i t,f = βφ i t-1,f + (1 -β)x i t,f x 1 * t,f , i = 1, . . . , I, (10)
where the smoothing parameter β is set to achieve a good tradeoff between low noise PSD variance and fast tracking of speech variation. The noise frames and speech frames are classified based on the minimum statistics [START_REF] Li | Estimation of relative transfer function in the presence of stationary noise based on segmental power spectral density matrix subtraction[END_REF] of the PSD of x 1 t,f , i.e. the first entry of φ 1 t,f . If the frames are well classified then the noise frames only include negligible speech power, due to the sparsity and non-stationarity of speech; the speech frames include noise power similar to the noise frames, due to the stationarity of noise. Therefore, inter-frame spectral subtraction can be performed as follows: for each speech frame, the cross-and auto-PSD of its nearest noise frame is subtracted from its cross-and auto-PSD, then its noise-free cross-and auto-PSD is obtained and denoted as φi t,f .

Instead of using x i t,f , we use φi t,f to construct (3). Correspondingly, we have a new formula (4), which is still valid, since it is equivalent to, with noise removed, taking the crossand auto-PSD between both sides of the initial formula (4) and x 1 t,f . In the RLS process, only the speech frames (after spectral subtraction) are used, and the noise frames are skipped. A speech frame with a preceding noise frame is initialized with the latest speech frame.

3) Consistency test: In practice, a DP-RTF estimate can sometimes be unreliable. Possible reasons are that in a small frame region, (i) the CTF is time-varying due to a fast movement of the speakers, (ii) multiple speakers are present, (iii) only noise is present due to a wrong noise-speech classification, or (iv) only reverberation is present at the end of a speech occurrence. In [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF], a consistency test was proposed to tackle this problem: If a small frame region indeed corresponds to one active speaker, the DP-RTFs estimated using different reference channels are consistent, otherwise the DP-RTFs are biased, with inconsistent bias values. In the present work, we use the first and second channels as references, we obtain the DP-RTF estimates ci t,f (with c1 t,f = 1) and ci t,f (with c2 t,f = 1), respectively. Then ci t,f and ci t,f /c 1 t,f are two estimates of the same DP-RTF a i 0,f /a 1 0,f . To measure the similarity between these two estimates, we define the vectors

c i 1,t,f = (1, ci t,f ) and c i 2,t,f = (1, ci t,f /c 1 t,f )
, where the first entries are the DP-RTFs corresponding to a 1 0,f /a 1 0,f = 1. The similarity is the cosine of the angle between the two unit vectors:

d i t,f = |c i H 1,t,f c i 2,t,f | c i H 1,t,f c i 1,t,f c i H 2,t,f c i 2,t,f , (11)
where H denotes conjugate transpose. If

d i t,f ∈ [0, 1]
is larger than a threshold (which is fixed to 0.75 in this work) then the two estimates are consistent, otherwise they are simply ignored. Then, the two estimates are averaged and normalized as done in [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF], resulting in a final complex-valued feature ĉi t,f whose module lies in the interval [0, 1]. Finally, at frame t, we obtain a set of features

C t = {{ĉ i t,f } i∈I f } F -1 f =0
, where I f ⊆ {2, . . . , I} denotes the set of microphone indices that pass the consistency test. Note that I f is empty if frame t is a noise frame at frequency f , or if no channel passes the consistency test. Each one of these features is assumed to be associated with only one speaker.

III. LOCALIZATION OF MULTIPLE MOVING SPEAKERS

In this section we describe the proposed frame-wise online multiple-speaker localizer. We start by briefly presenting the underlying complex Gaussian mixture model, followed by the recursive estimation of its parameters.

A. Generative Model for Multiple-Speaker Localization

In order to associate DP-RTF features from C t with speakers and to localize each active speaker, we adopt the generative model proposed in [START_REF] Dorfan | Tree-based recursive expectationmaximization algorithm for localization of acoustic sources[END_REF]. Let D = { θ1 , . . . , θd , . . . , θD } be a set of D candidate source directions, e.g. azimuth angles. An observed feature ĉi t,f (cf. Section II), when emitted by a sound source located along the direction θd , is assumed to be drawn from a complex-Gaussian distribution with mean

c i,d f and variance σ 2 , i.e. ĉi t,f |d ∼ N c (c i,d f , σ 2 ). The mean c i,d f
is the predicted feature at frequency f for channel i, and is precomputed based on direct-path propagation along azimuth θd to the microphones. The variance σ 2 is empirically set as a constant value. The marginal density of an observed feature ĉi t,f (taking into account all candidate directions) is a CGMM with each component corresponding to a candidate direction:

p(ĉ i t,f |D) = D d=1 w d N c (ĉ i t,f ; c i,d f , σ 2 ), (12)
where w d ≥ 0 is the prior probability (component weight) of the d-th component, with

D d=1 w d = 1.
Let us denote the vector of weights with w = (w 1 , . . . , w D ) . Note that this vector is the only free parameter of the model.

Assuming that the observations in C t are independent, the corresponding (normalized) negative log-likelihood function, as a function of w d , is given by:

L t = - 1 |C t | ĉi t,f ∈Ct log D d=1 w d N c (ĉ i t,f ; c i,d f , σ 2 ) , (13)
where |C t | denotes the cardinality of C t . Once L t is minimized, each weight w d represents the probability that a speaker is active in the direction θd . Therefore, sound source localization amounts to the minimization of L t . In addition, taking into account the fact that the number of actual active speakers is much lower than the number of candidate directions, an entropy term was proposed in [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] as a regularizer to impose a sparse solution for w d . The entropy is defined as

H = - D d=1 w d log(w d ). (14)
The concave-convex procedure [START_REF] Yuille | The concave-convex procedure[END_REF] was adopted in [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF], to minimize the objective function L + γH w.r.t. w, where L is the normalized negative log-likelihood of the DP-RTF features of all frames, i.e. batch mode optimization, and the positive scalar γ was used to control the tradeoff between likelihood minimization and imposing sparsity over the weights. In the batch mode, the weight vector w is shared across all frames. Hence this method is not suitable for moving speakers.

B. Recursive Parameter Estimation

We now describe a recursive method for updating the weight vector from w t-1 to w t , i.e. from frame t -1 to frame t, using the DP-RTF features at t. This can be formulated as the following online optimization problem [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF]:

w t = argmin w χ(w, w t-1 ) + η(L t + γH), (15) s
w d = 1, (16)
where χ(a, b) is a distance between a and b. The positive scalar factor η controls the parameter update rate. To minimize [START_REF] Pavlidi | Real-time multiple sound source localization and counting using a circular microphone array[END_REF], the derivative of the objective function w.r.t w is set to zero, yielding a set of equations with no closed-form solution. To speed up the computation, it is assumed that w t is close to w t-1 , thence the derivative of L t + γH at w can be approximated with the derivative of L t + γH at w t-1 .

This assumption is reasonable when parameter evolution is not too fast. As a result, when the distance χ(w, w t-1 ) is Euclidean, the objective function leads to gradient descent with a step length equal to η. Nevertheless, the constraints (16) lead to an inefficient gradient descent procedure. To obtain an efficient solver, we exploit the fact that the weights w d are probability masses, hence we replace the Euclidean distance with the more suitable Kullback-Leibler divergence, i.e. χ(w, w t-1 ) = D d=1 w d log w d w t-1,d , which results in the exponentiated gradient algorithm [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF].

The partial derivatives of L t and H w.r.t w d at the point w t-1,d are computed with, respectively:

∂(L t ) ∂w d w t-1,d = - 1 |C t | ĉi t,f ∈Ct N c (ĉ i t,f ; c i,d f , σ 2 ) D d =1 w t-1,d N c (ĉ i t,f ; c i,d f , σ 2 ) , ∂H ∂w d w t-1,d = -(1 + log(w t-1,d )), ∀d ∈ {1 . . . D}. (17)
Then, the exponentiated gradient,

r t-1,d = e -η ∂(-L t ) ∂w d w t-1,d +γ ∂H ∂w d w t-1,d , ∀d ∈ {1 . . . D}, (18)
is used to update the weights with:

w t,d = r t-1,d w t-1,d D d =1 r t-1,d w t-1,d , ∀d ∈ {1 . . . D}. ( 19
)
It is clear from [START_REF] Ban | Exploiting the complementarity of audio and visual data in multi-speaker tracking[END_REF] that the parameter constraints ( 16) are necessarily satisfied. The exponentiated gradient algorithm sequentially evaluates ( 17), ( 18) and ( 19) at each frame. At the first frame, the weights are initialized with the uniform distribution, namely w 1,d = 1 D . When C t is empty, such as during a silent period, the parameters are recursively updated with

w t,d = (1 -η )w t-1,d + η 1 D .
The weight w t as a function of θd , i.e. w t,d , exhibits a handful of peaks that could correspond to active speakers. The use of an entropy regularization term was shown to both suppress small spurious peaks, present without using the regularization term, and to sharpen the peaks corresponding to actual active speakers, thus allowing to better localize true speakers and to eliminate erroneous ones. In the case of moving speakers, a peak should shift along time from a direction θd to a nearby direction. Spatial smoothing of the weight function raises the weight values around a peak, which results in smoother peak jumps. In our experiments, spatial smoothing is carried out with w t,d = (w t,d +0.02w t,d-1 +0.02w t,d+1 )/1.04, where the smoothing factor 0.02 is empirically chosen in order to smooth peak motion from one frame to the next, while avoiding the peaks to collapse. One may think that spatial smoothing and entropy regularization neutralize each other, but in practice it was found that their combination is beneficial.

C. Peak Selection and Frame-wise Speaker Localization

Frame-wise localization and counting of active speakers could be carried out by selecting the peaks of w t ( θd ) larger than a predefined threshold [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF], [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF]. However, peak selection does not exploit the temporal dependencies of moving speakers. Moreover, peak selection can be a risky process since a too high or too low threshold value may lead to undesirable missed detection or false alarm rates. In order to avoid these problems, we adopt a weighted-data Bayesian framework: all the candidate directions and the associated weights are used as observations by the multiple speaker tracking method described in Section IV below. The localization results obtained with peak selection are compared with the localization results obtained with the proposed tracker in Section V.

IV. MULTIPLE SPEAKER TRACKING

Let N be the maximum number of speakers that can be simultaneously active at any time t, and let n be the speaker index. Moreover, let n = 0 denote no speaker. We now introduce the main variables and their notations. Upper case letters denote random variables while lower case letters denote their realizations.

Let S tn be a latent (or state) variable associated with speaker n at frame t, and let S t = (S t1 , . . . , S tn , . . . , S tN ). S tn is composed of two parts: the speaker direction and the speaker velocity. In this work, speaker direction is defined by an azimuth θ tn . To avoid phase (circular) ambiguity we describe the direction with the unit vector U tn = (cos(θ tn ), sin(θ tn )) . Moreover, let V tn ∈ R be the angular velocity. Altogether we define a realization of the state variable as s tn = [u tn ; v tn ] where the notation [•; •] stands for vertical vector concatenation.

Let O t = (O t1 , . . . , O td , . . . , O tD ) be the observed variables at frame t. Each realization o td of O td is composed of a candidate location, or azimuth θtd ∈ D, and a weight w td . The weight w td is the probability that there is an active speaker in the direction θtd , namely [START_REF] Pavlidi | Real-time multiple sound source localization and counting using a circular microphone array[END_REF]. As above, let the azimuth be described by a unit vector b td = (cos( θtd ), sin( θtd )) . In summary we have o td = [b td ; w td ]. Moreover, let Z td be a (latent) assignment variable associated with each observed variable O td , such that Z td = n means that the observation indexed by d at frame t is assigned to active speaker n ∈ {0, . . . , N }. Note that Z td = 0 is a "fake" assignmentthe corresponding observation is assigned to an audio source that is either background noise or any other source that has not yet been identified as an active speaker.

The problem at hand can now be cast into the estimation of the filtering distribution p(s t , z t |o 1:t ), and further inference of s t and z t . In this work we make two hypotheses, namely (i) that the observations at frame t only depend on the assignment and state variables at t, and (ii) that the prior distribution of the assignment variables is independent of all the other variables. By applying the Bayes rule together with these hypotheses, and ignoring terms that do not depend on s t and z t , the filtering distribution is proportional to:

p(s t , z t |o 1:t ) ∝ p(o t |s t , z t )p(z t )p(s t |o 1:t-1 ), (20)
which contains three terms: the observation model, the prior distribution of the assignment variable and the predictive distribution over the sources state. We now characterize each one of these three terms.

1) Audio observation model: The audio observation model describes the distribution of the observations given speakers state and assignment. We assume the different observations are independent, conditioned on speakers state and assignment, which can be written as:

p(o t |z t , s t ) = D d=1 p(o td |z t , s t ). ( 21
)
Since the weights describe the confidence associated with each observed azimuth, we adopt the weighted-data GMM model of [START_REF] Gebru | Em algorithms for weighted-data clustering with application to audio-visual scene analysis[END_REF]:

p(b td |Z td = n, s tn ; w td ) = N (b td ; Ms tn , 1 w td Σ) if n ∈ {1, . . . , N } U(vol(G)) if n = 0 , (22)
where the matrix M = [I 2×2 , 0 2×1 ] projects the state variable onto the space of source directions and Σ is a covariance matrix (set empirically to a fixed value in the present study).

Note that the weight plays the role of a precision: The higher the weight w td , the more reliable the source direction b td . The case Z td = 0 follows a uniform distribution over the volume of the observation space.

2) Prior distribution: The prior distribution of the assignment variable is independent over observations and is assumed to be uniformly distributed over all the speakers (including the case n = 0), hence:

p(z t ) = D d=1 p(Z td = n) with π dn = p(Z td = n) = 1 N + 1 . (23)
3) Predictive distribution: The predictive distribution describes the relationship between the state s t and the past observations up to frame t, o 1:t-1 . To calculate this distribution, we first marginalize p(s t |o 1:t-1 ) over s t-1 , writing:

p(s t |o 1:t-1 ) = p(s t |s t-1 )p(s t-1 |o 1:t-1 )ds t-1 , (24)
where the two terms under the integral are the state dynamics and the marginal filtering distribution of the state variable at frame t -1, respectively. We model the state dynamics as a linear-Gaussian first-order Markov process, independent over the speakers, i.e. :

p(s t |s t-1 ) = N n=1 N (s tn ; D t-1,n s t-1,n , Λ tn ), (25)
where Λ tn is the dynamics' covariance matrix and D t-1,n is the state transition matrix. Given the estimated azimuth θ t-1,n and angular velocity v t-1,n at frame t -1, we have the following relation:

cos(θ tn ) sin(θ tn ) = cos(θ t-1,n + v t-1,n ∆t) sin(θ t-1,n + v t-1,n ∆t) , ( 26
)
where ∆t is the time increment between two consecutive frames. Expanding [START_REF] Cevher | Acoustic multitarget tracking using direction-of-arrival batches[END_REF] and assuming that the angular displacement v t-1,n ∆t is small, the state transition matrix can be written as:

D t-1,n =   1 0 -sin(θ t-1,n )∆t 0 1 cos(θ t-1,n )∆t 0 0 1   . (27)
In the following D t-1,n is written as D, only to lighten the equations.

A. Variational Expectation Maximization Algorithm

At this point, the standard solution to the calculation of the filtering distribution consists of using EM methodology. EM alternates between evaluating the expected complete-data loglikelihood and maximizing this expectation with respect to the model parameters. More precisely, the expectation writes:

J(Θ, Θ o ) = E p(zt,st|o1:t,Θ o ) [log p(z t , s t , o 1:t |Θ)] , (28)
where Θ o denotes the current parameter estimates and Θ denotes the new estimates, obtained via maximization of [START_REF] Ma | Tracking an unknown time-varying number of speakers using tdoa measurements: A random finite set approach[END_REF]. However, given the hybrid combinatorial-and-continuous nature of the latent space, such solution is intractable in practice, due to combinatorial explosion. We thus propose to use of a variational approximation to solve the problem efficiently. We inspire from [START_REF] Ba | An online variational bayesian model for multi-person tracking from cluttered scenes[END_REF] and propose the following factorization:

p(z t , s t |o 1:t ) ≈ q(z t , s t ) = q(z t ) N n=0 q(s tn ). (29)
The optimal solution is then given by two E-steps, an E-S step for each individual state variable S tn and an E-Z step for the assignment variable Z t :

log q(s tn ) = E q(zt) m =n q(stm) [log p(z t , s t |o 1:t )], (30)
log q(z t ) = E q(st) [log p(z t , s t |o 1:t )]. (31)
It is easy to see that in order to compute [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF] and [START_REF] Kivinen | Exponentiated gradient versus gradient descent for linear predictors[END_REF], two elements are needed: the predictive distribution ( 24) and the marginal filtering distribution at t -1, p(s t-1 |o 1:t-1 ).

Remarkably, as a consequence of the factorization (29), we can replace p(s t-1 |o 1:t-1 ) with q(s t-1 ) = N n=1 q(s t-1,n ) in [START_REF] Fallon | Acoustic source localization and tracking of a time-varying number of speakers[END_REF] and compute the predictive distribution as follows:

p(s t |o 1:t-1 ) ≈ p(s t |s t-1 ) N n=1 q(s t-1,n )ds t-1 . (32)
This predictive distribution factorizes across speakers. Moreover, one prominent feature of the proposed variational approximation is that, if the posterior distribution at time t -1 q(s t-1,n ) is assumed to be a Gaussian, say

q(s t-1,n ) = N (s t-1,n ; µ t-1,n , Γ t-1,n ), (33)
then (the approximation of) the predictive distribution ( 32) is a Gaussian. More specifically, the derivation of (32) leads to:

p(s tn |o 1:t-1 ) = N (s tn ; Dµ t-1,n , DΓ t-1,n D + Λ tn ). (34)
Moreover, as we will see in the E-S-step below, the posterior distribution at time t, q(s tn ), is also a Gaussian.

1) E-S step: The computation of the variational posterior distribution q(s tn ), for all currently tracked speakers, is carried out by developing [START_REF] Li | Online localization of multiple moving speakers in reverberant environments[END_REF] as follows. We first exploit ( 20), ( 21), ( 23) and [START_REF] Yuille | The concave-convex procedure[END_REF] to rewrite log p(z t , s t |o 1:t ) in (30) as a sum of individual log-probabilities. Then we eliminate all terms not depending on s tn and we evaluate the expectation of the remaining terms. Because the terms not depending on s tn were disregarded, the expectation is computed only with respect to q(z t ). This nicely makes the computation of q(s tn ) independent of the structure of q(s tm ) for m = n. Eventually, this yields a Gaussian distribution:

q(s tn ) = N (s tn ; µ tn , Γ tn ), (35)
with the following parameters:

Γ tn = D d=1 α tdn w td M Σ -1 M + Λ tn + DΓ t-1,n D -1 -1 , (36)
µ tn = Γ tn M Σ -1 D d=1 α tdn w td b td + Λ tn + DΓ t-1,n D -1 Dµ t-1,n , (37)
where α tdn = q(Z td = n) is the variational posterior distribution of the assignment variable, which will be detailed in Section IV-A2. Importantly, the first two entries of µ tn in (37) represent the estimated azimuth of speaker n. The "final" azimuth estimate at frame t is thus given by this subvector at the end of the VEM iterations. Since we use a unit-vector representation, we normalize this vector at each iteration of the algorithm. Finally, note that since we have shown that q(s t-1,n ) being Gaussian leads to q(s tn ) being Gaussian as well, it is sufficient to assume that q(s 1n ) is Gaussian, namely at t = 1: q(s 1n ) = N (s 1n ; µ 1n , Γ 1n ).

2) E-Z step: Developing (31) with the same principles as above, one can easily find that the variational posterior distribution of the assignment variable factorizes as:

q(z t ) = D d=1 q(z td ). (38)
In addition, we obtain a closed-form expression for q(z td ):

α tdn = q(Z td = n) = ρ tdn π dn N i=0 ρ tdi π di , (39)
where π dn was defined in [START_REF] Gebru | Audio-visual speaker diarization based on spatiotemporal Bayesian fusion[END_REF], and ρ tdn is given by:

ρ tdn =        N (b td ; Mµ tn , 1 w td Σ) ×e -1 2 tr w td M Σ -1 MΓtn if 1 ≤ n ≤ N U(vol(G))
if n = 0.

(40)

3) M-step: Once the two expectation steps are executed, we maximize J in [START_REF] Ma | Tracking an unknown time-varying number of speakers using tdoa measurements: A random finite set approach[END_REF] with respect to the model parameters, i.e. the covariance matrix of the state dynamics Λ tn . By exploiting again the proposed variational approximation, the dependency of J on Λ tn can be written as:

J(Λ tn ) = E q(stn) log N (s tn ; Dµ t-1,n , DΓ t-1,n D + Λ tn ) ,
which can be further developed as:

J(Λ tn ) = log |DΓ t-1,n D + Λ tn | + Tr (DΓ t-1,n D + Λ tn ) -1 × (41) ((µ tn -Dµ t-1,n )(µ tn -Dµ t-1,n ) + Γ tn ) .
By equating to zero the gradient of (41) w.r.t. Λ tn , we obtain:

Λ tn = Γ tn -DΓ t-1,n D +(µ tn -Dµ t-1,n )(µ tn -Dµ t-1,n ) . (42)

B. Speaker-Birth Process

A birth process is used to initialize new tracks, i.e. speakers that become active. We take inspiration from the birth process for visual tracking proposed in [START_REF] Ba | An online variational bayesian model for multi-person tracking from cluttered scenes[END_REF] and adapt it to audio tracking. The general principle is the following. In a short period of time, say from t -L to t, with L being small (typically 3), we assume that at most one new (yet untracked) speaker becomes active. For each frame from t -L to t, among the observations assigned to n = 0 we select the one with the highest weight, and thus obtain an observation sequence õt-L:t . We then compute the marginal likelihood of this sequence according to our model, 0 = p(õ t-L:t ). If these observations have been generated by a speaker that has not been detected yet (hypothese H 1 ), then they are assumed to be consistent with the model, i.e. exhibit smooth trajectories, and 0 will be high; otherwise, i.e. if they have been generated by background noise (hypothese H 0 ), they will be more randomly spread over the range of possible observations, and 0 will be low. Giving birth to a new speaker track amounts to setting a threshold 1 and deciding between the two hypotheses:

0 H1 > < H0 1 . ( 43
)
This process is applied continuously over time to detect new speakers. This includes speaker track initialization at t = 1.

Note that initially all the assignment variables are set to n = 0 (background noise), namely Z 1d = 0, ∀d.

As for the computaiotn of p(õ t-L:t ), we first rewrite it as the marginalization of the joint probability of the selected observations and the state trajectory ŝt-L:t of a potential speaker:

0 = p(õ t-L:t , ŝt-L:t )dŝ t-L:t , (44)
which, under the proposed model, is given by: All the terms in the above equation have been defined during the derivation of our model except the marginal prior distribution of the state p(ŝ t-L ), and all these terms are Gaussian. For the track-birth process, we just want to test if the trajectory of observations from t-L to t is coherent, and we can define here p(ŝ t-L ) as a non-informative distribution, such as a uniform distribution. In practice we choose a Gaussian distribution with a very large covariance, to ensure a closed-form solution to (45). Due to room limitation, we do not present more details. Let us just mention that in practice we set L = 3, which enables efficient speaker birth detection.

0 = (45) t i=t-L+1 p(õ i |ŝ i )p(ŝ i |ŝ i-1 ) p(õ t-L |ŝ t-L )p(ŝ t-L )dŝ t

C. Speaker Activity Detection

A very interesting feature of the proposed model is that, once speaker tracks have been estimated, the posterior distribution of the assignment variables Z t can be used for speech activity detection, i.e. who are the active speakers at each frame, a task also referred to as speaker diarization in the multi-speaker context. This can be formalized as testing for each frame t and each speaker n between the two following hypotheses: H 1 : Speaker n is active at frame t, and H 0 : Speaker n is silent at frame t. In the present work, this is done by computing the following weighted sum of weights, averaged over a small number of frames L to take into account speaker activity inertia, and comparing with a threshold δ, a test formally written as:

t i=t-L +1 D d=1 α idn w d i H1 > < H0
δ.

(46)

Overall, the variational EM tracking algorithm is described in Algorithm 2.

V. EXPERIMENTS

A. Experimental setups 1) Datasets: We tested and empirically validated our method with the LOCATA and the Kinovis multiple speaker tracking (Kinovis-MST) datasets. The LOCATA (a IEEE-AASP challenge for sound source localization and tracking) [START_REF] Löllmann | The LOCATA challenge data corpus for acoustic source localization and tracking[END_REF] data were recorded in the Computing Laboratory of the Department of Computer Science of Humboldt University Berlin. The room size is 7.1 m × 9.8 m × 3 m, with a reverberation time T 60 ≈ 0.55 s. We report the results of the development corpus for tasks #3 and #5 with a single moving speaker, and for tasks #4 and #6 with two moving speakers, each task comprising three recorded sequences. 2There are twelve microphones arranged such as to form a spherical array and placed on the head of a NAO robot. We used two microphone configurations: four quasi-planar microphones, located on the top of the head, numbered 5, 8, 11, 12, and eight microphones numbered 1, 3, 4, 5, 8, 10, 11, 12. An optical motion capture system was used to provide ground-truth positions of the robot and of the speakers. The participants speak continuously during the entire recordings. However, speech pauses are inevitable and these pauses may last several seconds. Each participant has a headmounted microphone. We applied the voice activity detector [START_REF] Li | Voice activity detection based on statistical likelihood ratio with adaptive thresholding[END_REF] to these microphone signals to obtain ground-truth voice activity information of each participant. The signal-to-noise ratio (SNR) is approximatively 23.4 dB The Kinovis-MST dataset was recorded in the Kinovis multiple-camera laboratory at INRIA Grenoble. 3 The room size is 10.19 m × 9.87 m × 5.6 m, with T 60 ≈ 0.53 s. A v5 NAO robot with four microphones [START_REF] Li | Reverberant sound localization with a robot head based on direct-path relative transfer function[END_REF] was used. The geometric layout of the microphones is similar to the one of the robot used in LOCATA. The speakers were moving around the robot with a speaker-to-robot distance ranging between 1.5 m and 3.5 m. As with LOCATA, a motion capture system was used to obtain ground-truth trajectories of the moving participants and the location of the robot. Ten sequences were recorded with up to three participants, for a total length of about 357 s. The robot's head has built-in fans located nearby the microphones, hence the recordings contain a notable amount of stationary and spatially correlated noise with an SNR of approximatively 2.7 dB [START_REF] Li | Reverberant sound localization with a robot head based on direct-path relative transfer function[END_REF]. The participants behave more naturally than in the LOCATA scenarios, i.e. they take speech turns in a natural multi-party dialog. When one participant is silent, he/she manually hides the infrared marker located on his head to make it invisible to the motion capture system. This provides ground-truth speech activity information for each participant. This dataset and the associated annotations allow us to test the proposed tracking algorithm when the number of active speakers varies over time.

2) Parameter setting: For both datasets, we perform 360 •wide azimuth estimation and tracking: D = 72 azimuth directions at every 5 • in [-175 • , 180 • ] are used as candidate directions. The CGMM mean c i,d f is the head-related transfer function (HRTF) ratio between two microphones, which are precomputed based on the direct-path propagation model for each candidate direction. In the Kinovis-MST dataset, the HRTFs have been measured to compute the CGMM means. For LOCATA, the TDOAs are computed based on the coordinate of microphones, which are then used to compute the phase of the CGMM means, while the magnitude of the CGMM means are set to a constant, e.g. 0.5, for all the frequencies. All the recorded signals are resampled to 16 kHz. The STFT uses the Hamming window with length of 16 ms and shift of 8 ms. The CTF length is Q = 8 frames. The RLS forgetting factor λ is computed using ρ = 1. The smoothing factor β is set to 0.9. The exponentiated gradient update factor is η = 0.07. The smoothing factor η is set to 0.065. The entropy regularization factor is γ = 0.1. For the tracker, the covariance matrix is set to be isotropic Σ = 0.03I 2 . The threshold giving birth to a new identity is τ 1 = 0.75 and L = 3. To decide whether a person is speaking or is silent, L = 3 frames are used, with a threshold δ = 0.15. At each time instance, the VEM algorithm Fig. 3: ROC curve for the LOCATA dataset. has 5 iterations. Corresponding to the STFT frame shift, i.e. 8 ms, the frame rate of the proposed system is 125 frames per second.

3) Comparison with Baseline Methods: The proposed method is evaluated both in "frame-wise localization" mode and in "tracker" mode. In the first mode, the frame-wise online localization module of Section III is applied without the tracker of Section IV. Instead, it is followed by the peak selection process described in [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF]. This method is referred to as DP-RTF using EG (DPRTF-EG). In tracker mode, DPRTF-EG is directly followed by the proposed VEM tracker, without peak selection. It is then simply referred to as VEMtracker. In that case, the directions of active speakers are given by the state variable, and the continuity of the speaker tracks is given by the assignment variable. We compare DPRTF-EG with several baseline methods:

• The standard beamforming-based localization method called SRP using phase transform (PHAT) (SRP-PHAT) [START_REF] Dibiase | Robust localization in reverberant rooms[END_REF]. The same STFT configuration and candidate directions are used for SRP-PHAT and for the proposed method. The steering vector for each candidate direction is derived from the HRTFs and TDOAs for the Kinovis-MST and LOCATA datasets, respectively. The frame-wise SRP is recursively smoothed with a smoothing factor set to 0.065. • A method combining PRP features, CGMM model and parameter update using REM [START_REF] Schwartz | Speaker tracking using recursive EM algorithms[END_REF], referred to as PRP-REM. We also combine the DPRTF features and CGMM with REM (referred to as DPRTF-REM). This is to evaluate the proposed DP-RTF feature w.r.t. PRP, and the EG-based online parameters update method w.r.t. REM. For both baselines, the STFT and CGMM settings are the same as for the proposed method. The updating factor of REM is set to 0.065.

4) Evaluation Metrics:

The detected speakers should be assigned to the actual speakers for performance evaluation. This is done using a greedy matching algorithm. First the azimuth difference for all possible detected-actual speaker pairs are computed, then the detected-actual speaker pair with the smallest difference is picked out as a matched pair. This procedure is iterated until the detected or actual speakers are all picked out. For each matched pair, the detected speaker is then considered to be successfully localized if the azimuth difference is not larger than 15 • . The absolute error is calculated for the successfully localized sources. The mean absolute error (MAE) is computed by averaging the absolute error of all speakers and frames. For the unsuccessful localizations, we count the miss detection (MD) (speaker active but not detected) and false alarms (FAs) (speaker detected but not active). Then the MD and FA rates are computed, using all the frames, as the percentage of the total MDs and FAs out of the total number of actual speakers, respectively. In addition to these localization metrics, we also count the identity switches (IDs) to evaluate the tracking continuity. ID is an absolute number. It represents the number of the identity changes in the tracks for a whole test sequence.

The computation time is measured with the real-time factor (RF), which is the processing time of a method divided by the length of the processed signal. Note that all the methods are implemented in MATLAB.

B. Results for LOCATA Dataset

For convenience, both the spatial spectrum of SRP-PHAT and the CGMM component weights profile will be referred to as heatmaps. Fig. 2 shows an example of a result obtained with a LOCATA sequence. Two speakers are moving and continuously speaking with short pauses. The SRP-PHAT heatmap (Fig. 2 (b)) is cluttered due to the non ideal beampattern of the microphone array and to the influence of reverberation and noise. For most of the time, SRP-PHAT has prominent response power for the true speaker directions. Localization of the most dominant speaker can be made by selecting the direction with the largest response power. However, it is difficult to correctly count the number of active speakers and localize less dominant speakers, since there exist a number of spurious peaks. PRP-REM (Fig. 2 (c)) exhibits a clearer heatmap compared to SRP-PHAT, but there exist some spurious trajectories as well, since the PRP features are contaminated by reverberation. DPRTF-REM (Fig. 2 (d)) removes most of the spurious trajectories, which illustrates the robustness of the proposed DP-RTF feature against reverberation. From Fig. 2 (e), it can be seen that the proposed EG algorithm further removes the interferences by applying the entropy regularization. In addition, the peak evolution is smoother compared with Fig. 2 (d), which is mainly due to the use of the spatial smoothing. Fig. 2 (f) illustrates the result obtained with the proposed VEM tracker, with DPRTF-EG providing the observations. The proposed tracker gives smoother and cleaner results compared with the other methods. Even when the observations have a low weight, the tracker is still able to give the correct speaker trajectories. This is ensured by the second term in [START_REF] Li | Voice activity detection based on statistical likelihood ratio with adaptive thresholding[END_REF] which exploits the source dynamics model and continues to provide localization information even when w t,d (and/or α tdn ) becomes small. As a result, the tracker is able to preserve the identity of speakers in spite of the (short) speech pauses. In the presented sequence example, the estimated speaker identities are quite consistent with the ground truth.

To empirically evaluate the quality of the heatmaps provided by the localization methods, we computed the receiver operating characteristic (ROC) curve (MD rate versus FA rate) for the LOCATA dataset by varying the peak selection threshold, for each tested method, Fig. 3. For the ROC curve, the closer to the left-bottom the better. As already mentioned, in addition to using four microphones, we also tested an eight-microphone configuration, which is referred to as DPRTF-EG-8ch. By analyzing the ROC curves, one notices that the methods based on DP-RTF perform better than SRP-PHAT and than PRP-REM, which is consistent with the heatmaps of Fig. 2: SRP-PHAT and PRP-REM are more sensitive to the presence of reverberations than the proposed methods. The performance of both DPRTF-REM and DPRTF-EG cannot be easily discriminated using the ROC curves. DPRTF-EG-8ch performs slightly better than DPRTF-EG, which means that the performance of the proposed method can be slightly improved by increasing the number of microphones. One may conclude that the proposed method is well suited when only a small number of microphones are available. With all methods, the FA rate can be trivially decreased to be close to 0 by increasing the peak selection threshold. However, the MD rate cannot be decreased to 0 even with a very small peak-selection threshold, since some speech frames that are actually present cannot be detected as the heatmap peaks due to the influence of noise For each curve, a good balance between FA rate and MD rate is achieved at the left-bottom corner, which can be detected as the point with the minimum distance to the origin. The average localization results corresponding to this optimal left-bottom point are summarized in Table I for each tested method. It can be seen that, besides MD and FA, the DPRTF-based methods achieve smaller MAE than SRP-PHAT and PRP-REM, since the proposed DP-RTF features are robust against reverberation and thus leads to smaller biases for the heatmap peaks. DPRTF-EG has a higher MD rate than DPRTF-REM, while it also has lower FA rate, and a lower MAE, due to the effect of entropy regularization. With eight microphones, i.e. DPRTF-EG-8ch, MD is 1% smaller than the MD of DPRTF-EG, since the use of a non coplanar microphone setup provides more accurate localization than a coplanar setup. The proposed tracker performs the best in terms of MD and of FA. For the four-microphone configuration, the tracker slightly reduces FA compared to DPRTF-EG. It also reduces the MD score since some correct speaker trajectories can be recovered even when the observations have (very) low weights, as explained above. In addition, the MAE is noticeably reduced when more microphones are used by the VEM tracker, which is not the case with the DPRTF-EG localizer. This phenomenon indicates that, compared with the localizer, the tracker is able to better exploit additional information available with extra microphones, namely to revise the speaker trajectory estimation, since the state dynamics of the tracker helps correcting the possibly inaccurate additional localization information. The proposed tracker achieves quite consistent speaker ID estimation. For the whole LOCATA dataset, only ten identity switches were observed when using DPRTF-EG, and this number is reduced to six when using DPRTF-EG-8ch. The remaining identity switches are mainly due to speakers with crossing trajectories, a hard case for multiple audio-source tracking.

As for the computation time, SRP-PHAT has the smallest RF. Based on the fact that the RFs of DPRTF-REM and DPRTF-EG are identical, we can conclude that the REM algorithm and the proposed EG algorithm have comparable computational complexities. The RFs of PRP-REM, DPRTF-REM (or DPRTF-EG) and DPRTF-EG-8ch are different due to different computational complexities for feature estimation, more precisely due to the different dimensions of the vector to be estimated. The CTF identification used for DP-RTF estimation solves an RLS problem with the unknown CTF vector ãf ∈ C (IQ-1)×1 . Remind that I and Q denote the number of microphones and the CTF length, respectively. In the present work, we have set I = 4/Q = 8 for DPRTF-REM (or DPRTF-EG), I = 8/Q = 8 for DPRTF-EG-8ch. PRP is defined based on the narrow-band assumption, or equivalently based on the CTF with Q = 1, thence we have I = 4/Q = 1 for PRP-REM. The proposed localization method, i.e. DPRTF-EG with four microphones, has an RF smaller than one, which means it can be run in real time. The RF for the proposed tracker (VEM) is computed by the sum of the localization time and of the tracking time. For acoustic tracking, the tracker observes an direction of arrival (DOA) estimate every 8 ms. However, an 8 ms speaker motion is small. Thus in practice, the tracker uses one DOA estimate per 32 ms intervals, which leads to an RF of 2.05 for the four-channel (4ch) case and 4.05 for the eight-channel (8ch) case. The RF can be further improved by using less DOA estimates.

C. Results for Kinovis-MST Dataset

Fig. 4 shows an example of result for a Kinovis-MST sequence. Three participants are moving and intermittently speaking. It can be seen that, for many frames, the response power of SRP-PHAT and the CGMM component weights of PRP-REM corresponding to the true active speakers are not prominent, compared to the spurious trajectories. Again, DPRTF-REM and DPRTF-EG provide much better heatmaps, though they also miss some speaking frames, e.g. at the beginning of Speaker 3's trajectory (in blue). The possible reasons are i) the NAO robot (v5) has a relative strong egonoise [START_REF] Li | Reverberant sound localization with a robot head based on direct-path relative transfer function[END_REF], and thus the signal-to-noise ratio of the recorded signals is relative low, and ii) the speakers are moving with a varying source-to-robot distance and the direct-path speech is contaminated by more reverberations when the speakers are distant. Overall, DPRTF-REM and DPRTF-EG are able to monitor the moving, appearance, and disappearance of active speakers for most of the time, with a small time lag due to the temporal smoothing. This kind of recording/scenario is very challenging for the tracking method, especially for speaker identity preservation, since the participants are intermittently speaking and moving. In a general manner, the proposed tracker achieves relatively good results, as illustrated in Fig. 4 (f). The tracked trajectories are smooth and clean. If the true trajectory of one speaker has an approximately constant direction, the tracker is able to reidentify the speaker even after a long silence thanks to the above-mentioned combination of observations and dynamics in [START_REF] Li | Voice activity detection based on statistical likelihood ratio with adaptive thresholding[END_REF], e.g. Speaker 1's trajectory in red. In the case that the speaker changes his/her movement when he/she is silent, the track can be lost. When the person speaks again, it is indeed difficult to re-idendify him/her based on the dynamics estimated before the silence period. The tracker may then prefer to give birth to a new speaker. This is illustrated by the black trajectory turning into green, and the blue trajectory turning into cyan in Fig. 4. Note that the silence periods are here much longer than in the LOCATA example of Fig. 2. II summarizes the localization and tracking results for the optimal bottom-left point of the ROC curves. It can be seen that, for the four localization methods, MAEs are quite close, namely the heatmap peaks have similar biases. Compared with the results for the LOCATA dataset, the advantage of the proposed tracker is more significant for the Kinovis-MST dataset. In particular, the FA rate is reduced by 7.8% relatively to DPRTF-EG, and is similar to the FA rate obtained with the LOCATA dataset. This means that the dynamic model associated with the tracker can efficiently reduce the influence of incorrect source localizations caused by noise and by complex source movements. The identity switches are mainly caused by speakers changing their direction of movement while during silent periods, as discussed above. Compared to the LOCATA dataset, DPRTF-EG has smaller RF, since the Kinovis-MST dataset is noisier and more noise frames are skipped in the RLS algorithm.

D. Dicussion

The experimental results obtained with the two datasets clearly show the effectiveness of the proposed method based on DP-RTF estimation, multiple speaker localization and variational tracking. To improve robustness, temporal smoothing is used, which leads to localization/tracking latency. This latency causes MD and FA observed at both the beginning and the end of continuous speech segments. However, it can be observed from the examples shown in Fig. 2 and Fig. 4 that the latency is not that severe. The Kinovis-MST dataset is more challenging than the LOCATA dataset for speaker localization/tracking, due to its lower SNR and the presence of casual speaking style. Even though, the proposed methods achieve a comparable FA rate with the two datasets. Concerning the MD score, when applied to Kinovis-MST, the method yields larger MD rates than when applied to LOCATA. This is due to the large number of TF bins dominated by a high SNR score, present in the Kinovis-MST recordings. The tracker's dynamics attenuate the influence of these TF bins to a limited extend.

VI. CONCLUSION

In this paper, we proposed and combined i) a recursive DP-RTF feature estimation method, ii) an online multiplespeaker localization method, and iii) an multiple-speaker tracking method. The resulting framework provides online speaker counting, localization and consistent tracking (i.e. preserving speaker identity over a track in spite of intermittent speech production). The three algorithms are computationally efficient. In particular the tracking algorithm implemented in variational Bayesian framework yields a tractable solver under the form of VEM. Experiments with two datasets, recorded in realistic environment, verify that the proposed method is robust against reverberation and noise. Moreover, the tracker is able to efficiently track multiple moving speakers, detect whether they are speech or they are silent, as long as the motion associated with silent people is smooth. However, the tracking of the person from silent to active remains a difficult task. The combination of the proposed method with speaker identification will be addressed in the future.

The proposed VEM tracker can be easily adapted to work in tandem with any frame-wise localizer providing source location estimates and/or corresponding weights (and if no weights are provided by the localizer, the tracker can be applied with all weights set to one). This makes the proposed tracker very flexible, and easily reusable by the audio processing research community.

  compute the exp. gradient: {r t-1,d } D d=1 Eq. (10) updates the GMM weights: {w td } D d=1 Azimuth: {b td = (cos( θd ), sin( θd )) } D d=1 Observed variables: {o td = [b td ; w td ]} D d=1 EM iteration: E-Z: Eq. (39), assignments {α tdn } D,N d=1,n=1 E-S: Eq. (36), posterior covariances {Γ tn } N n=1 Eq. (37), posterior means {µ tn } N n=1 M-step: Eq. (42), covariances of state dynamics {Λ tn } N n=1 Output: Azimuth directions {s tn = µ tn } N n=1 and angular velocities of tracked speakers Candidate directions: { θd } D d=1

Fig. 2 :

 2 Fig. 2: Results of speaker localization and tracking for Recording 1 / Task 6 of LOCATA data. (a) Ground truth trajectory and voice activity (red for speaker 1, black for speaker 2). Intervals in the trajectories are speaking pauses. (b)-(e) One-dimensional heat maps as a function of time for the four tested localization methods. (f) Results for the proposed VEM-based tracker. Black and red colors demonstrate a successful tracking, i.e. continuity of the tracks despite of speech pauses.

Fig. 4 :

 4 Fig. 4: Results of speaker localization and tracking for one sequence of the Kinovis-MST dataset. (a) Ground truth trajectory and voice activity (red for speaker 1, black for speaker 2, blue for speaker 3). (b)-(e) One-dimensional heat maps as a function of time for the four tested localization methods. (f) Results for the proposed VEM-based tracker.

Fig. 5 :

 5 Fig. 5: ROC curve for the Kinovis-MST dataset.

Fig 5

 5 Fig 5 show the ROC curves for the Kinovis-MST dataset. Compared to the ROC curves for the LOCATA dataset, all the four localization methods have a worse ROC curve, especially along the MD rate axis, for the reasons mentioned above. TableIIsummarizes the localization and tracking results for the optimal bottom-left point of the ROC curves. It can be seen that, for the four localization methods, MAEs are quite close, namely the heatmap peaks have similar biases. Compared with the results for the LOCATA dataset, the advantage of the proposed tracker is more significant for the Kinovis-MST dataset. In particular, the FA rate is reduced by 7.8% relatively to DPRTF-EG, and is similar to the FA rate obtained with the LOCATA dataset. This means that the dynamic model associated with the tracker can efficiently reduce the influence of incorrect source localizations caused by noise and by complex source movements. The identity switches are mainly caused by speakers changing their direction of movement while during silent periods, as discussed above. Compared to the LOCATA dataset, DPRTF-EG has smaller RF, since the Kinovis-MST dataset is noisier and more noise frames are skipped in the RLS algorithm.

TABLE I :

 I Localization and tracking results for the LOCATA data.

		MD rate (%) FA rate (%) MAE ( • ) IDs RF
	SRP-PHAT	39.2	18.6	5.2	-	0.06
	PRP-REM	30.9	19.6	5.0	-	0.30
	DPRTF-REM	23.3	15.2	4.6	-	0.97
	DPRTF-EG	23.9	13.0	4.0	-	0.97
	DPRTF-EG-8ch	22.7	13.2	4.1	-	3.03
	VEM + EG	22.7	12.4	4.1	10 2.05
	VEM + EG-8ch	22.9	11.0	3.2	6 4.05

TABLE II :

 II Localization and tracking results for the Kinovis-MST dataset.

		MD rate (%) FA rate (%) MAE ( • ) IDs RF
	SRP-PHAT	60.0	37.1	5.5	-	0.07
	PRP-REM	40.3	23.1	5.1	-	0.32
	DPRTF-REM	37.6	22.0	5.5	-	0.73
	DPRTF-EG	31.4	19.5	5.3	-	0.73
	VEM + EG	31.1	11.7	4.9	11 2.12

The results obtained with the proposed method were officially submitted to the LOCATA challenge and they will be available soon at https://locata. lms.tf.fau.de/.

The Kinovis-MST dataset is publicly available at: https://team.inria.fr/ perception/the-kinovis-mst-dataset/

This work was supported by the ERC Advanced Grant VHIA #340113.