Charles Blandin
email: charles.blandin@inria.fr

Alexey Ozerov
email: alexey.ozerov@inria.fr

Emmanuel Vincent
email: emmanuel.vincent@inria.fr

Multi-source TDOA estimation in reverberant audio using angular spectra and clustering 1

Keywords: Multiple source localization, TDOA estimation, angular spectrum, clustering

We consider the problem of estimating the time differences of arrival (TDOAs) of multiple sources from a two-channel reverberant audio signal. While several clustering-based or angular spectrum-based methods have been proposed in the literature, only relatively small-scale experimental evaluations restricted to either category of methods have been carried out so far. We design and conduct the first large-scale experimental evaluation of these methods and investigate a two-step procedure combining angular spectra and clustering. In addition, we introduce and evaluate five new TDOA estimation methods inspired from signal-to-noise-ratio (SNR) weighting and probabilistic multi-source modeling techniques that have been successful for anechoic TDOA estimation and audio source separation. The results show that clustering-based methods do not improve upon angular spectrum-based methods. For 5 cm microphone spacing, the best TDOA estimation performance is achieved by one of the proposed SNR-based angular spectrum methods. For larger spacing, a variant of the generalized cross-correlation with phase transform (GCC-PHAT) method performs best.

Introduction

Recorded audio signals often result from the mixture of several sound sources. The problem of source localization consists of estimating the spatial positions of the sources and has many applications such as video-conferencing, surveillance, or source separation [1,[START_REF]Blind Speech Separation[END_REF]. When the signal is recorded by an array of sensors, this problem is often addressed by estimating the Time Difference Of Arrival (TDOA) of each source for each pair of sensors [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF]. Localization performance then mostly depends on the accuracy of the estimated TDOAs [1]. In the following, we focus on TDOA estimation of two or more sources for a given pair of sensors in a reverberant environment. The estimation of the number of sources and the assessment of the resulting localization performance for different array geometries are left for later study.

Most TDOA estimation methods use the Short Time Fourier Transform (STFT) of the signal [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF][START_REF] Yılmaz | Blind separation of speech mixtures via timefrequency masking[END_REF][START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]. Let us denote by x(t, f ) = [x 1 (t, f ), x 2 (t, f )] T and s n (t, f ), n = 1, . . . , N , the STFTs of the observed signals and the n-th source signal, where t = 1, . . . , T and f = 1, . . . , F are, respectively, time frame and frequency bin indices. With these notations, the source mixing process can be modeled as [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF] x(t, f ) =

N n=1 d(f, τ n )s n (t, f ) + b(t, f ), (1)
where

d(f, τ n ) = [1, e -2iπf τn ] T (2)
is the so-called steering vector associated with the n-th source of TDOA τ n (in seconds), and b(t, f ) models the reverberant part of the signal and additive noise, if any. We distinguish three general approaches for TDOA estimation. The simplest one consists of computing the TDOA of the mixture signal locally in each time-frequency bin and localizing the peaks of the resulting histogram [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF][START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF][START_REF] Mouba | A source localization/separation/respatialization system based on unsupervised classification of interaural cues[END_REF][START_REF] Loesch | Source number estimation and clustering for underdetermined blind source separation[END_REF]. This approach is restricted to closely spaced or binaural microphones, since in the case of widely spaced microphones local TDOA computation becomes ambiguous due to spatial aliasing [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF]. The second approach consists of iteratively estimating the time-frequency bins associated to each source and the corresponding TDOAs by means of some clustering algorithm [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF][START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF][START_REF] Mandel | An EM algorithm for localizing multiple sound sources in reverberant environments[END_REF]. This approach can be used for any microphone spacing but is sensitive to the initialization of the parameters (i.e., clusters and TDOAs). The third approach [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF][START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF][START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF][START_REF] Chami | A phase-based dual microphone method to count and locate audio sources in reverberant rooms[END_REF][START_REF] Liu | Localization of multiple sound sources with two microphones[END_REF] consists of building for each time-frequency bin a function of TDOA that is likely to exhibit a large value for true TDOAs, and pooling it across the time-frequency plane so as to obtain a so-called angular spectrum. The source TDOAs are then estimated as the highest peaks of this angular spectrum. This approach works for any microphone spacing and does not need any initial guess of the TDOAs. However, the possible presence of secondary peaks can alter the estimation of the TDOAs. In this paper, we focus on the two latter approaches, namely clustering and angular spectrum, which are always applicable.

Two main gaps are found in the literature. First, only relatively smallscale experimental evaluations restricted to either approach have been carried out so far [START_REF] Brutti | Comparison between different sound source localization techniques based on a real data collection[END_REF][START_REF] Loesch | Comparison of different algorithms for acoustic source localization[END_REF][START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF] and little is known about the variation of performance with respect to the microphone spacing, the reverberation time, the number of sources, the chosen angular spectrum pooling function or the chosen clustering initialization method. Second, existing methods typically rely on the assumption that the sources are disjoint in the time-frequency plane and affect the same weight to all time-frequency bins, regardless of the fact that the associated spatial information is less accurate in the presence of overlapping sources or reverberation.

This article aims to fill in these gaps. We design and conduct a largescale evaluation of angular spectrum-based and clustering-based methods on 1482 different configurations and investigate the use of the former for the initialization of the latter. In addition, we introduce and evaluate five new TDOA estimation methods inspired from signal-to-noise ratio (SNR) weighting or probabilistic modeling techniques that have been successful for anechoic TDOA estimation [START_REF] Arberet | A robust method to count and locate audio sources in a multichannel underdetermined mixture[END_REF][START_REF] Deville | Temporal and time-frequency correlation-based blind source separation methods. Part I: Determined and underdetermined linear instantaneous mixtures[END_REF][START_REF] Bofill | Identifying single source data for mixing matrix estimation in instantaneous blind source separation[END_REF], histogram-based reverberant TDOA estimation [START_REF] Loesch | Source number estimation and clustering for underdetermined blind source separation[END_REF] or audio source separation [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF][START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], but have not yet been explored for angular spectrum-based or clustering-based reverberant TDOA estimation. The proposed methods account for the presence of diffuse noise or interfering sources in each time-frequency bin and rely prioritarily on the time-frequency bins resulting from the direct sound of a single source. The code of all methods and the experimental data are available at http://bass-db.gforge.inria.fr/bss locate/.

The rest of the paper is organized as follows. In section 2, a short review of existing angular spectrum-based methods is given and the proposed SNR-based methods are introduced. Existing and proposed clustering-based methods are presented in section 3. The experimental evaluation is detailed in section 4, and conclusions are drawn in section 5.

Angular spectrum-based methods

The principle of angular spectrum-based methods is to construct a function φ(τ ) of TDOA τ whose peaks indicate the TDOAs of the sources. This is commonly achieved as follows. A local angular spectrum (or coherence measure) φ(t, f, τ ) is computed in each time-frequency bin (t, f ) for all discrete values of τ lying on a uniform grid in the range of possible TDOAs. This function is chosen so that it is likely to exhibit large values for the TDOAs of the sources which are active in this time-frequency bin. In order to robustify the estimation process and to overcome the spatial aliasing ambiguity occuring at high frequencies, the function φ(t, f, τ ) is summed over all frequencies. Then, it is reduced to a single dimension to obtain the angular spectrum from which the TDOAs are estimated. This is typically done by summing it over all time frames [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF]:

φ sum (τ ) = T t=1 F f =1 φ(t, f, τ ). (3)
A limitation of this approach is that it makes it difficult to localize a source that is active only within few time frames due to the integration of irrelevant information when the source is inactive. This can be addressed by taking the maximum over all time frames instead [START_REF] Chami | A phase-based dual microphone method to count and locate audio sources in reverberant rooms[END_REF]:

φ max (τ ) = max t F f =1 φ(t, f, τ ). (4)
Existing methods differ by the definition of the local angular spectrum function φ(t, f, τ ) and the choice of the pooling function, i.e., "sum" (3) or "max" (4). Alternative pooling functions consisting of taking the p-th largest value over all time frames [START_REF] Loesch | Adaptive segmentation and separation of determined convolutive mixtures under dynamic conditions[END_REF] or the mean of the p largest values decreased the average performance compared to the "max" in our experiments and are not considered hereafter.

Popular existing methods

Existing methods typically extract the spatial information in time-frequency bin (t, f ) from the empirical covariance matrix R xx (t, f ) of the input signal, which can be computed in the neighborhood of each time-frequency bin (t, f ) as [START_REF] Vincent | Underdetermined instantaneous audio source separation via local Gaussian modeling[END_REF]

R xx (t, f ) = t ′ ,f ′ w(t ′ -t, f ′ -f )x(t ′ , f ′ )x(t ′ , f ′ ) H t ′ ,f ′ w(t ′ -t, f ′ -f ) , ( 5
)
where w is a time-frequency windowing function of length L f × L t defining the size and the shape of the neighborhood, and (•) H denotes the Hermitian transposition operator.

The generalized cross-correlation with phase transform (GCC-PHAT) method [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF] is certainly the most popular angular spectrum-based method. Based on the assumption that the direct sound of one source predominates in each time-frequency bin, the TDOA of this source τ is estimated from the phase difference between the two channels represented by the argument of R xx (t, f ) 1,2 2 . Indeed, this phase difference is expected to be close to 2πf τ modulo 2π. The local angular spectrum is then defined as

φ GCC (t, f, τ ) = ℜ R xx (t, f ) 1,2 | R xx (t, f ) 1,2 | e -2iπf τ , (6)
where ℜ(z) denotes the real part of a complex number z. In [START_REF] Loesch | Adaptive segmentation and separation of determined convolutive mixtures under dynamic conditions[END_REF], it is proposed to use a non-linear function ρ of φ GCC (t, f, τ ) defined by ρ(u) = 1tanh(α √ 1u) in order to emphasize the large values of φ GCC (t, f, τ ).

In multiple signal classification (MUSIC) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], under the same assumption that one source is predominant in each time-frequency bin (t, f ), the local angular spectrum is computed as a measure of fit between the steering vector d(f, τ ) and the first principal component v(t, f ) of R xx (t, f ):

φ MUSIC (t, f, τ ) = 1 - 1 2 d(f, τ ) H v(t, f ) 2 -1 , (7)
where d(f, τ ) is defined by (2).

Nesta et al. [START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF] propose a method relaxing the assumption of one predominant source in each time-frequency bin. The time-frequency plane is split into time-frequency blocks, and it is assumed that there are at most two predominant sources in each block. Then, Independent Component Analysis (ICA) is applied in each time-frequency block (t, f ) to obtain two amplitudenormalized mixing coefficients r 1 (t, f ) and r 2 (t, f ) (see [START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF] for details) that are likely to be close to e -2iπf τ 1 and e -2iπf τ 2 up to a permutation, where τ 1 and τ 2 are the TDOAs of the two predominant sources in the considered block. The local angular spectrum of this method called cumulative state coherence transform (cSCT) is given by

φ cSCT (t, f, τ ) = 2 j=1 ρ 1 2 e -2iπf τ -r j (t, f ) , (8)
where ρ(u) = 1tanh(αu).

Proposed methods

All existing methods, except the cSCT, use the assumption that in each time-frequency bin, only one source is predominant. For cSCT, this assumption is replaced by that of at most two predominant sources. Both assumptions do not hold exactly for most audio data [START_REF] Vincent | Underdetermined instantaneous audio source separation via local Gaussian modeling[END_REF][START_REF] Arberet | Blind spectral-GMM estimation for underdetermined instantaneous audio source separation[END_REF]. Thus, some of the estimated local angular spectra φ(t, f, τ ) do not represent any "true" TDOA possibly leading to poor estimation of the TDOAs. A bounded timefrequency weighting function based on interchannel correlation has been proposed in [START_REF] Roman | Speech segregation based on sound localization[END_REF] to give more weight to the time-frequency bins involving a single predominant source. However, this function leads to overestimate the weight at low frequencies where interchannel correlation is large regardless of the number of active sources. Inspired by [START_REF] Arberet | A robust method to count and locate audio sources in a multichannel underdetermined mixture[END_REF], where it was done for anechoic mixtures, we propose to use the SNR as an unbounded measure to determine whether the information contained in a time-frequency bin results from a single source. We propose three methods to estimate the SNR below.

SNR estimation by beamforming

In each time-frequency bin, we define the SNR in the direction corresponding to the TDOA τ by the ratio between the sound power in this direction and the residual power. We estimate the power in the direction corresponding to the TDOA τ by the Minimum Variance Distortionless Response (MVDR) beamformer [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF]:

P (t, f, τ ) = d(f, τ ) H R xx (t, f ) -1 d(f, τ ) -1 , (9)
where R xx (t, f ) and d(f, τ ) are computed by ( 5) and (2), respectively. Then, we compute the residual power by subtracting the estimated power in the direction from an estimate of the total power: 1 2 tr R xx (t, f ) -P (t, f, τ ). Finally, we define the SNR in this direction as Our preliminary experiments showed that this function overestimates the SNR at low frequencies, where phase differences are small, regardless of the number of active sources. As a consequence, the values of φ MVDR (t, f, τ ) at low frequencies mask the values at higher frequencies (see Fig. 1, left).

φ MVDR (t, f, τ ) = d(f, τ ) H R xx (t, f ) -1 d(f, τ ) -1 1 2 tr R xx (t, f ) -d(f, τ ) H R xx (t, f ) -1 d(f, τ ) -1 . (10) φ MVDR (t=0, f, τ) f (kHz) -4 - 2

SNR estimation under a diffuse noise model

To address this problem, we propose to jointly estimate the sound power in direction τ and the power of the residual signal by using a diffuse noise model. We assume that in each time-frequency bin only one source s(t, f ) of TDOA τ is predominant, i.e., the mixing model (1) becomes

x(t, f ) = d(f, τ )s(t, f ) + b(t, f ). ( 11
)
Moreover, we assume that s(t, f ) and b(t, f ) follow independent zero-mean Gaussian distributions with, respectively, variance v s (t, f, τ ) and covariance

v b (t, f, τ )Ψ(f ). The variances v s (t, f, τ ), v b (t, f, τ ) > 0 represent the source
and noise power, and Ψ(f ) is the covariance matrix of a gain-normalized diffuse noise [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF][START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF]:

Ψ(f ) = 1 sinc(2πf d c ) sinc(2πf d c ) 1 , (12)
with d being the distance between the two microphones (in meters), c the speed of sound (in meters/sec.), and sinc(y) = sin(y) y the cardinal sine function. With these assumptions, it can be shown from ( 11) that the covariance matrix of the mixture in the time-frequency bin (t, f ) and for TDOA τ can be expressed as

R xx (t, f, τ ) = v s (t, f, τ )d(f, τ )d H (f, τ ) + v b (t, f, τ )Ψ(f ), (13)
and the log-likelihood of x(t, f ) can be written as:

log p(x(t, f )) = log N (x(t, f ); 0, R xx (t, f, τ )) -tr R -1 xx (t, f, τ ) R xx (t, f ) -log det (πR xx (t, f, τ )) . ( 14
)
Using a closed form solution from [START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF], we estimate v s (t, f, τ ) and v b (t, f, τ ) in the maximum likelihood sense as

v s (t, f, τ ) v b (t, f, τ ) = (diag(Λ 1 ) , diag(Λ 2 )) -1 diag(A -1 R xx (t, f )(A H ) -1 ), (15)
where (Y , Z) denotes the concatenation of matrices (or vectors) Y and Z, diag(Y) denotes the column vector of diagonal entries of matrix Y, A is the matrix whose columns are the eigenvectors of d(f, τ )d H (f, τ )Ψ -1 (f ), and 15) is negative [START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF]. Finally, the SNR in the time-frequency bin (t, f ) and for TDOA τ is computed as

Λ 1 , Λ 2 are equal respectively to A -1 d(f, τ )d H (f, τ )(A H ) -1 and A -1 Ψ(f )(A H ) -1 . Non-negativity of v s (t, f, τ ) and v b (t, f, τ ) is imposed by setting v s (t, f, τ ) to zero and v b (t, f, τ ) to 1 2 tr Ψ -1 (f ) R xx (t, f ) when v s (t, f, τ ) or v b (t, f, τ ) resulting from (
φ DNM (t, f, τ ) = v s (t, f, τ ) v b (t, f, τ ) . ( 16
)
This method effectively addresses the SNR overestimation problem at low frequencies but results in slightly wider peaks (see Fig. 1, center).

SNR estimation by frequency weighted beamforming

To take advantage of both the precision of φ MVDR (see Fig. 1, left) and the global shape of φ DNM (see Fig. 1, center), we combine these two methods by expressing a relationship between them, assuming that the input signal consists of a single source of TDOA τ = 0 and diffuse noise. In other words, we replace in [START_REF] Loesch | Source number estimation and clustering for underdetermined blind source separation[END_REF] the empirical covariance matrix R xx (t, f ) by the covariance matrix R xx (t, f, τ ) in ( 13), and we get

φ MVDR (t, f, τ ) = 1 + 2v s (t, f, τ )/v b (t, f, τ ) + sinc(2πf d c ) 1 -sinc(2πf d c ) . ( 17
)
By inverting [START_REF] Loesch | Comparison of different algorithms for acoustic source localization[END_REF], we estimate the SNR as:

φ MVDRW (t, f, τ ) = v s (t, f, τ ) v b (t, f, τ ) = w d (f )φ MVDR (t, f, τ ) + w d (f ) -1, ( 18
)
where

w d (f ) = 1 2 (1 -sinc(2πf d c )
) is a weighting factor depending on the frequency and the microphone spacing d. This factor reduces the impact of low frequencies, as it can be seen from Figure 2 and from Figure 1 (top right) below 1 KHz. It can be also noted from Figure 1 that φ MVDRW has sharper peaks than both φ MVDR and φ DNM .

Clustering-based methods

In angular spectrum-based methods a measure of source activity (e.g., the SNR) can be exploited, but the estimation of this measure does not rely on any estimates of the TDOAs. In contrast, the principle of clusteringbased methods consists of estimating this measure given some estimates of the TDOAs, then reestimating the TDOAs relying on this measure, and so on. This results in iterating between the following two steps:

• re-estimate the contribution of each source to each time-frequency bin (clusters), given current estimates of the TDOAs,

• re-estimate the TDOAs, given current estimates of the clusters.

Popular existing methods

Sawada et al. [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF] propose to perform the clustering in a hard manner, by associating each time-frequency bin to the closest source n tf in the sense of the Euclidean distance between its steering vector and the phase and amplitude-normalized observation. Each TDOA τ n is then reestimated from the time-frequency bins of the corresponding cluster only.

It is also possible to perform the clustering in a soft manner in a probabilistic setting. Araki et al. [START_REF] Araki | Stereo source separation and source counting with MAP estimation with Dirichlet prior considering spatial aliasing problem[END_REF] assume that the phase difference between the two channels follows a wrapped Gaussian distribution for each source and each time-frequency bin. This results in a Gaussian mixture model (GMM), whose parameters are estimated via the Expectation-Maximization (EM) algorithm. The E-step achieves soft clustering by estimating the probability of n tf . This method also includes a sparsity-inducing prior on the GMM mixture weights enabling the estimation of the number of sources.

Pointing that the Euclidean distance and the Gaussian distribution over the phase differences do not accurately represent the effect of multiple sources and reverberation, Izumi et al. [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF] adopt the following model instead:

x(t, f ) = s n tf (t, f )d(f, τ n tf ) + b(t, f ). ( 19
)
The source coefficients s n tf (t, f ) are considered as deterministic parameters and the noise term b(t, f ) is assumed to be diffuse and modeled as a zeromean Gaussian random vector of covariance matrix v b Ψ(f ), where Ψ(f ) is given by ( 12) and v b is a constant parameter. Parameter estimation is achieved via the EM algorithm again.

Proposed methods 3.2.1. EM algorithm with one predominant source in each time-frequency bin

In line with [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF], we assume that in each time-frequency bin (t, f ), there is only one predominant source n tf and a diffuse noise b(t, f ). However, in contrast to [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF], we assume that both s n tf (t, f ) and b(t, f ) follow independent zero-mean Gaussian distributions with covariances v s n tf (t, f ) and v b n tf (t, f )Ψ(f ), where v s n tf (t, f ) and v b n tf (t, f ) represent respectively the source and the noise variances, and Ψ(f ) is defined by [START_REF] Mandel | An EM algorithm for localizing multiple sound sources in reverberant environments[END_REF]. Given the predominant source index n tf , the mixture covariance matrix is

R xx,n tf (t, f ) = Ξ x (f, τ n tf , v s n tf (t, f ), v b n tf (t, f )),
where

Ξ x (f, τ, v s , v b ) v s d(f, τ )d H (f, τ ) + v b Ψ(f ). ( 20
)
The set of parameters to be estimated is

θ = {τ n } N n=1 , v s n (t, f ), v b n (t, f ) N,T,F n,t,f =1 . ( 21
)
Under the above assumptions and assuming a uniform prior over the source indices, the observation x(t, f ) follows the Gaussian mixture model (GMM) distribution

p(x(t, f )|θ) = n 1 N p(x(t, f )|τ n , v s n (t, f ), v b n (t, f )), (22)
where

p(x(t, f )|τ, v s , v b ) = N (x(t, f ); 0, Ξ x (f, τ, v s , v b )), (23)
and N (x(t, f ); 0, Ξ x (f, τ, v s , v b )) is defined by [START_REF] Chami | A phase-based dual microphone method to count and locate audio sources in reverberant rooms[END_REF].

We use an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] to estimate the parameters θ in the maximum likelihood sense, considering the set of predominant source indices {n tf } t,f as latent data. The updates of the resulting algorithm and some hints for its derivation are given in Appendix A.

EM algorithm with multiple sources in each time-frequency bin

As it was already mentioned in section 2.2, the assumption of only one predominant source does not hold exactly for most audio data. Several works [START_REF] Vincent | Underdetermined instantaneous audio source separation via local Gaussian modeling[END_REF][START_REF] Arberet | Blind spectral-GMM estimation for underdetermined instantaneous audio source separation[END_REF][START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF][START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] have shown that relaxing this assumption can be very beneficial for audio source separation. Thus, we here investigate whether such an approach could be beneficial for multi-source localization. The model presented below is mostly inspired by the models proposed in [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF][START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]. We consider that all sources can be present in each time-frequency bin and model the mixing process as

x(t, f ) = D(f, τ )s(t, f ) + b(t, f ), (24)
where τ = [τ 1 , . . . , τ N ] is the vector of TDOAs,

D(f, τ ) = [d(f, τ 1 ), . . . , d(f, τ N )], s(t, f ) = [s 1 (t, f ), . . . , s N (t, f )] T and b(t, f ) = [b 1 (t, f ), b 2 (t, f )] T .
We assume that b(t, f ) and s(t, f ) follow zero-mean Gaussian distributions with covariance matrices respectively equal to v b (t, f )Ψ(f ) and R ss (t, f ), where R ss (t, f ) is the diagonal matrix with n-th diagonal entry being equal to v s n (t, f ). With these assumptions the observation x(t, f ) follows a zero-mean Gaussian distribution with covariance matrix

R xx (t, f ) = D(f, τ )R ss (t, f )D H (f, τ ) + v b (t, f )Ψ(f ). ( 25
)
The set of parameters to be estimated is

θ = {τ n } N n=1 , {v s n (t, f )} N,T,F n,t,f =1 , v b (t, f ) T,F t,f . (26)
To estimate the parameters in the maximum likelihood sense we derive an EM algorithm, considering the sources {s(t, f )} t,f as latent data. The algorithm is summarized in Appendix B. We consider four variants of the algorithm:

1. "EM-multi (TF noise)": the noise variances v b (t, f ) are unconstrained, 2. "EM-multi (F noise)": in line with [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], the noise variances v b (t, f ) are constrained to be constant over time, i.e., v b (t, f ) = v b (f ), 3. "EM-multi (const noise)": in line with [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF], the noise variances v b (t, f ) are constrained to be constant over time and frequency, i.e., v b (t, f ) = v b , 4. "EM-multi (no noise)": in line with [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF], the noise variances v b (t, f ) are fixed to a small positive value ε fix .

Evaluation

We performed a large-scale evaluation of the methods presented in this article. The data, the evaluation measures and the code of the tested methods are available on http://bass-db.gforge.inria.fr/bss locate/.

Data

The experimental evaluation was carried out on a large number of mixtures of male speech, female speech and music sources taken from the database of the 2008 Signal Separation Evaluation Campaign (SiSEC) [START_REF] Vincent | The 2008 signal separation evaluation campaign: A community-based approach to large-scale evaluation[END_REF] "underdetermined speech and music mixtures" task. Mixing filters were simulated with the Roomsimove Toolbox3 for a rectangular room of dimensions 4.45 m × 3.55 m × 2.5 m and omnidirectional microphones. We considered all possible combinations of the following parameters:

• Number of sources N : from 2 to 6.

• Reverberation time RT 604 : 50 ms, 100 ms, 150 ms, 250 ms, 500 ms, 750 ms.

• Microphone spacing d: 5 cm, 15 cm, 30 cm, 1 m.

• Distance between the sources and the center of the microphone pair: 20 cm, 50 cm, 1 m, 2 m.

• Angular position of the sources: between 3 and 5 randomly generated scenarios according to the number of sources, i.e., 5 scenarios for 2 sources, 4 for 3 -4 sources, and 3 for 5 -6 sources. Moreover, the scenarios were generated with the following restrictions: the Directions Of Arrival (DOAs) cannot be smaller (greater) than 30 degrees (150 degrees) and the absolute difference between DOAs of any pair of sources cannot be smaller than 15 degrees for 2-5 sources and 10 degrees for 6 sources.

• Three source types (female speech, male speech and music).

We only kept situations for which the distance between the sources and the microphone pair is lower than 0.8d to ensure that the far-field assumption holds, so that the relation between a TDOA τ and its corresponding DOA η (in degrees) can be expressed by

τ = d c cos 2π η 360 . ( 27
)
This resulted in a total of 4446 mixtures. All the signals were sampled at 16 kHz.

Evaluation measures 4.2.1. Localization

Evaluation in terms of localization consists of measuring the capacity of the methods to estimate the true TDOAs with some tolerance. Since we do not attempt to estimate the number of sources N , we apply the considered methods for all possible numbers of sources J from 1 to 20. For a given J, we select the J highest peaks of the angular spectra and run the clustering algorithms with J clusters. We then evaluate the list of J estimated TDOAs compared to the list of N true TDOAs in terms of recall, precision and Fmeasure. An estimated TDOA τ is considered as correctly estimated if it is close enough to the true TDOA τ , in the sense that the corresponding DOAs η and η, computed by inverting [START_REF] Roman | Speech segregation based on sound localization[END_REF], differ by less than 5 degrees modulo 180 degrees. Denoting by I J the number of correctly estimated TDOAs, recall, precision and F-measure are defined by [START_REF] Van Rijsbergen | Information retrieval[END_REF] Recall(J) =

I J N , (28)
Precision(J) = I J J , (29)
F-measure(J) = 2 × Recall(J) × Precision(J) Recall(J) + Precision(J) . ( 30
)
Figure 3 shows the average recall, precision and F-measure of several angular spectrum-based methods over all mixtures with N = 3 sources. The results are represented as a function of the assumed number of sources J. As expected, the recall increases and the precision decreases with J, while the F-measure is maximum for some J = J opt . J opt equals to the number of sources N for some method, but not for all. Thus, J opt depends on the method A and on the number of sources N . We observed a similar behavior for other methods and other numbers of sources. We also noticed that it depends on the microphone spacing d but that it is insensitive to the other parameters. Thus, for each method A, each number of sources N and each microphone spacing d, we compute J opt (A, N, d) that maximizes the average F-measure over the corresponding mixtures. In the rest of this section, we assume that J is fixed to J opt (A, N, d) for all mixtures.

Accuracy

When the TDOAs are correctly estimated, one would like to know how accurate these estimates are. Let the TDOAs of all the mixtures considered be enumerated through as {τ i } I i=1 . Let I ⊂ {1, . . . , I} denote the subset of TDOAs that are correctly estimated by all the methods under comparison. We define the average accuracy as follows:

Accuracy = 1 i∈I N -1 i i∈I N -1 i (η i -ηi ) mod 180 , (31)
where η i and ηi are the DOAs computed from the true τ i and estimated τ i TDOAs using [START_REF] Roman | Speech segregation based on sound localization[END_REF], and N i denotes the number of sources in the corresponding mixture. Note that this measures the accuracy of only those TDOAs which are correctly estimated by all methods under consideration instead of each method individually.

Parameters

All the methods evaluated below had exactly the same front-end to compute the empirical covariance matrices R xx (t, f ) (5), except the cSCT [START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF], for which we used an implementation provided by the author of [START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF]. The STFT was computed with half-overlapping sine windows of length 1024. The time-frequency windowing function w in (5) was the outer product of two Hanning windows. Its size was set to L f = L t = 1 for GCC-PHAT and GCC-NONLIN, L f = 15 and L t = 3 for the other angular spectrum-based methods and L f = L t = 3 for all clustering-based methods, since this gave the best results in our preliminary experiments.

Angular spectrum-based methods

We consider four state-of-the-art local angular spectra: GCC-PHAT [START_REF] Knapp | The generalized cross-correlation method for estimation of time delay[END_REF], a version of GCC-PHAT with a non-linear function [START_REF] Loesch | Adaptive segmentation and separation of determined convolutive mixtures under dynamic conditions[END_REF] (denoted as GCC-NONLIN), MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] and cSCT [START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF], and the three proposed SNRbased local angular spectra: MVDR, DNM and MVDRW (see Sec. 2.2). All seven local angular spectra were evaluated with both φ sum (•) (3) and φ max (•) (4) pooling functions, except the cSCT [START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF] 5 . Note also that for all the methods, in line with [START_REF] Loesch | Adaptive segmentation and separation of determined convolutive mixtures under dynamic conditions[END_REF], we do not select peaks that differ by less than 5 degrees from a higher peak.

The results in terms of average recall, precision, F-measure and accuracy are reported in Table 1. Since the average accuracy is computed over the TDOAs that are correctly estimated by all the compared methods, as explained in section 4.2.2, we have selected for this comparison the methods that performed better than 0.7 in terms of average F-measure. This was done in order to avoid a comparison over a very small subset of TDOAs.

Analyzing the results in terms of F-measure, we see that for the "sum" pooling function MVDR and MVDRW outperform all other methods, which confirms our findings in [START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF]. Using the "max" pooling function instead of "sum" improves MVDR and MVDRW. However, it improves even more GCC-PHAT and GCC-NONLIN, which perform best in the end. The results in terms of average accuracy appear correlated with the average F-measure, and GCC-PHATmax is the most accurate on average. We then investigate the behaviour of the average F-measure as a function of the microphone spacing, the reverberation time and the number of sources. For this comparison we chose five methods among the best ones, namely GCC-NONLINsum, GCC-NONLINmax, cSCTsum, MVDRWsum and MVDRWmax. The results are shown in Figure 4. While the reverberation time and the number of sources have little influence on the ranking of different methods, the microphone spacing has. Note that in case of the "sum" pooling function the advantage of MVDRW against GCC-NONLIN was essentially due to its better performance for the 5 cm microphone spacing. Using the "max" pooling function improves GCC-NONLIN for all microphone spacings, but it improves MVDRW only for small spacings, while leading to a big performance degradation for large spacings. However, our proposed MVDRWmax method outperforms GCC-NONLINmax by 0.03 Fmeasure for 5 cm microphone spacing. This setting is useful in practice, since it matches the size of handheld devices such as portable audio recorders and smartphones.

Clustering-based methods

The evaluation of clustering-based methods, as compared to that of angular spectrum-based methods, is much more computationally expensive within our evaluation framework. Indeed, as explained in Sec. 4.2.1, each clustering-based method must be run 20 times (for J from 1 to 20) with several EM iterations. Since we have not noticed any influence of signal type on the ranking of angular spectrum-based methods, we retain only the speech mixtures for this evaluation. We consider two state-of-the-art methods by Sawada et al. [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF] and Izumi et al. [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF], and the four of the proposed methods: the "EM-predom" method introduced in Sec. 3.2.1 and the variants of the "EM-multi" method introduced in Sec. 3.2.2. To evaluate the impact of initialization on the each method was run twice: (i) with randomly initialized TDOAs and (ii) with initial TDOAs estimated by GCC-NONLINmax, i.e., the best angular spectrum-based method. For each mixture, the random initial TDOAs were drawn from the uniform distribution on the interval [-d/c, d/c]. For each method the source and the noise variances v s n (t, f ) and v b n (t, f ) were initialized by some constant values v s,init and v b,init adjusted during preliminary tests. Each method was run for a maximum of 100 iterations until the estimated TDOAs did not change from one iteration to the next.

The results in terms of average evaluation measures, including those of the initialization (noted by "None"), are summarized in Table 2. All methods fail when randomly initialized: the best results achieved by the proposed "EM-predom" method remain far below those of GCC-NONLINmax and all other methods except that Sawada et al. do not improve over mere random guess. The initialization provided by GCC-NONLINmax greatly improves the average F-measure for all methods. However, GCC-NONLINmax alone is the most accurate in the end: the four "EM-multi" methods provide exactly the same results, while Sawada et al. and "EM-predom" result in a small performance degradation. This illustrates the fact that these methods remain stuck in the local optimum provided by GCC-NONLINmax, so that improved clustering techniques are needed to overcome this performance ceiling.

We choose the four distinct methods after proper initialization, namely GCC-NONLINmax, Sawada et al. [START_REF] Sawada | Grouping separated frequency components by estimating propagation model parameters in frequency-domain blind source separation[END_REF], Izumi et al. [START_REF] Izumi | Sparseness-based 2ch BSS using the EM algorithm in reverberant environment[END_REF] and "EM-predom", to investigate the behaviour of the average F-measure as a function of different parameters. The results are shown in Figure 5. We see that the clustering-based methods, as compared to GCC-NONLINmax used for the initialization, decrease the average F-measure for all considered conditions. Moreover, the drop in performance is mostly pronounced for small microphone spacings.

Conclusion

In this paper we introduced several multi-source TDOA estimation methods based on angular spectra and clustering. The common motivation behind all these methods is to go beyond the assumption of a single predominant source per time-frequency bin. This is achieved by either using the SNR as an unbounded measure of source activity or by exploiting directly the multi-source hypothesis within a probabilistic model. We evaluated the proposed and the five most popular state-of-the-art methods on 1482 different configurations. To our knowledge, this is by far the largest-scale evaluation of multi-source TDOA estimation methods to date. Among all angular spectrum-based methods, the best performance was achieved by the proposed MVDRWmax method for 5 cm microphone spacing and by GCC-NONLINmax for larger spacings. The former setting is useful in practice, since it matches the size of handheld devices such as portable audio recorders and smartphones.

All clustering-based methods were evaluated with both random TDOA initialization and an initialization provided by the best angular spectrumbased method. First, we showed that it is very important to initialize the clustering algorithms with "good" TDOA estimates, since random initialization leads to really poor performance. Second, we observed that none of the clustering-based methods was able to improve the TDOA estimates compared to the best angular spectrum-based method.

Our findings show that SNR weighting and probabilistic modeling techniques that have been successful for anechoic TDOA estimation and audio source separation bring little or no improvement for reverberant TDOA estimation compared to GCC-PHAT, so that more specific approaches are needed to solve this problem in the future. The role of the "max" pooling function, its success and its potential limits should be better understood, so as to be able to propose more powerful non-linear pooling functions. The resulting localization performance should be assessed for different array geometries and the problem of estimating the number of sources remains open, with few methods proposed so far [START_REF] Araki | Stereo source separation and source counting with MAP estimation with Dirichlet prior considering spatial aliasing problem[END_REF][START_REF] Chami | A phase-based dual microphone method to count and locate audio sources in reverberant rooms[END_REF][START_REF] Loesch | Adaptive segmentation and separation of determined convolutive mixtures under dynamic conditions[END_REF]. Finally, we still believe in the potential of clustering-based methods provided that the considered i.i.d. source variance models are replaced by more structured audio-specific models as in, e.g., [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. An alternative way would be to try injecting the "max" pooling function within the probabilistic model behind clustering. assumptions of Section 3.2.2 is equal (up to some additive constant) to:

Q(θ, θ ′ ) c = - n,
G s (t, f ) = R ′ ss (t, f )D H (f, τ ′ )(R ′ xx (t, f )) -1 , (B.5) R ′ xx (t, f ) = D(f, τ ′ )R ′ ss (t, f )D H (f, τ ′ ) + v b′ (t, f )Ψ(f ). (B.6)
One iteration of the EM parameter updates optimizing (A.1) for the auxiliary function (B.1) consists of the following steps:

1. Compute R xx (t, f ), R xs (t, f ) and R ss (t, f ) using [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] with J (t, f ) ⊂ {1, . . . , T } × {1, . . . , F } denoting the subset of timefrequency bins where v b (t, f ) is considered constant and #J (t, f ) denoting the number of elements in this subset. This formulation allows the implemention of the first three variants of the algorithm mentioned at the end of section 3.2.2. To implement the fourth variant one just needs to skip the noise variance update. 5. Set θ ′ = θ.

Figure 1 :

 1 Figure 1: Local angular spectra φ(t, f, τ ) plotted at time t = 0 (top) and global angular spectra φ sum (τ ) pooled over all time frames (bottom) for φ MVDR (left), φ DNM (center) and φ MVDRW (right), computed from the mixture of three female speech sources placed at 50 cm from the center of the microphone pair, with d = 15 cm and a reverberation time of 500 ms. Dash-dotted lines indicate the true TDOAs.

Figure 2 :

 2 Figure 2: Weighting factor w d (f ) for microphone spacings d = 5 cm, 15 cm, 30 cm and 1 m.

Figure 3 :

 3 Figure 3: Average recall (left), precision (center) and F-measure (right) as functions of J for several angular spectrum-based methods and for all mixtures with N = 3 sources.

Figure 4 :

 4 Figure 4: Angular spectrum-based methods: average F-measure as a function of the microphone spacing d (left), the reverberation time RT60 (center), and the number of sources N (right).

Figure 5 :

 5 Figure 5: Clustering-based methods: average F-measure as a function of microphone spacing d (left), reverberation time RT60 (center), and number of sources N (right). All methods are initialized by the TDOAs estimated by GCC-NONLINmax.

Table 1 :

 1 Angular spectrum-based methods: average results for all 4446 mixtures. The average accuracy was computed over the 7755 TDOAs correctly estimated by the methods with average F-measure greater than 0.7 (from a total of 16614).

		Recall	Precision	F-measure	Accuracy
	Local angular	(0 to 1)	(0 to 1)	(0 to 1)	(degrees)
	spectrum function	φ sum	φ max	φ sum	φ max	φ sum	φ max	φ sum	φ max
	GCC-PHAT [3]	0.61	0.84	0.72	0.85	0.65	0.84	-	0.42
	GCC-NONLIN [24]	0.67	0.85	0.73	0.85	0.69	0.85	-	0.45
	MUSIC [5]	0.65	0.64	0.33	0.30	0.43	0.41	-	-
	cSCT [13]	0.67	0.82	0.74	0.83	0.70	0.82	0.83	0.53
	MVDR	0.73	0.80	0.75	0.80	0.74	0.80	0.56	0.48
	DNM	0.70	0.72	0.73	0.71	0.71	0.71	0.63	0.61
	MVDRW	0.73	0.80	0.75	0.80	0.74	0.80	0.56	0.48

Table 2 :

 2 Clustering-based methods: average results for 1482 female speech mixtures. Notations: "rand" means initializing clustering by random TDOAs, "init" means initializing it by TDOAs estimated with GCC-NONLINmax, and "None" means evaluating the initialization without running any clustering algorithm. The average accuracy was computed over the 4384 TDOAs correctly estimated by the methods with average F-measure greater than 0.7 (from a total of 5538).

		Recall	Precision	F-measure	Accuracy
	Clustering algorithm	(0 to 1)	(0 to 1)	(0 to 1)	(degrees)
		rand	init	rand	init	rand	init	rand	init
	None	0.53	0.90	0.24	0.90	0.32	0.90	-	0.56
	Sawada et al. [11]	0.66	0.88	0.76	0.87	0.49	0.87	-	0.75
	Izumi et al. [6]	0.29	0.30	0.48	0.51	0.30	0.35	-	-
	EM-predom	0.66	0.84	0.81	0.85	0.52	0.85	-	0.66
	EM-multi (TF noise)	0.94	0.90	0.50	0.90	0.32	0.90	-	0.56
	EM-multi (F noise)	0.94	0.90	0.50	0.90	0.32	0.90	-	0.56
	EM-multi (const noise)	0.94	0.90	0.48	0.90	0.31	0.90	-	0.56
	EM-multi (no noise)	0.94	0.90	0.47	0.90	0.32	0.90	-	0.56

  , f, τ , θ ′ ) tr Ψ -1 (f ) R xx (t, f ) -Ψ -1 (f )D(f, τ ) R H xs (t, f ) -Ψ -1 (f ) R xs (t, f )D H (f, τ ) + Ψ -1 (f )D(f, τ ) R ss (t, f )D H (f, τ ) , (B.2) R xs (t, f ) = R xx (t, f )G H s (t, f ), (B.3) R ss (t, f ) = G s (t, f ) R xs (t, f ) + I N -G s (t, f )D(f, τ ′ ) R ′ ss (t, f ), (B.4) with R xx (t, f ) defined by[START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], and

	t,f	log v s n (t, f ) +	R ss (t, f ) n,n v s n (t, f )	-	t,f	2 log v b (t, f ) +	M(t, f, τ , θ ′ ) v b (t, f )	,
							(B.1)
	where
	M(t

  , (B.3) and (B.4). 2. Update the source variances as v s n (t, f ) = R ss (t, f ) n,n . 3. Update the TDOAs as τ n = τ * n , where: ′ ) is defined by (B.2). More precisely, each TDOA τ n (n = 1, . . . , N ) is updated in turn, while keeping the other TDOAs {τ m } m =n fixed 7 .

	4. Update noise variances:
	v b (t, f ) =	1 2 • #J (t, f )	( t, f )∈J (t,f )	M( t, f , τ , θ ′ ),	(B.8)
	τ * n = arg min τn∈Γ	t,f	1 v b (t, f )	M(t, f, τ , θ ′ ),	(B.7)

and M(t, f, τ , θ

This work was supported in part by the ECHANGE project, funded by ANR, and by the Quaero Programme, funded by OSEO.

Here b Rxx(t, f )i,j denotes the (i, j)-th element of matrix b Rxx(t, f ).

E. Vincent, D. Campbell, Roomsimove, a Matlab toolbox for the computation of simulated room impulse reponses for moving sources, http://www.irisa.fr/metiss/members/ evincent/software.

[START_REF] Yılmaz | Blind separation of speech mixtures via timefrequency masking[END_REF] RT60 is the time required for reflections of a direct sound to decay by 60 dB below the level of the direct sound.

Since we have only an implementation of the cSCT[START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF] given us by the author of[START_REF] Nesta | Cumulative state coherence transform for a robust two-channel multiple source localization[END_REF] without the corresponding sources, we were not able injecting the "max" pooling function into this implementation.

Note that in practice this needs to be done only once for all iterations.

Both the usage of the noise variances v b (t, f ) in the updates of TDOAs (B.7) and the alternating nature of these updates do not guarantee the maximization of the auxiliary function, as in (A.1). However, they guarantee its non-decrease, i.e., Q(θ (l) , θ (l-1) ) ≥ Q(θ (l-1) , θ (l-1) ). Thus, the resulting algorithm is rather a Generalized EM (GEM) algorithm[START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Updating TDOAs jointly, instead of alternatively, is possible as well, but it is avoided here, since it would lead to a computational complexity growing exponentially with the number of sources.

Acknowledgment

The authors would like to thank Francesco Nesta for sharing his cSCT implementation, and Nobutaka Ito for discussion during the course of this work.

Appendix A. Updates for the EM algorithm with one predominant source in each time-frequency bin An EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is an iterative algorithm consisting in updating the parameters θ (l) at every iteration l as follows:

where Q(θ, θ ′ ) is the so-called auxiliary function. Under the assumptions of Section 3.2.1 this function is equal (up to some additive constant) to

where, for any parameter z from θ, z ′ denotes the corresponding parameter from θ ′ and p(x(t, f )|τ, v s , v b ) is defined by [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF].

As before, the TDOAs are only estimated on the discrete uniform grid in the range of possible TDOAs. Let us denote this grid by Γ. One iteration the EM parameter updates optimizing (A.1) for the auxiliary function (A.2) consists of the following steps:

1. Estimate vs n (t, f, τ ) and vb n (t, f, τ ) in the maximum likelihood sense for all time-frequency bins and all possible TDOAs τ ∈ Γ using (15) 6 . 2. Compute the posterior cluster probabilities γ n (t, f ) as in (A.4) and normalize them so that n γ n (t, f ) = 1. 3. Update the TDOAs τ n as follows:

Appendix B. Updates for the EM algorithm with multiple sources in each time-frequency bin

Performing some derivations analogous to those from [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] and [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], it can be shown that the auxiliary function Q(θ, θ ′ ) for the EM algorithm under


Benjamin Charlier
email: benjamin.charlier@math.univ-toulouse.fr

Necessary and sufficient condition for the existence of a Fréchet mean on the circle

Keywords: circular data, Fréchet mean, uniqueness. AMS classification: 62H11

Let (S 1 , d S 1 ) be the unit circle in R 2 endowed with the arclength distance. We give a sufficient and necessary condition for a general probability measure µ to admit a well defined Fréchet mean on (S 1 , d S 1 ). We derive a new sufficient condition of existence P (α, ϕ) with no restriction on the support of the measure. Then, we study the convergence of the empirical Fréchet mean to the Fréchet mean and we give an algorithm to compute it.

In many fields of interest, results of an experiment are objects taking values in non-Euclidean spaces. A rather general framework to model such data is Riemannian geometry and more particularly quotient manifolds. As an illustration, in biology or geology, directional data are often used, see e.g. [START_REF] Kanti | Directional statistics[END_REF] or [START_REF] Fisher | Statistical analysis of circular data[END_REF] and references therein. In this case, observations take their values in the circle or a sphere, that is, an Euclidean space quotiented by the action of scaling.

The usual definitions of basic statistical concepts were developed in an Euclidean framework. Therefore, these definitions must be adapted for random variables with values in non-Euclidean spaces such as manifolds. To describe the localization of a probability distribution, one needs to define a central value such as a mean or a median. There has been multiple attempts to give a definition of a mean in non Euclidean space, see among many others [START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds. i[END_REF][START_REF] Manuel Corcuera | Riemannian barycentres and geodesic convexity[END_REF][START_REF] Le | Locating Fréchet means with application to shape spaces[END_REF][START_REF] Le | Estimation of Riemannian barycentres[END_REF][START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF][START_REF] Émery | Sur le barycentre d'une probabilité dans une variété[END_REF][START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF] or [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF].

In this paper, we consider the so-called Fréchet mean, see [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF][START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF][START_REF] Le | Locating Fréchet means with application to shape spaces[END_REF] or [START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds. i[END_REF] and references therein. We are particularly interested in the study of its uniqueness. The Fréchet mean is defined on general metric spaces by extending the fact that Euclidean mean minimizes the sum of square the distance to the data, see equation (1.3) below. To study the well definiteness of the Fréchet mean on a manifold, two facts must be taken into account: non uniqueness of geodesics from one point to another (existence of a cut locus) and the effect of curvature, see e.g. [START_REF] Manuel Corcuera | Riemannian barycentres and geodesic convexity[END_REF] for further discussion. Due to the cut-locus, the distance function is no longer convex and finding conditions to ensure the uniqueness of the Fréchet mean is not obvious. Two main directions have been explored in the literature: bounding the support of the measure in [START_REF] Samuel | Spherical averages and applications to spherical splines and interpolation[END_REF] for the n-spheres and in [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF][START_REF] Kendall | Wiley Series in Probability and Statistics[END_REF][START_REF] Le | Locating Fréchet means with application to shape spaces[END_REF][START_REF] Afsari | Riemannian L p center of mass: existence, uniqueness, and convexity[END_REF] for manifolds, or consider special cases of absolutely continuous radial distributions, see [START_REF] Kaziska | The karcher mean of a class of symmetric distributions on the circle[END_REF] for the unit circle and or [START_REF] Kendall | Wiley Series in Probability and Statistics[END_REF][START_REF] Le | On the consistency of procrustean mean shapes[END_REF] for projective spaces. In a sense, these two conditions control the concentration of the probability measure. The philosophy behind these works is to ensure a convexity property of the Fréchet functional given by equation (1.2) below, see e.g. the introduction of [START_REF] Afsari | Riemannian L p center of mass: existence, uniqueness, and convexity[END_REF] for a review of the above cited papers.

Fréchet mean on the circle

A standard way to extend the definition of the Euclidean mean in non-Euclidean metric space is to use the minimization property of the Euclidean mean. This definition is usually credited to M. Fréchet in [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF] although some authors credit it to E. Cartan, see e.g. [START_REF] Kendall | Wiley Series in Probability and Statistics[END_REF]. Let S 1 be the unit circle of the plane, S 1 = {x 2 1 + x 2 2 = 1, (x 1 , x 2 ) ∈ R 2 }. endowed with the arclength distance given for all x = (x 1 , x 2 ), p = (p 1 , p 2 ) ∈ S 1 by

d S 1 (x, p) = 2 arcsin x -p 2 , (1.1)
where xp = (x 1p 1 ) 2 + (x 2p 2 ) 2 is the Euclidean norm in R 2 . In the whole paper, µ is a probability measure on S 1 and the Fréchet functional is defined for all p ∈ S 1 by

F µ (p) = 1 2 S 1 d 2 (x, p)dµ(x). (1.2)
The Fréchet functional is Lipschitz since by the triangle inequality we have

|F µ (p 1 ) -F µ (p 2 )| ≤ 2πd S 1 (p
F µ (p). (1.3)
The argmins of F µ are also called Riemannian center of mass [START_REF] Oller | Intrinsic analysis of statistical estimation[END_REF] or intrinsic mean [START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds. i[END_REF] as the (S 1 , d S 1 ) is a simple one dimensional compact Riemannian manifold. The advantage of dealing with a simple object such as the circle is that curvature problems disappear and we only face the cutlocus problem. In this sense, it allows us to completely understand its effect on the non-convexity of the distance function d S 1 , and to give a complete answer about the problem of uniqueness. In what follows, we fully characterize probability measures that admit a well defined Fréchet mean on the circle (S 1 , d S 1 ). In particularly a necessary and sufficient condition is given in Theorem 4.1, which links the existence of a Fréchet mean for a measure µ to the comparison between the distribution µ and the uniform measure λ on S 1 . The surprising fact is that λ appears as a benchmark to discriminate measures having a well defined Fréchet mean. The uniform measure λ is the 'worst' possible case as all points of the circle is a Fréchet mean, indeed the Fréchet functional (1.2) is constant and equals to π 3

.

In opposition to what have been done before we do not try to ensure convexity property on the Fréchet functional. Indeed, the definition of the Fréchet mean relies on the global optimization problem (1.3) which is, in general, non convex. The advantage of our approach is that we do not need to restrict the support or suppose restrictive conditions of symmetry on the density. As the geometry of flat manifold is simple, we can derive explicit form on the Fréchet functional and its derivative which can be hard to compute in non-flat manifolds such as n-dimensional spheres.

Organization of the paper

In Section 2, we introduce notations that will be used throughout the paper. In Section 3, we give explicit expressions for the Fréchet functional and its derivative and we discuss some properties of critical points of the Fréchet functional. Section 4 contains the main result with the necessary and sufficient condition of Theorem 4.1 for the existence of the Fréchet mean for a general measure. We also propose a new sufficient criterion P (α, ϕ) that ensures the well definiteness of the Fréchet mean. In Section 5, we study the convergence of the empirical Fréchet mean to the Fréchet mean, and describe an algorithm to compute the empirical Fréchet mean.

Notations

In what follows, 1 A denotes the indicator function of the set A ⊂ R and the notation b a f (t)dµ p 0 (t) stands for the Lebesgue integral

[a,b[ f (t)dµ p 0 (t) if a ≤ b and [b,a[ f (t)dµ p 0 (t) if b > a.

Normal coordinates

Given a base point p ∈ S 1 , there is a canonical chart called the exponential map defined from T p S 1 ≃ R, the tangent space of S 1 at p, to S 1 and denoted by

e p : R -→ S 1 θ -→ e p (θ) = R θ p, where R θ = cos(θ) -sin(θ) sin(θ) cos(θ)
This map is onto but not one to one as it is 2π-periodic. Such parametrizations are called normal coordinates systems centered at p and θ p 1 p 2 is nothing else but the coordinate of p 2 read in a normal coordinate system centered in p 1 . To simplify the notations, we will omit the exponent p 1 if no confusion is possible and we will write θ p 1 p 2 = θ p 2 . The cut locus of a point p 0 ∈ S 1 is denoted by p0 and is equal to the opposite point (in R 2 ) of p 0 , that is p0 = -p 0 . In a normal coordinates system centered at p ∈ S 1 , the coordinate of p0 is

θ p p0 = θ p p 0 -π if 0 ≤ θ p p 0 < π or θ p p0 = θ p p 0 + π if π ≤ θ p p 0 < 0.

Distance function and probability measures and the Fréchet functional

The arclength distance between two points p 1 , p 2 ∈ S 1 was defined in (1.1). Given normal coordinates θ p 1 , θ p 2 ∈ [-π, π[ of these points in the chart centered at an arbitrary point ,

d S 1 (p 1 , p 2 ) = d R/2πZ (θ p 1 , θ p 2 ) := min{|θ p 1 -θ p 2 + 2πk| , k ∈ Z}.
This means that the circle S 1 is locally isometric to the real line R. Unless specified, µ will denote a general probability measure on (S 1 , B(S 1 )) where B(S 1 ) is the Borel set of S 1 ⊂ R 2 . Given a point p 0 ∈ S 1 , µ p 0 is the image measure of µ through e -1 p 0 : S 1 -→ [-π, π[. This is a measure on R with a support in [-π, π[ which is defined by

µ p 0 (A) = µ • e p 0 (A ∩ [-π, π[), for all A ∈ B(R), (2.1)
where B(R) is the Borel set in R. The usual Euclidean mean/expectation and variance of µ p 0 are denoted m(µ

p 0 ) = R tdµ p 0 (t).
Finally, let us introduce for any p 0 ∈ S 1 , the map F µp 0 : [-π, π[-→ R given by

F µp 0 (θ) := F µ (e p 0 (θ)) = 1 2 π -π d 2 R/2πZ (t, θ)dµ p 0 (t) = 1 2 θ-π -π (t + 2π -θ) 2 dµ p 0 (t) + π θ-π (t -θ) 2 dµ p 0 (t), if 0 ≤ θ < π, θ+π -π (t -θ) 2 dµ p 0 (t) + π θ+π (t -2π -θ) 2 dµ p 0 (t), if -π ≤ θ < 0.
3 The Fréchet functional on the Circle ) and fix an arbitrary p 0 ∈ S 1 . Then, F µ : S 1 -→ R is differentiable in following sense :

1. Let p ∈ S 1 be a point with a cut locus of µ-measure 0, i.e µ({-p}) = 0. Then F µp 0 is continuously differentiable at θ p 0 p and we have

d dθ F µp 0 (θ p 0 p ) = θ p 0 p -2πµ p 0 ([-π, -π + θ p 0 p [) -m(µ p 0 ), if 0 ≤ θ p 0 p < π, θ p 0 p + 2πµ p 0 ([π + θ p 0 p , π[) -m(µ p 0 ), if -π ≤ θ p 0 p < 0. (3.1) 2. The function d dθ F µp 0 is left continuous on [-π, π[. Then we extend the definition of the derivative of F µ by setting for all θ ∈] -π, π[ d dθ F µp 0 (θ) := lim ε→0 - d dθ F µp 0 (θ + ε), (3.2
)

and d dθ F µp 0 (-π) = lim ε→0 -d dθ F µp 0 (π + ε). 3. Let p ∈ S 1 be a point with a cut locus of positive measure, i.e µ({-p}) > 0. Then, p is a cusp point of F µ in the sense that lim ε→0 -d dθ F µp 0 (θ p 0 p + ε) -lim ε→0 + d dθ F µp 0 (θ p 0 p + ε) = -µ({-p}).
Note that the left-continuity comes from our convention on the exponential map which is defined on

[-π, π[. If a measure µ is such that µ({p}) = 0 for all p ∈ S 1 then F µ is of class C 1 on [-π, π[.
Differentiability issues appear when the measure µ has atoms, see Figure 1.

Proof. For convenience we omit in this proof the superscript p 0 by writing θ p = θ p 0 p for all p ∈ S 1 . In the coordinate system centered at p 0 we have for all

θ p ∈ [-π, π[ F µp 0 (θ p ) = 1 2 R t 2 dµ p 0 (t) -θ p m(µ p 0 ) + 1 2 θ 2 p + 2π g + µp 0 (θ p )1 [0,π[ (θ p ) + g - µp 0 (θ p )1 [-π,0[ , (θ p ) (3.3) where g + µp 0 (θ) = -π+θ -π (π + t -θ)dµ p 0 (t) and g - µp 0 (θ) = π θ+π (π -t + θ)dµ p 0 (t).
Hence, to prove Proposition 3.1, we just have to study the derivative of g + µp 0 and g -

µp 0 . For all θ p ∈]0, π[ and ε ∈ R such that θ p + ε ∈]0, π[ we have, 1 ε g + µp 0 (θ p + ε) -g + µp 0 (θ p ) = 1 ε -π+θp+ε -π+θp (π + t -θ p )dµ p 0 (t) - -π+θp+ε -π dµ p 0 (t) (3.4)
The limit from the left in equation

(3.4) is lim ε→0 -1 ε (g + µp 0 (θ p + ε) -g + µp 0 (θ p )) = -µ p 0 ([-π, -π + θ p [) when 0 < θ p < π. The (left) derivative at θ p = 0 is given by lim ε→0 -1 ε g - µp 0 (θ p + ε) -g + µp 0 (θ p ) = 0. Similarly, if -π ≤ θ p < 0, we have lim ε→0 -1 ε (g - µp 0 (θ p + ε) -g - µp 0 (θ p )) = µ p 0 ([π + θ p , π[
) and statement 2 is proved.

To prove statement 1, suppose that the cut locus of p is of µ-measure 0. In this case, the limit from the left and from the right in equation (3.4) are equal as lim ε→0

1 ε -π+θp+ε -π+θp (π + t -θ p )dµ p 0 (t) = 0 since µ p 0 ({θ p -π}) = 0. Thence, formula (3.3) yields d dθ g + µp 0 (θ p ) = -µ p 0 ([-π, -π + θ p [), if θ p ∈ [0, π[ and d dθ g - µp 0 (θ p ) = µ p 0 ([π + θ p , π[), if θ p ∈ [-π, 0[. Finally, suppose that p ∈ S 1 has a cut locus of positive measure. If 0 ≤ θ p < π, it means that µ p 0 ({θ p -π}) > 0 and then, d dθ F µp 0 (θ p ) -lim ε→0 + d dθ F µp 0 (θ p + ε) = -lim ε→0 + µ p 0 ([-π + θ p , -π + θ p + ε[) = -µ p 0 ({-π +θ p }).
The case -π ≤ θ p < 0 is similar and the proof of statement 3 is completed.

Local minimum of the Fréchet functional

The critical points of the Fréchet functional are the points at which the derivative of F µ , in the sens of Proposition 3.1, is 0. As a immediate consequence of equation (3.1) we have

d dθ F µp 0 (θ p 0 p ) = -m(µ p ) (3.5)
Thus, the critical points are precisely the exponential barycenters (i.e. points p ∈ S 1 at which m(µ p ) = 0). This fact was already shown in [START_REF] Émery | Sur le barycentre d'une probabilité dans une variété[END_REF] or [START_REF] Le | Locating Fréchet means with application to shape spaces[END_REF] for general manifolds. Note that a critical point of the Fréchet functional is not in general an extremum (see an illustration Figure 1) but it is worth noticing that the minima of the Fréchet functional are regular points in the sens of the following result, Corollary 3.1. Let µ be a probability measure on S 1 . The cut locus of a (local or global) minimum of F µ is of µ-measure 0.

Proof. Let p ∈ S 1 be a point satisfying µ({-p}) > 0. For any p 0 ∈ S 1 , Statement 3 of Proposition 3.1 ensure that the derivative d dθ F µp 0 of the Frćhet functional has a negative jump at θ p 0 p . Hence, the signs of lim ε→0 -d dθ F µp 0 (θ p 0 pm + ε), lim ε→0 + d dθ F µp 0 (θ p 0 pm + ε) is either (+, +), (+, -) or (-, -). This means that θ p 0 p cannot be a minimum of F µp 0 since it would correspond to the case (-, +).

Remark 3.1. Note that assumptions of Corollary 1 in [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: basic tools for geometric measurements[END_REF] and Theorem 1 in [START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds. i[END_REF] contain a condition of the form µ({p}) to ensure (classical) differentiability of the Fréchet functional at its minimum. In the case of the circle, Corollary 3.1 shows us that the Fréchet functional is automatically differentiable at its minima.

4 Necessary and sufficient condition for the existence of the Fréchet mean

Main result

Theorem 4.1. Let µ be a probability measure and p * ∈ S 1 be a critical point of F µ . Then, the following propositions are equivalent, 1. p * is a well defined Fréchet mean of (S 1 , µ) .

2. For all 0 < θ < π

θ 0 λ([-π, -π + t[) -µ p * ([-π, -π + t[)dt > 0,
and for all -π ≤ θ < 0,

0 θ λ(]π + t, π[) -µ p * (]π + t, π[)dt > 0,
where λ is the uniform measure on [-π, π[ and µ p * is defined in (2.1).

Theorem 4.1 gives a necessary and sufficient condition for the existence of the Fréchet mean of a general measure µ on the circle S 1 . This condition is given in terms of comparison between the µ-measure and the uniform measure λ of balls centered at the cut locus of a global minimum. The important point is that the µ-measure of a (small) neighborhood of the cut locus of p * cannot be larger than the uniform measure of this neighborhood.

As µ is a probability measure, the functions

t -→ λ([-π, -π + t[) -µ p * ([-π, -π + t[) and t -→ λ(]π -t, π[) -µ p * (]π -t, π[) do not need to be always positive for t ∈ [-π, 0[ and t ∈ [0, π[ respectively.
An example where this quantity is always positive is when µ admits a density which is a decreasing function of the distance to a point p * , see [START_REF] Kaziska | The karcher mean of a class of symmetric distributions on the circle[END_REF]. In this case, the density is radially distributed around its mode p * which is, by Theorem 4.1, the Fréchet mean of µ. Many classical probability distributions used in circular data analysis follow this framework: von Mises distribution, wrapped normal distribution, geodesic normal distribution [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: basic tools for geometric measurements[END_REF] , etc... Another well-known example of distributions that admit a well defined Fréchet mean is distributions with support included in an hemisphere, see [START_REF] Samuel | Spherical averages and applications to spherical splines and interpolation[END_REF]. More precisely, suppose that there exists a point p ∈ S 1 with µ({p

∈ S 1 , d S 1 (p, p) ≤ π 2 }) = 1 and µ({p ∈ S 1 , d S 1 (p, p) < π 2 }) > 0. In this case, Statement 2 of Theorem 4.1 holds since one can show that the minimum p * of F µ is in {p ∈ S 1 , d S 1 (p, p) < π 2 } and that F µ p * (θ)-F µ p * (0) > 1 2 (θ+π-2θ p) 2 for θ ∈ [-π, θ p-π 2 [, F µ p * (θ)-F µ p * (0) = θ 2 2 for θ ∈ [θ p-π 2 , θ p+ π 2 [ and F µ p * (θ) -F µ p * (0) > 1 2 (θ -π -2θ p) 2 for all θ ∈ [θ p + π 2 , π[.
The case of equality corresponds to distributions with support in the boundary of the hemisphere, that is

µ p * = (1 -ε)δ θ p * p -π 2 + εδ θ p * p + π 2 with ε = 1 π θ p * p + 1 2
and in this case, there are two global argmins at 0 and 2θ p ± π, see Figure 2. (a) As already noticed, there exists at least one global argmin p * ∈ S 1 of F µ since F µ is a continuous function defined on the compact set S 1 . Moreover, Proposition 3.1 and Corollary 3.1 ensure that p * is a regular critical point of F µ , i.e. a zero of the derivative with a cut locus of µ-measure 0.

-3 -2 -1 0 1 2 3 1.4 1.6 1.8 2 2.2 2.4 (b) Figure 2: Let µ θ = (1 -ε)δ θ-π 2 + εδ θ+ π 2 with ε = θ π + 1 2 . (a) F µ θ with θ = 0. (b) F µ θ with θ = 3π
In the normal coordinate system centered at p * the functional G µ p * (θ) = F µ p * (θ) -F µ p * (0) has a particularly simple expression thanks to equation (3.5). Indeed, we have

G µ p * (θ) = 1 2 θ 2 + 2π1 [0,π[ (θ) θ-π -π (π + t -θ)dµ p * (t) + 2π1 [-π,0[ (θ) π θ+π (π -t + θ)dµ p * (t).
It is clear that Statement 1 is equivalent to the fact that the unique zero of G µ p * is 0. Thence, Theorem 5.1 will be an easy consequence of Lemma 4.1 below since we have λ([-π, -π + t[) = t 2π for all 0 ≤ t < π. Lemma 4.1. Let µ be a probability measure on S 1 and p * ∈ S 1 be an argmin of F µ . Then for any

θ ∈ [-π, π[ G µ p * (θ) = 2π        θ 0 t 2π -µ p * ([-π, -π + t[)dt, if 0 ≤ θ < π, 0 θ -t 2π -µ p * ([π + t, π[)dt, if -π ≤ θ < 0
Proof. The probability measure µ can be decomposed as follow,

µ = aµ d + (1 -a)µ δ , 0 ≤ a ≤ 1, (4.1)
where µ d is a probability measure such that µ d ({p}) = 0 for all p ∈ S 1 and µ δ = +∞ j=1 ω j δ p j where +∞ j=0 ω j = 1 and the p j 's are in S 1 . Hence, we consider the two cases separately : first, when the measure is non atomic, and then, when it is purely atomic. The general case follows immediately in view of equation (4.1).

First, assume that µ is an atomless measure of S 1 . Proposition 3.1 ensures that F µ is continuously differentiable everywhere and the real function F µ p * is of class C 1 on [-π, π[. Formula (3.1) and the fundamental theorem of calculus gives for all θ ∈ [-π, π[

G µ p * (θ) = θ 0 tdt -2π θ 0 µ p * ([-π, -π + t[)dt, if 0 ≤ θ < π, 0 θ µ p * ([π + t, π[)dt, if -π ≤ θ < 0. (4.2)
Consider now the case where µ is a purely atomic measure. First, we treat the case where the number of mass of Dirac in the sum is finite, i.e µ = n j=1 ω j δ p j , n ∈ N . Recall that F µ p * is a Lipschitz function on [-π, π[. Proposition 3.1 ensures that the derivative is piecewise continuous and formula (3.1) holds for all θ ∈ [-π, π[\{θ pj } n j=1 , i.e points that have a cut locus of µ-measure 0. Hence for all θ ∈ [-π, π[, equation (4.2) holds too.

To treat the case where µ = +∞ j=1 ω j δ x j we proceed by approximation. Let φ(n) = {j ∈ N | ω j ≥ 1 2 n } and remark that Card(φ(n)) < +∞ for all n ∈ N since

+∞ j=1 ω j = 1. Then if ν n p * = 1 c(n) j∈φ(n) ω j δ x j , where c(n) = j∈φ(n) ω j is a normalizing constant, we have for all θ ∈ [-π, π[, G ν n p * (θ) = θ 0 tdt -2π θ 0 ν n p * ([-π, -π + t[)dt, if 0 ≤ θ < π, 0 θ ν n p * ([π + t, π[)dt, if -π ≤ θ < 0.
The sequence (ν n p * ) n≥1 converges to µ in total variation. By the dominated convergence Theorem for all

θ p ∈ [-π, π[, G ν n p * (θ) converge as n → ∞ to (4.2).

The criterion P (α, ϕ)

Although the necessary and sufficient condition of Theorem 4.1 is a key step to understand the problem of non uniqueness of the Fréchet mean on S 1 , it is of little practice interest: we have to know a priori a critical point p * of the Fréchet functional. In this section, we derive sufficient conditions of existence with no restriction on the support of the probability measure and that are easily usable. Let us introduce the following definition : Definition 4.1. Let f : S 1 -→ R + be a probability density, p ∈ S 1 , α ∈]0, 1] and ϕ ∈]0, π[. We say that f satisfies the property where f p = f • e p : [-π, π[-→ R + . Moreover, we say that f ∈ P (α, ϕ) if there is a p ∈ S 1 such that f satisfies P (p, α, ϕ).

P (p, α, ϕ) if for all |θ| ≥ ϕ f p (θ) ≤ 1 -α 2π , (4.3)
The parameters α and ϕ control the concentration of µ around p. The idea is to control the mass lying in the complementary of the ball B(p, ϕ), see Figure 3. We have the following properties : Lemma 4.2. Let f : S 1 -→ R be a probability density on the circle. Then

1. P (p, α 1 , ϕ 1 ) =⇒ P (p, α 2 , ϕ 2 ) if α 1 ≥ α 2 and ϕ 1 ≤ ϕ 2 . 2. Let p 1 , p 2 ∈ S 1 and ϕ < π 2 . If d S 1 (p 1 , p 2 ) < π -ϕ then P (p 1 , α, ϕ) =⇒ P (p 2 , α, ϕ + d S 1 (p 1 , p 2 )). 3. If f satisfies P (p, α, ϕ) then |m(µ p )| ≤ ϕ + 1-α 4π (π -ϕ) 2 .
Proof. The first proposition is obvious in view of the Definition 4.1.

To prove the second claim suppose that 0 < θ p 1 p 2 = -θ p 2 p 1 ≤ πϕ (the other case is similar) and write

f p 2 (θ) = f p 1 (θ + θ p 1 p 2 )1 [-π,π-θ p 1 p 2 [ (θ) + f p 1 (θ + θ p 1 p 2 -2π)1 [π-θ p 1 p 2 ,π[ (θ).
In particular, it implies that π > πθ p 1 p 2 ≥ ϕ and since P (p 1 , α, ϕ) holds, we have

f p 2 (θ) ≤ 1-α 2π , if θ ≥ ϕ -θ p 1 p 2 or θ ≤ -ϕ -θ p 1 p 2 . This is equivalent to the fact that P (p 2 , α, min{| -ϕ -θ p 1 p 2 |, |ϕ -θ p 1 p 2 |}) = P (p 2 , α, ϕ + θ p 1
p 2 ) holds. The case ϕπ ≤ θ p 1 p 2 ≤ 0 is similar and we have P (p 2 , α, ϕθ p 1 p 2 ). Finally recall that |θ p 1 p 2 | = d S 1 (p 1 , p 2 ) and the property is proved. To show the last claim, we only need to consider the case where µ p has its support on [0, π[.

Indeed, µ p = ωµ - p + (1 -ω)µ + p where µ - p ([-π, 0[) = µ + p ([0, π[) = 1 and 0 ≤ ω ≤ 1. It yields that m(µ p ) = R td(ωµ - p + (1 -ω)µ + p ) = R + td(-ωµ - p + (1 -ω)µ + p ) ≤ R + tdµ + p . Then, if the density f p of µ p has its support in [0, π[ we have m(µ p ) ≤ ϕ 1 - π ϕ f p (t)dt + π ϕ tf p (t)dt ≤ ϕ + 1 -α 2π π ϕ (t -ϕ)dt,
which gives the result.

If the density f is sufficiently concentrated around a critical point p * of F µ then, this point is the Fréchet mean of µ. More precisely we have the following result, Proposition 4.1. Let µ be a probability measure with density f : S 1 -→ R + and p * ∈ S 1 be a critical point of F µ . If f satisfies P (p * , α, ϕ) with α ∈]0, 1] and 0 < ϕ < ϕ α = π √ α 1+ √ α then, µ admits a well defined Fréchet mean at p * . For all α ∈]0, 1], we have ϕ 0 = 0 ≤ ϕ α < π 2 = ϕ 1 . Note that if α = 1 then µ has its support included in the ball B(p * , π

2 ) = {p ∈ S, d S 1 (p * , p) ≤ π 2 } and when α < 1, the support of µ can be the entire circle S 1 .

Proof. Let G µ p * (θ) = F µ p * (θ) -F µ p * (0) for all θ ∈ [-π, π[. As the measure µ admits a density f , G µ p * is twice differentiable and equation (3.1) 2 and the second term is bounded from below by

implies that d 2 dθ 2 G µ p * (θ) = 1 -2πf(-π + θ), if 0 ≤ θ < π, and d 2 dθ 2 G µ p * (θ) = 1 -2πf(π + θ), if -π ≤ θ < 0. Since f ∈ P (p * , α, ϕ), the function G µ p * is convex on [-π + ϕ, π -ϕ]
G µ p * (θ) = G µ p * (π -ϕ) + θ π-ϕ t -2πµ p * ([-π, -π + t[)dt. Since f ∈ P (p * , α, ϕ) we have G µ p * (π -ϕ) ≥ α 2 (π -ϕ)
π π-ϕ t -2πν([-π, -π + t[)dt where ν = 1 2 (δ -ϕ + δ ϕ ). It yields for θ ∈ [π -ϕ, π[, G µ p * (θ) ≥ 1 2 (α(π -ϕ) 2 -ϕ 2 ). (4.4)
The right hand side of (4.4) is strictly positive if

ϕ < ϕ α = π √ α 1+ √ α . Similarly, the same condition implies G µ p * (θ) > 0 for θ ∈ [-π, -π + ϕ[.
We are now able to define a functional class of densities that admit a well defined Fréchet mean without restriction on the support of the measure. where α δ be the square of the root of (5

-6δ + δ 2 )X 3 + (1 -δ 2 )X 2 -(2δ + 1)X -1 that lies in ]0, 1] and ϕ α = π √ α 1+ √ α .
Then µ admits a well defined Fréchet means.

Firstly, remark that there is no need to know a critical point a priori. Secondly, the parameter δ controls the concentration of f via the inequality α δ < α and ϕ ≤ δϕ α . There is a tradeoff between α and the possible value of ϕ: the smaller α is (i.e the less f is concentrated) the smaller ϕ must be (i.e we need to control the value of the density on a bigger interval). In Tabular 1 we give examples of numerical values. Note that the column corresponding to δ = 0 is given as a reference only as the set P (α δ , δϕ α ) is empty for this values of δ. δ = 0

1 10 1 5 1 3 1 2
α δ ≤ 0.39 0.46 0.54 0.69 1 δϕ α δ ≥ 0 0.12 0.26 0.47

π 4
Table 1: Some values of α δ and δϕ α δ depending on δ ∈]0, 1 2 [.

Proof. We show that under the hypothesis of the Theorem 4.2, there is a critical point p * of F µ satisfying d S 1 (p, p * ) ≤ (1δ)ϕ α where p ∈ S 1 is a point satisfiying f ∈ P (p, α, ϕ). Thence, by Lemma 4.2, f belongs to P (p * , α, δϕ α + (1δ)ϕ α ) = P (p * , α, ϕ α ) and Proposition 4.1 will ensure that p * is the Fréchet mean of µ.

In the rest of the proof, we show that there is a p * ∈ S 1 such that d dθ F µp (θ p p * ) = 0 with d S 1 (p, p * ) ≤ (1δ)ϕ α . To this end, suppose that m(µ p ) ≥ 0 (the case m(µ p ) < 0 is similar). Equation (3.5) implies that d dθ F µp (0) = -m(µ p ) ≤ 0, and let us check that, under the hypothesis of Theorem 4.2, we have d dθ F µp ((1δ)ϕ α ) ≤ 0. Since d dθ F µp is a continuous function, the intermediate value Theorem will ensure the existence of a critical point p * such that

|θ p p * | ≤ (1 -δ)ϕ α . Equation (3.1) gives d dθ F µp ((1 -δ)ϕ α ) = (1 -δ)ϕ α -2πµ p ([-π, -π + (1 -δ)ϕ α [) -m(µ p ).
We have

-2πµ p ([-π, -π+(1-δ)ϕ α [) ≥ (α-1)(1-δ)ϕ α since f ∈ P (p, α, δϕ α ) and |-π + (1 -δ)ϕ α | ≥ δϕ α . Moreover, -m(µ) ≥ δϕ α -1-α 4π (π -δϕ α ) 2 by Lemma 4.2 Statement 3. It gives, d dθ F µp ((1 -δ)ϕ α ) ≥ π (5 -6δ + δ 2 )α √ α + (1 -δ 2 )α -(2δ + 1) √ α -1 1 + √ α
This quantity is positive as soon as 1 ≥ α > α δ , where √ α δ is the root of the polynomial X -→

(5

-6δ + δ 2 )X 3 + (1 -δ 2 )X 2 -(2δ + 1)X -1 that lies in ]0, 1].
It is easy to see that α 1 2 = 1 and numerical experiments show that the (increasing) function δ -→ α δ takes its value in ]0.39, 1[ for δ ∈]0, 1 2 [. [START_REF] Émery | Sur le barycentre d'une probabilité dans une variété[END_REF] Fréchet mean of an empirical measure

Existence

Let X 1 , . . . , X n be independent and identically distributed random variables with value in (S 1 , d S 1 ) and of probability distribution µ. The empirical measure is defined as usual by

µ n = 1 n n i=1
δ X i and we note p * n the empirical Fréchet mean defined as the unique argmin of

F µ n (p) = 1 2n n i=1 d 2 S 1 (p, X i ), p ∈ S 1 .
In [START_REF] Ziezold | On expected figures and a strong law of large numbers for random elements in quasi-metric spaces[END_REF] a strong law of large number is given for the empirical Fréchet mean in a semi metric space which is the case of (S 1 , d S 1 ). In particular, if p * n exists for each n ∈ N, the empirical Fréchet mean is a consistent estimator of the Fréchet mean. Indeed, the empirical Fréchet mean is well defined almost surely for a wide class of probability measures as the following fact from [START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds. i[END_REF] Remark 2.6 shows, Lemma 5.1. Let µ be a non atomic probability measure on the circle, i.e satisfying µ({p}) = 0 for all p ∈ S 1 . Then for all n ∈ N the empirical Fréchet mean exists almost surely.

Hence, the empirical Fréchet mean p * n of a probability measure µ can be computed even if µ does not possess a well defined Fréchet mean.

Consistency

If the Fréchet mean p * of µ is well defined, we study the rate of convergence of the empirical Fréchet mean p * n to p * . Proposition 5.1. Let µ be a measure with density f : S 1 -→ R that admits a well defined Fréchet mean p * . Then, there exists a strictly increasing function ρ :]0, π[-→]0, +∞[ such that for all p ∈ S 1 ,

p = p * F µ (p) -F µ (p * ) ≥ ρ(d S 1 (p, p * )).
If p * n denotes the empirical Fréchet mean, we have for all x > 0 We now focus on the proof of the concentration inequality (5.1) which is divided in two steps: first we show the uniform convergence in probability of F µ n p to F µp and then, we deduce the convergence of their argmins by using the lower bound given by the function ρ. Equation (3.1) implies that

P ρ(d S 1 (p * n , p * )) ≥ C(s) x n ≤ 2e -x . ( 5
2 sup θ∈[-π,π[ F µ n p (θ) -F µp (θ) ≤ 2π sup θ∈[-π,π[ d dθ F µ n p (θ) - d dθ F µp (θ) + 2 F µ n p (0) -F µp (0) ≤ 4π 2 m(µ p ) -m(µ n p ) + 2 m 2 (µ p ) -m 2 (µ n p ) (5.2) + 4π 2 sup θ∈[-π,π[ µ p ([-π, θ[) -µ n p ([-π, θ[) (5.3)
where m 2 (ν) = R t 2 dν(t), for a measure ν on R. The term (5.3) can be controlled in probability using the Dvoretzky-Kiefer-Wolfowitz inequality (see e.g [START_REF] Massart | The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality[END_REF]), and we have for all x > 0,

P 4π 2 sup θ∈[-π,π[ µ p ([-π, θ[) -µ n p ([-π, θ[) ≥ 4π 2 x n ≤ 2e -x .
For the terms of (5.2) which involve the first and second moment of µ p and µ n p , we use an Hoeffding type inequality which gives for all x > 0, P 4π 2 m(µ p )m(µ n p ) + 2 m 2 (µ p )m 2 (µ n p ) ≥ s(4π 2 + 2s) x n ≤ 2e -x , where s = max{|x -y| , x, y ∈ support µ} is the diameter of the support of µ. Combining these two concentration inequalities we have for all x > 0,

P 2 sup θ∈R F µp (θ) -F µ n p (θ) ≥ (4π 2 + 4π 2 s + 2s) x n ≤ 2e -x . (5.4)
We now use a classical inequality in M-estimation,

F µp (θ p * n ) -F µp (θ p * ) ≤ 2 sup θ∈[-π,π[ F µ n p (θ) - F µp (θ)
. By Lemma 5.2, there exists an increasing function ρ : R + -→ R + such that F µp (θ p * n ) -F µp (θ p * ) ≥ ρ(d S 1 (θ p * n , θ p * )). Plugging this in equation (5.4) we have,

P ρ(d S 1 (θ p * n , θ p * )) ≥ (4π 2 + 4π 2 s + 2s) x n ≤ 2e -x ,
and the proof of Proposition 5.1 is completed.

The function ρ that appears in the statement of Proposition 5.1 can be explicitly computed if the density f ∈ P (α, ϕ). The parameter α ∈]0, 1] can be interpreted as a measure of the convexity of F µp on the interval [-ϕ, ϕ]. For example, if α = 1 and ϕ

= ϕ α = π 2 , then µ has its support contained in [-π 2 , π 2 ] and F µp is quadratic on [-π 2 , π 2
]. Proposition 5.2. Let µ be a probability measure with density f : S 1 -→ R + satisfying the hypothesis of Theorem 4.2. Note p * the Fréchet mean of µ and for all x > 0 we have

P d S 1 (p * n , p * ) ≥ B(α, ϕ) x n 1 4 ≤ 2e -x , where B(α, ϕ) = C max π 2 γ(α,ϕ) , 2 α with γ(α, ϕ) = 1 2 (α(π -ϕ) 2 -ϕ 2 ) and C = 4π(2π 2 + π + 1).
Proof. This result is a direct consequence of Proposition 5.

[-π + ϕ α , π -ϕ α ], F µ p * (θ) -F µ p * (0) = 1 2 θ 2 d 2 dθ 2 F µ p * ( θ) ≥ α 2 θ 2 .
For all θ ∈ [-π, -ϕ α [∪]ϕ α , π[ we have by inequality (4.4),

F µ p * (θ) -F µ p * (0) ≥ 1 2 (α(π -ϕ α ) 2 -ϕ 2 α ) = γ(α, ϕ α ) > 0. Then, let ρ(t) = t 2 min{ γ(α,ϕα) π 2 , α 2 } ≥ t 2 min{ γ(α,ϕ) π 2 , α 2
} for any ϕ ≤ ϕ α .

Computation

Computation of the Fréchet mean of a general probability measure may not be an easy task as it is a global optimization problem. In practice the Fréchet functional is not a convex function and a gradient descent algorithm will only give a local minimum which depends on the initialization point. where n = 10 and µ n = n i=1 δ X i where the X i i.i.d of uniform law λ. The red points are the local minima computed with the described algorithm. (c) In blue: plot of F µ n where n = 4. In green: the derivative of F µ n . The critical points are given by the intersection between the green curve and the x-axis in black.

The Fréchet functional of an empirical measure is a continuous piecewise quadratic function as it can be written

F µ n p (θ) = 1 2n n i=1 θ -θ p X i +2π(1 [π+θ,π[ (θ p X i )1 [-π,0[ (θ)-1 [-π,θ-π[ (θ p X i )1 [0,π[ (θ)) 2 .
This formula together with Corollary 3.1 implies that the regular critical points (i.e point with no mass at their cut locus) of F µ are precisely the local minima of F µ . Moreover, the cumulative distribution function of µ n is, here, piecewise constant with jumps of size 1 n and we have

µ n p ([-π, t[) = 1 n n i=1 1 [-π,t[ (θ p X i ) = 1 n Card{θ p X i < t}.
Thence, the derivative of the empirical Fréchet functional is piecewise linear and to find the critical points amounts to solve n affine relations given by equation (3.1), see also Figure 4(c) for an illustration. Note, that in practice, there are less than n solutions, see Figure 4(a) and 4(b).

The following algorithm takes as input the values {X i } n i=1 and returns the Fréchet mean of µ n = 1 n n i=1 δ X i .

Initialization

Step: Choose an arbitrarily point p ∈ S 1 . Compute the coordinates {θ p X i } n i=1 and reorder them is increasing order. We denote τ - 0 = -π ≤ τ - 1 ≤ τ - 2 ≤ . . . ≤ τ - n 1 < 0 = τ - n 1 +1 the n 1 negative sorted terms and τ + 0 = π > τ + 1 ≥ . . .

≥ τ + n 2 ≥ 0 = τ + n 2 +1
the n 2 = nn 1 positive sorted terms.

Compute the mean m(µ n p ) = This algorithm can be extended to more general measures that the empirical one. The approach will be the same: find the critical points of the Fréchet functional with formula 3.1 and compare the values of the local minima. Unfortunately, there may be some computational issues as general cumulative distribution function will be not piecewise constant anymore.

Conclusion

It is not straightforward to extend criterion such as the one given in Theorem 4.1 to more general spaces, e.g. for the n dimensional sphere S n . Recall that the circle S 1 is a flat space in the sense that it is locally isometric to the Euclidean space R. Then, the only phenomenon that induces uniqueness issues of the Fréchet mean is the presence of a cut locus. The criterion presented in this note relies on an explicit formula for the gradient of the Fréchet mean. Curvature has an extra effect on the metric and makes difficult to derive exact computation on the Fréchet functional and its gradient. Moreover, it is not clear if the role played by the uniform measure as a benchmark in the well definiteness of the Fréchet mean in S 1 can be extended to n-spheres or non flat manifolds.

1

 1 Statistics for non-Euclidean data.

Figure 1 : Let µ = 1 6 δ p 1 + 2 3 δ p * + 1 6 δ p 2 with p 1 = R 2π 3 p * and p 2 = R -2π 3 p

 1323 Figure 1: Let µ = 1 6 δ p 1 + 2 3 δ p * + 1 6 δ p 2 with p 1 = R 2π 3 p * and p 2 = R -2π 3 p * . In blue: F µ p * . In green:

10 4. 2

 102 Proof of Theorem 4.1

Figure 3 :

 3 Figure 3: An example of distribution satisfying P (α, ϕ). Plot of the density f p in blue with the bounds P (p, 0.1, 2) in yellow and P (p, 0.5, 1.6) in red.

  and has a unique minimum at 0. Let us show that 0 is the only argmin of F µ p * on [-π, π[. If θ ∈ [πϕ, π[, we have thanks to Lemma 4.1

Theorem 4 . 2 .

 42 Let 0 < δ < 1 2 be a parameter of concentration and µ a probability measure with density f ∈ P (α, ϕ) (see Definition 4.1) with α δ ≤ α ≤ 1 and ϕ ≤ δϕ α

. 1 )Lemma 5 . 2 .

 152 where s = max{|x -y| , x, y ∈ support(µ)} and C(s) = (4π 2 + 4π 2 s + 2s) ≤ 4π(2π 2 + π + 1). The function ρ in the statement of the Proposition 5.1 determines the rate of convergence of p * n to p * . Proof. The first claim about the lower bound ρ is a direct consequence of the following Lemma: Let f : [-π, π[-→ R + be a continuous function on [-π, π[ (see Section 3). If f vanishes at a unique point θ 0 ∈ [-π, π[, then there exists a strictly increasing function ρ:]0, π[-→]0, +∞[ such that for all θ ∈ [-π, π[\{0}, f (θ) ≥ ρ(d S 1 (θ 0 , θ)).Proof. As f is the restriction on [-π, π[ of a continuous periodic function we can assume, without loss of generality, that θ 0 = 0. Then, define for all θ ∈]0, dt, where g(t) = min f (-t), f (t), min |τ |>t f (τ ) , for all t ∈ [-π, π[.

Figure 4 :

 4 Figure 4: (a) and (b) Plots of F µ n where n = 10 and µ n = ni=1 δ X i where the X i i.i.d of uniform law λ. The red points are the local minima computed with the described algorithm. (c) In blue: plot of F µ n where n = 4. In green: the derivative of F µ n . The critical points are given by the intersection between the green curve and the x-axis in black.

  To guaranty the injectivity, we choose to restrict the domain of definition R of e p to [-π, π[. Thus, for all p 1 , p 2 ∈ S 1 there is a now uniqueθ p 1 p 2 ∈ [-π, π[ satisfying e p 1 (θ p 1 p 2 ) = p 2 and

	e p : [-π, π[-→ S 1	and	e -1 p : S 1 -→ [-π, π[,	for all p ∈ S 1 .

  3.1 The derivative of the Fréchet functionalA function f : [-π, π[-→ R is said left continuous on [-π, π[ if it is left continuous everywhere on ]π, π[ and with lim ε→0 -f (π + ε) = f (-π).Similarly, f is said to be continuous on [-π, π[ if it is left and right continuous on [-π, π[. We provide here an explicit expression of the derivative of F µ , Proposition 3.1. Let µ be a probability measure on (S 1 , d S 1

  1 and we only have to find a strictly increasing function ρ :[0, π] -→ R + satisfying for all θ ∈ [-π, π[, F µ (p) -F µ (p * ) ≥ ρ(d S 1 (p, p * )).As µ p * admits a density f p * , the Fréchet functional F µ p * is twice differentiable and f ∈ P (p * , α, ϕ α ), see proof of Theorem 4.2. Thus, for all θ ∈ [-π + ϕ α , ϕ απ], a second order Taylor expansion of F µ p * at 0 ensures that for some θ ∈

  1 n (τ - 1 + . . . + τ - n 1 + τ + 1 + . . . + τ + n 2 ) and initialize θ p p * to 0, says. # The first step compares all the local minima in [0, π[ Step 1: For i from 0 to n 1 do # θ p p * ,new is the candidate to be a critical point between τ - i and τ - ,new is a critical point and then test its value. If better, keep it. if τ - i + π ≤ θ p p * ,new ≤ τ - i+1 + π and F µ n p (θ p p * ,new ) ≤ F µ n p (θ p p * ) then θ p p * := θ p p * ,new end if end for. # The Step 2 is the same as Step 1 but for local minima in [-π, 0[ Step 2: For i = 0 to n 2 do Let θ p p * ,new = -2π i n + m(µ n p ) if τ + i+1π ≤ θ p p * ,new ≤ τ + iπ and F µ n p (θ p p * ,new ) ≤ F µ n p (θ p p * ) then θ p p * := θ p p * ,new end if end for. # The value of θ p p * is the best argmin Output: Return p * = e p (θ p p

	Let θ p p * ,new = 2π i n + m(µ n p ) # verify if θ p p	i+1

* * ).

Acknowledgements. The author would like to thank Dominique Bakry for help and encouragement during the writing of this paper and Jérémie Bigot for his careful reading of the manuscript. The author acknowledges the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-JCJC-SIMI1 DEMOS.