Tracking Multiple Audio Sources with the von Mises Distribution and Variational EM

TL;DR

This paper introduces a novel variational EM approach for tracking multiple moving audio sources using the von Mises distribution, enabling efficient and accurate source trajectory estimation in dynamic environments.

Contribution

It proposes a Bayesian filtering framework with a variational approximation for circular data, and introduces a new source birth method for improved initialization and detection.

Findings

Effective tracking of moving sources demonstrated on LOCATA dataset.

The variational EM algorithm is computationally efficient and accurate.

The source birth method improves detection and trajectory estimation.

Abstract

In this paper we address the problem of simultaneously tracking several moving audio sources, namely the problem of estimating source trajectories from a sequence of observed features. We propose to use the von Mises distribution to model audio-source directions of arrival with circular random variables. This leads to a Bayesian filtering formulation which is intractable because of the combinatorial explosion of associating observed variables with latent variables, over time. We propose a variational approximation of the filtering distribution. We infer a variational expectation-maximization algorithm that is both computationally tractable and time efficient. We propose an audio-source birth method that favors smooth source trajectories and which is used both to initialize the number of active sources and to detect new sources. We perform experiments with the recently released LOCATA dataset comprising two moving sources and a moving microphone array mounted onto a robot.

Tables

Table 1. TABLE I: Method evaluation with the LOCATA dataset.

Method	MD (%)	FA (%)	MAE (°)
vM-PHD [13]	33.4	9.5	4.5
GM-ZO [16]	27.0	10.8	4.7
GM-FO [16]	22.3	6.3	3.2
vM-VEM (proposed)	23.9	5.9	2.6

Equations

$$\displaystyle p(\mbox{\boldmath$s$}_{t},\mbox{\boldmath$z$}_{t}|\mbox{\boldmath$y$}_{1:t})\propto p(\mbox{\boldmath$y$}_{t}|\mbox{\boldmath$z$}_{t},\mbox{\boldmath$s$}_{t})p(\mbox{\boldmath$z$}_{t})p(\mbox{\boldmath$s$}_{t}|\mbox{\boldmath$y$}_{1:t-1}),$$

$$p(\mbox{\boldmath$y$}_{t}|\mbox{\boldmath$z$}_{t},\mbox{\boldmath$s$}_{t})=\prod_{m=1}^{M_{t}}p(y_{tm}|\mbox{\boldmath$z$}_{t},\mbox{\boldmath$s$}_{t}).$$

$$p(y_{tm}|Z_{tm}=n,s_{tn})=\{\begin{array}{ll}\mathcal{M}(y_{tm};s_{tn},{\kappa }_y{\omega }_{tm})& n\ne 0\\ \mathcal{U}(y_{tm})& n=0\end{array},$$

$$p(\mbox{\boldmath$z$}_{t})=\prod_{m=1}^{M_{t}}p(Z_{tm}=n),$$

$$\displaystyle p(\mbox{\boldmath$s$}_{t}|\mbox{\boldmath$y$}_{1:t-1})$$

$$p(\mbox{\boldmath$s$}_{t}|\mbox{\boldmath$s$}_{t-1})=\prod_{n=1}^{N}\mathcal{M}(s_{tn};s_{t-1,n},\kappa_{d}),$$

$$p(\mbox{\boldmath$s$}_{t},\mbox{\boldmath$z$}_{t}|\mbox{\boldmath$y$}_{1:t})\approx q(\mbox{\boldmath$s$}_{t})q(\mbox{\boldmath$z$}_{t}).$$

$$\displaystyle q(\mbox{\boldmath$z$}_{t})$$

$$\displaystyle q(\mbox{\boldmath$s$}_{t})$$

$$Q(\Theta,\tilde{\Theta})=\mathrm{E}_{q(\mbox{\boldmath$s$}_{t})q(\mbox{\boldmath$z$}_{t})}\Big{[}\log p(\mbox{\boldmath$y$}_{t},\mbox{\boldmath$s$}_{t},\mbox{\boldmath$z$}_{t}|\mbox{\boldmath$y$}_{1:t-1},\Theta,\tilde{\Theta})\Big{]}.$$

$$q(\mbox{\boldmath$z$}_{t})=\prod_{m=1}^{M_{t}}q(z_{tm}),\quad q(\mbox{\boldmath$s$}_{t})=\prod_{n=1}^{N}q(s_{tn}),$$

$$p(s_{tn}|\mbox{\boldmath$y$}_{1:t-1})=\int p(s_{tn}|s_{t-1,n})p(s_{t-1,n}|\mbox{\boldmath$y$}_{1:t-1})\textrm{d}s_{t-1,n}.$$

$$p(s_{tn}|\mbox{\boldmath$y$}_{1:t-1})\approx\mathcal{M}(s_{tn};\mu_{t-1,n},\tilde{\kappa}_{t-1,n}),$$

$${\tilde{\kappa }}_{t-1,n}=A^{-1}(A({\kappa }_{t-1,n})A({\kappa }_d)),$$

$${\mu }_{tn}$$

$$\left(\frac{{\kappa }_y{\sum }_{m=1}^{M_t}{\alpha }_{tmn}{\omega }_{tm}\sin (y_{tm})+{\tilde{\kappa }}_{t-1,n}\sin ({\mu }_{t-1,n})}{{\kappa }_y{\sum }_{m=1}^{M_t}{\alpha }_{tmn}{\omega }_{tm}\cos (y_{tm})+{\tilde{\kappa }}_{t-1,n}\cos ({\mu }_{t-1,n})}\right),$$

$${\kappa }_{tn}$$

$$+2({\kappa }_y{)}^2\sum _{m=1}^{M_t}\sum _{l=m+1}^{M_t}{\alpha }_{tmn}{\omega }_{tm}{\alpha }_{tln}{w}_{tl}\cos (y_{tm}-y_{tl})$$

$${+2{\kappa }_y{\tilde{\kappa }}_{t-1,n}\sum _{m=1}^{M_t}({\alpha }_{tmn}{\omega }_{tm}\cos (y_{tm}-{\mu }_{t-1,n})))}^{1/2},$$

$${\alpha }_{tmn}=q(z_{tm}=n)=\frac{{\pi }_n{\beta }_{tmn}}{{\sum }_{l=0}^N{\pi }_l{\beta }_{tml}}$$

$${\beta }_{tmn}=\{\begin{array}{ll}{\omega }_{tm}{\kappa }_{y}A({\omega }_{tm}{\kappa }_y)\cos (y_{tm}-{\mu }_{tn})& n\ne 0\\ 1/(2\pi )& n=0,\end{array}$$

$${\tau }_j=p({\hat{y}}_{t-L:t}^j)=\int p({\hat{y}}_{t-L:t}^j,s_{t-L:t})\text{d}{s}_{t-L:t}.$$

$${\tau }_j=\prod _{l=0}^L\frac{I_0({\overline{\kappa }}_{t-l}^j)}{2\pi I_0({\kappa }_y{\hat{\omega }}_{t-l}^j)I_0({\hat{\kappa }}_{t-l}^j)},$$

$${\overline{\kappa }}_{t-l}^j=\sqrt{({\hat{\kappa }}_{t-l}^j{)}^2+({\kappa }_y{\hat{\omega }}_{t-l}^j{)}^2+{\hat{\kappa }}_{t-l}^j{\kappa }_y{\hat{\omega }}_{t-l}^j\cos ({\hat{y}}_{t-l}^j-{\hat{\mu }}_{t-l}^j)},$$

$${\hat{\mu }}_{t-l+1}^j={\tan }^{-1}\left(\frac{{\hat{\kappa }}_{t-l}^j\sin ({\hat{\mu }}_{t-l}^j)+{\kappa }_y{\hat{\omega }}_{t-l}^j\sin ({\hat{y}}_{t-l}^j)}{{\hat{\kappa }}_{t-l}^j\cos ({\hat{\mu }}_{t-l}^j)+{\kappa }_y{\hat{\omega }}_{t-l}^j\cos ({\hat{y}}_{t-l}^j)}\right),$$

$${\hat{\kappa }}_{t-l+1}^j=A^{-1}(A({\tilde{\kappa }}_{t-l}^j)A({\kappa }_d)).$$

$$\sum _{t^{\prime }=t-D}^t\sum _{m=1}^{M_t}{\alpha }_{t^{\prime }mn}{\omega }_{t^{\prime }m}\genfrac{}{}{0ex}{}{\genfrac{}{}{0ex}{}{active}{>}}{\genfrac{}{}{0ex}{}{<}{silent}}\delta$$

$$q(\mbox{\boldmath$s$}_{t})\propto\exp\Big{(}\mathrm{E}_{q(\mbox{\boldmath$z$}_{t})}\log p(\mbox{\boldmath$s$}_{t},\mbox{\boldmath$z$}_{t}|\mbox{\boldmath$y$}_{1:t})\Big{)}.$$

$$q(\mbox{\boldmath$s$}_{t})\propto\exp\Big{(}\mathrm{E}_{q(\mbox{\boldmath$z$}_{t})}\log p(\mbox{\boldmath$y$}_{t}|\mbox{\boldmath$s$}_{t},\mbox{\boldmath$z$}_{t})\Big{)}p(\mbox{\boldmath$s$}_{t}|\mbox{\boldmath$y$}_{1:t-1}).$$

$$\displaystyle\mathrm{E}_{q(\mbox{\boldmath$z$}_{t})}\log p(\mbox{\boldmath$y$}_{t}|\mbox{\boldmath$s$}_{t},\mbox{\boldmath$z$}_{t})$$

$$\displaystyle=\mathrm{E}_{q(\mbox{\boldmath$z$}_{t})}\sum_{m=1}^{M_{t}}\log p(y_{tm}|\mbox{\boldmath$s$}_{t},z_{tm})$$

$$\displaystyle=\sum_{m=1}^{M_{t}}\mathrm{E}_{q(z_{tm})}\log p(y_{tm}|\mbox{\boldmath$s$}_{t},z_{tm})$$

$$\displaystyle=\sum_{m=1}^{M_{t}}\sum_{n=0}^{N}q(z_{tm}=n)\log p(y_{tm}|\mbox{\boldmath$s$}_{t},z_{tm}=n)$$

$$=\sum _{m=1}^{M_t}\sum _{n=0}^N{\alpha }_{tnm}\log p(y_{tm}|s_{tn},z_{tm}=n)$$

$$=\sum _{m=1}^{M_t}\sum _{n=0}^N{\alpha }_{tnm}\log \mathcal{M}(y_{tm};s_{tn},{\omega }_{tm}{\kappa }_y)$$

$$\displaystyle\stackrel{{\scriptstyle\mbox{\boldmath$s$}_{t}}}{{=}}\sum_{m=1}^{M_{t}}\sum_{n=0}^{N}\alpha_{tnm}\omega_{tm}\kappa_{y}\cos(y_{tm}-s_{tn}),$$

$$\displaystyle q(\mbox{\boldmath$s$}_{t})\propto$$

Full text

I Introduction

We address the problem of tracking several moving audio sources. Audio tracking is useful for audio-source separation, spatial filtering, speaker diarization, speech enhancement and speech recognition, which in turn are essential methodologies, e.g. home assistants. Audio-source tracking is difficult because audio signals are adversely affected by noise, reverberation and interferences between acoustic signals.

Single-source tracking methods are based on observing time differences of arrival between microphones. Since the mapping between TDOAs and the source locations is non-linear, sequential Monte Carlo approaches are used, e.g. [1, 2, 3]. Alternatively, directions of arrival can be used. The problem is cast into a linear dynamic model, e.g. [4]. In this case source directions should however be modeled as circular random variables, e.g. the wrapped Gaussian distribution [5], or the von Mises distribution [6, 7].

Multiple-source tracking is more challenging: (i) the number of active sources is unknown and varies over time, (ii) several DOAs need be detected, and (iii) DOA-to-source assignments must be estimated. An unknown number of sources was addressed using random finite sets [8]. Since the probability density function (pdf) is computationally intractable, its first-order approximation can be propagated in time using the probability hypothesis density (PHD) filter [8, 9]. In [10] the PHD filter was applied to audio recordings to track multiple sources from TDOA estimates. In [11] the wrapped Gaussian distribution is incorporated within a PHD filter.

The von Mises-Fisher distribution was used in [12] to build a factorial filter.

A mixture of von Mises distributions was combined with a PHD filter in [13]. The main drawback of PHD filters is that explicit observation-to-source associations are not established. Instead, post-processing techniques are required for track labelling [14].

A variational approximation of the multiple target tracking was addressed in [15]: bservation-to-target associations are discrete latent variables which are estimated with an variational expectation maximization (VEM) solver. Moreover, the problem of tracking a varying number of targets is addressed via track-birth and track-death processes. The variational approximation of [15] was recently extended to track multiple speakers with audio [16] and audio-visual data [17].

This paper builds on [7, 15, 16] and proposes to use the von Mises distribution to model the DOAs of multiple acoustic sources with circular random variables.

The Bayesian filtering formulation for the multi-source tracking problem is intractable over time, due to the combinatorial nature of the unknown association between observed variables and latent variables. We propose a variational approximation of the filtering distribution. A novel mathematical framework is therefore proposed in order to deal with a mixture of von Mises distributions. The contribution of this paper is therefore a novel VEM algorithm that is both computationally tractable and time efficient.

Moreover, we propose an audio-source birth method that favors smooth source trajectories and which is used both to initialize the number of active sources and to detect new sources. We perform experiments with the recently released LOCATA dataset [18] comprising audio recordings of two moving sources from a moving microphone array in a real acoustic environment.

The paper is organized as follows. Section II describes the probabilistic model and Section III describes a variational approximation of the filtering distribution and the VEM algorithm. Section IV briefly describes the source birth method. Experiments and comparisons with other methods are described in Section V. Supplemental materials (mathematical derivations, software and videos) are available online.111https://team.inria.fr/perception/research/audiotrack-vonm/

II The Filtering Distribution

Let $N$ be the number of audio sources. Let \mbox{\boldmathy}_{t}=\{y_{t1},\dots y_{tm},\dots y_{tM_{t}}\} be the set of $M_{t}$ observed DOAs at time step $t$. Let \mbox{\boldmaths}_{t}=\{s_{t1},\dots,s_{tn},\dots s_{tN}\} be the set of $N$ latent DOAs, where $s_{tn}$ is the DOA of source $n$ and time $t$. Observed and source DOAs are realizations of random circular variables $Y$ and $S$, respectively, in the interval $]-\pi,\pi]$, i.e. azimuth directions. Let $Z_{tm}$ be a discrete association variable whose realizations take values in $\{0,1,\dots N\}$, i.e. $Z_{tm}=n$ means that observation $y_{tm}$ is assigned to source $n$ and $Z_{tm}=0$ means that the observation is “clutter”, hence assigned to none of the $N$ sources – we refer to “0” as a dummy source. For convenience, we also use the notation \mbox{\boldmathz}_{t}=\{z_{t1},\dots z_{tm},\dots z_{tM_{t}}\}.

Within a Bayesian model, multiple target tracking can be formulated as the estimation of the filtering distribution p(\mbox{\boldmaths}_{t},\mbox{\boldmathz}_{t}|\mbox{\boldmathy}_{1:t}), with the notation \mbox{\boldmathy}_{1:t}=(\mbox{\boldmathy}_{1},\dots\mbox{\boldmathy}_{t}). We assume that variables $s_{tn}$ follow a first-order Markov model, and that observations only depend on the current state and on the assignment variables. Moreover, we assume that the assignment variable does not depend on the previous observations. Under these assumptions the posterior, or filtering, pdf is given by:

[TABLE]

where p(\mbox{\boldmathy}_{t}|\mbox{\boldmathz}_{t},\mbox{\boldmaths}_{t}) is the observation likelihood, p(\mbox{\boldmathz}_{t}) is the prior pdf of the assignment variables and p(\mbox{\boldmaths}_{t}|\mbox{\boldmathy}_{1:t-1}) is the predictive pdf of the latent variables.

II-1 Observation likelihood

Assuming that observed DOAs are independently and identically distributed (i.i.d.), the observation likelihood can be written as:

[TABLE]

The likelihood that a DOA corresponds to a source is modeled by a von Mises distribution [7], whereas the likelihood that a DOA corresponds to a dummy source (e.g. noise) is modeled by a uniform distribution:

[TABLE]

where $\mathcal{M}(y\,;s,\kappa)=(2\pi I_{0}(\kappa))^{-1}\exp\{\kappa\cos(y-s)\}$ denotes the von Mises distribution with mean $s$ and concentration $\kappa$, $I_{p}(\cdot)$ denotes the modified Bessel function of the first kind of order $p$, $\kappa_{y}$ denotes the concentration of audio observations, $\omega_{tm}\in[0,1]$ is a confidence associated with each observation, and $\mathcal{U}(y_{tm})=(2\pi)^{-1}$ denotes the uniform distribution along the support of the unit circle.

II-2 Prior pdf of the assignment variables

Assuming that assignment variables are i.i.d., the joint prior pdf is given by:

[TABLE]

and we denote with $\pi_{n}=p(Z_{tm}=n)$, $\sum_{n=0}^{N}\pi_{n}=1$, the prior probability that source $n$ is associated with $y_{tm}$.

II-3 Predictive pdf of the latent variables

The predictive pdf extrapolates information inferred in the past to the current time step using a dynamic model for the source motion, i.e. DOA rotation:

[TABLE]

where $p(s_{t}|s_{t-1})$ denotes the prior pdf of the source motion and $p(s_{t-1}|y_{1:t-1})$ is the filtering pdf at $t-1$. The sources are assumed to move independently, and each source (DOA) follows a von Mises distribution:

[TABLE]

where $\kappa_{d}$ is the concentration of the state dynamics. $\Theta=\{\kappa_{y},\kappa_{d},\pi_{0},\ldots,\pi_{N}\}$ denotes the set of model parameters.

As already mentioned in Section I, the filtering distribution corresponds to a mixture model whose number of components grows exponentially along time, therefore solving (1) directly is computationally intractable. Below we infer a variational approximation of (1) which drastically reduces the explosion of the number of mixture components; consequently, it leads to a computationally tractable algorithm.

III Variational Approximation and Algorithm

Since solving (1) is computationally intractable, we propose to approximate the conditional independence between the latent and the assignment variables given all observations up to the current time step, t,, more precisely

[TABLE]

The proposed factorization leads to a VEM algorithm [19], where the posterior distribution of the two variables are found by two variational E-steps:

[TABLE]

where ($E[\cdot]$ is the expectation operator). The model parameters $\Theta$ are estimated by maximizing the expected complete-data log-likelihood:

[TABLE]

where $\tilde{\Theta}$ are the old parameters.

By combining the i.i.d. assumption, i.e. (2), with the variational factorization (7), we observe that the posterior pdf of the assignment variables and the posterior pdf of the latent variables can be factorized:

[TABLE]

and, therefore, the predictive pdf is separable:

[TABLE]

Moreover, assuming that the filtering pdf at $t-1$ follows a von Mises distribution, i.e. $q(s_{t-1,n})=\mathcal{M}(s_{t-1,n};\mu_{t-1,n},\kappa_{t-1,n})$, then the predictive pdf is approximately a von Mises distribution (see [7], [20, (3.5.43)]):

[TABLE]

where the predicted concentration parameter, $\tilde{\kappa}_{t-1,n}$, is:

[TABLE]

and where $A(a)=I_{1}(a)/I_{0}(a)$, and $A^{-1}(a)\approx(2a-a^{3})/(1-a^{2}).$ Using (8), (9) and (10), the filtering distribution is therefore obtained by iterating through three steps, i.e. the E-S, E-Z and M steps, provided below (detailed mathematical derivations can be found in the appendices).

III-1 E-S step

Inserting (1) and (12) in (9), $q(s_{tn})$ reduces to a von Mises distribution, $\mathcal{M}(s_{tn};\mu_{tn},\kappa_{tn})$. The mean $\mu_{tn}$ and concentration $\kappa_{tn}$ are given by:

[TABLE]

where $\alpha_{tmn}=q(Z_{tm}=n)$ denotes the variational posterior probability of the assignment variables. Therefore, the expressibility of the posterior distribution as a mixture of von Mises propagates over time, and only needs to be assumed at $t=1$. Please consult the supplementary materials for more details.

III-2 E-Z step

By computing the expectation over \mbox{\boldmaths}_{t} in (8), the following expression is obtained:

[TABLE]

where $\beta_{tmn}$ is given by (please consult the supplementary materials for a detailed derivation):

[TABLE]

III-3 M step

The parameter set $\Theta$ is evaluated by maximizing (10). The priors (4) are obtained using the conventional update rule [19]: $\pi_{n}\propto\sum_{m=1}^{M_{t}}\alpha_{tnm}$. The concentration parameters, $\kappa_{y}$ and $\kappa_{d}$, are evaluated using gradient descent (please consult the supplementary materials). Based on the E-S-step, E-Z-steo and M-step formulas above, the proposed VEM algorithm iterates until convergence at each time step, in order to estimate the posterior distributions and to update the estimated model parameters.

IV Audio-Source Birth Process

We now describe in detail the proposed birth process which is essential to initialize the number of audio sources as well as to detect new sources at any time. The birth process gathers all the DOAs that were not assigned to a source, i.e. assigned to $n=0$, at current time $t$ as well over the $L$ previous times ($L=2$ in all our experiments). From this set of DOAs we build DOA/observation sequences (one observation at each time $t$) and let $\hat{y}_{t-L:t}^{j}$ be such a sequence of DOAs, where $j$ is the sequence index. We consider the marginal likelihood:

[TABLE]

Using (12) and the harmonic sum theorem, the integral (17) becomes (please consult the supplementary materials):

[TABLE]

where $\hat{\omega}_{t}$ is the confidence associated with $\hat{y}_{t}$. The concentration parameters, $\overline{\kappa}_{t-l}^{j}$ and $\hat{\kappa}_{t-l+1}^{j}$, depend on the observations and are recursively computed for each sequence $j$:

[TABLE]

The sequence $j^{\ast}$ with the maximal marginal likelihood (18), namely $j^{\ast}=\operatornamewithlimits{argmax}_{j}(\tau_{j})$, is supposed to be generated from a not yet known audio source only if $\tau_{j^{\ast}}$ is larger than a threshold $\tau_{0}$: a new source $\tilde{n}$ is created in this case and $q(s_{t\tilde{n}})=\mathcal{M}(s_{t\tilde{n}};\hat{\mu}_{tj^{\ast}},\hat{\kappa}_{tj^{\ast}})$.

We note that, in practice, a source may become silent. In this case, the source is no longer associated with observations, and the proposed tracking algorithm relies solely on the source dynamics. If a source is silent for a long time the algorithm loses track of that source. If, after a while, the source becomes active again, a new track is initialized.

V Experimental Evaluation

The proposed method was evaluated using the audio recordings from Task 6 of the IEEE-AASP LOCATA222https://locata.lms.tf.fau.de/ challenge development dataset [18], which involves multiple moving sound sources, i.e. speakers, and a microphone array mounted onto the head of a biped humanoid robot. The LOCATA dataset consists of real-world recordings with ground-truth source locations provided by an optical tracking system. The size of the recording room is $7.1\times 9.8\times 3$ m, with $T_{60}\approx 0.55$ s. Task 6 contains three sequences of a total duration of $188.4$ s and two moving speakers. In our experiments we used four coplanar microphones, namely #5, #8, #11, and #12. The online sound-source localization method [16] was used to provide DOA estimates at each STFT frame, using a Hamming window of length 16 ms, with 8 ms shifts. The approach in  [16] requires a threshold, set to $0.3$ in our case, to select the number of significant active source, observed source DOAs, and the associated confidence values (see [21, 16]). The birth threshold, $\tau_{0}$, is set to 0.5 (Section IV).

To evaluate the method quantitatively, the estimated source trajectories are compared with the ground-truth trajectories over audio-active frames. Ground-truth audio-active frames are obtained using the voice activity detection (VAD) method of [22]. The permutation problem between the detected trajectories and the ground-truth trajectories is solved by means of a greedy gating algorithm: the error between all possible pairs of estimated and ground-truth trajectories is evaluated. Minimum-error pairs are selected for further comparison. A DOA estimate that is 15*∘* away from the ground-truth is treated as a false alarm detection. Sources that are not associated with a trajectory correspond to missed detections. For performance evaluation, the percentage of MDs and false alarms are evaluated over voice-active frames. The mean absolute error (MAE) the error between ground-truth DOAs and estimated DOA over all the active frames of all the speakers.

The observation-to-source assignment posteriors and the DOAs confidence weights are used to estimate voice-active frames:

[TABLE]

where $D=2$ and $\delta=0.025$ is a VAD threshold. Once an active source is detected, we output its trajectory.

The MAEs, MDs and FAs values, averaged over all recordings, are summarized in Table I. We compared the proposed von Mises VEM algorithm (vM-VEM) with three multi-speaker trackers: the von Mises PHD filter (vM-PHD) [13] and two versions the multiple speaker tracker of [16] based on Gaussians models (GM). [16] uses a first-order dynamic model whose effect is to smooth the estimated trajectories. We compared with both first-order (GM-FO) and zero-order (GM-ZO) dynamics. The proposed vM-VEM tracker yields the lowest false alarm (FA) rate of $5.9\%$ and MAE of $2.6$, and the second lowest MD rate of $23.9\%$. The GM-FO variant of [16] yields an MD rate of $22.3\%$ since it uses velocity information to smooth the trajectories. This illustrates the advantage of the von-Mises distribution to model directional data (DOA). The proposed von-Mises model uses a zero-order dynamics; nevertheless it achieves performance comparable with the Gaussian model that uses first-order dynamics.

The results for recordings #1 and #2 in Task 6 are shown in Fig. 1, using a sampling rate of 12 Hz for plotting. Note that the PHD-based filter method [13] has two caveats. First, observation-to-source assignments cannot be estimated (unless a post-processing step is performed), and second, the estimated source trajectories are not smooth. This stays in contrast with the proposed method which explicitly represents assignments with discrete latent variables and estimates them iteratively with VEM. Moreover, the proposed method yields smooth trajectories similar with those estimated by [16] and quite close to the ground truth.

VI Conclusion

We proposed a multiple audio-source tracking method using the von Mises distribution and we inferred a tractable solver based on a variational approximation of the posterior filtering distribution. Unlike the wrapped Gaussian distribution, the von Mises distribution explicitly models the circular variables associated with audio-source localization and tracking based on source DOAs. Using the recently released LOCATA dataset, we empirically showed that the proposed method compares favorably with two recent methods.

Appendix A Derivation of the E-S step

In order to obtain the formulae for the E-S step, we start from its definition in (9):

[TABLE]

We now use the decomposition in (1) to write:

[TABLE]

Let us now develop the expectation:

[TABLE]

where \stackrel{{\scriptstyle\mbox{\boldmaths}_{t}}}{{=}} denotes the equality up to an additive constant that does not depend on \mbox{\boldmaths}_{t}. Such a constant would become a multiplicative constant after the exponentiation in (21), and therefore can be ignored.

By replacing the developed expectation together with (12) we obtain:

[TABLE]

which can be rewritten as:

[TABLE]

(23) is important since it demonstrates that the a posteriori pdf of \mbox{\boldmaths}_{t} is separable on $n$ and therefore independent for each speaker. In addition, it allows us to rewrite the a posteriori pdf for each speaker, i.e., of $s_{tn}$ as a von Mises distribution by using the harmonic addition theorem, thus obtaining

[TABLE]

with $\mu_{tn}$ and $\kappa_{tn}$ defined as in (14) and (15).

Appendix B Derivation of the E-Z step

Similarly to the previous section, and in order to obtain the closed-form solution of the E-Z step, we start from its definition in (8):

[TABLE]

and we use the decomposition in (1),

[TABLE]

Since both the observation likelihood and the prior distribution are separable on $z_{tm}$, we can write:

[TABLE]

proving that the a posteriori pdf is also separable on $m$.

We can thus analyze the posterior of each $z_{tm}$ separately, by computing $q(z_{tm}=n)$:

[TABLE]

Let us first compute the expectation for $n\neq 0$:

[TABLE]

where for the last line we used the following variable change $\bar{s}=s_{tn}-\mu_{tn}$ and the definition of $I_{1}$ and $A$.

The case $n=0$ is even easier since the observation distribution is a uniform: $\mathrm{E}_{q(s_{tn})}\log p(y_{tm}|s_{tn},z_{tm}=n)=\mathrm{E}_{q(s_{tn})}-\log 2\pi=-\log(2\pi)$.

By using the fact that the prior distribution on $z_{tm}$ is denoted by $p(z_{tm}=n)=\pi_{n}$, we can now write the a posteriori distribution as $q(z_{tm}=n)\propto\pi_{n}\beta_{tmn}$ with:

[TABLE]

thus leading to the results in (16) and (3).

Appendix C Derivation of the M step

In order to derive the M step, we need first to compute the $Q$ function in (10),

[TABLE]

where each parameter is show below the corresponding term of the $Q$ function. Let us develop each term separately.

C-A Optimizing $\kappa_{y}$

[TABLE]

and by taking the derivative with respect to $\kappa_{y}$ we obtain:

[TABLE]

which corresponds to what was announced in the manuscript.

C-B Optimizing $\pi_{n}$’s

[TABLE]

This is the same formulae that is correct for any mixture model, and therefore the solution is standard and corresponds to the one reported in the manuscript.

C-C Optimizing $\kappa_{d}$

[TABLE]

where the dependency on $\kappa_{d}$ is implicit in $\tilde{\kappa}_{t-1,n}=A^{-1}(A(\kappa_{t-1,n})A(\kappa_{d}))$.

By taking the derivative with respect to $\kappa_{d}$ we obtain:

[TABLE]

with

[TABLE]

where $\tilde{A}(a)=\textrm{d}A^{-1}(a)/\textrm{d}a=(2-a^{2}+a^{4})/(1-a^{2})^{2}$.

By denoting the previous derivative as $B(\kappa_{d})=\frac{\partial\tilde{\kappa}_{t-1,n}}{\partial\kappa_{d}}$, we obtain the expression in the manuscript.

Appendix D Derivation of the birth probability

In this section we derive the expression for $\tau_{j}$ by computing the integral (17). Using the probabilistic model defined, we can write (the index $j$ is omitted):

[TABLE]

We will first marginalize $s_{t-L}$. To do that, we notice that if $p(s_{t-L})$ follows a von Mises with mean $\hat{\mu}_{t-L}$ and concentration $\hat{\kappa}_{t-L}$, then we can write:

[TABLE]

with

[TABLE]

where we used the harmonic addition theorem.

Now we can effectively compute the marginalization. The two terms involving $s_{t-L}$ are:

[TABLE]

with

[TABLE]

Therefore, the marginalization with respect to $s_{t-L}$ yields the following result:

[TABLE]

Since we have already seen that $p(s_{t-L+1})$ is also a von Mises distribution, we can use the same reasoning to marginalize with respecto to $s_{t-L+1}$. This strategy yields to the recursion presented in the main text.

Appendix E Results with errors