Image, Data, Roland Badeau
email: roland.badeau@telecom-paristech.fr

Stochastic reverberation model for uniform and non-diffuse acoustic fields Modèle stochastique de réverbération pour les champs acoustiques uniformes et non diffus

Keywords: Données, Signal Groupe Reverberation, Diffusion, Room impulse response, Stochastic models Réverbération, Diffusion, Réponse impulsionnelle de salle, Modèles stochastiques

In a recent research report, we introduced a general stochastic reverberation model that aims to represent the statistical properties of reverberation in a broad variety of acoustic environments. A simplified version of this model, dedicated to the particular case of diffuse (i.e. uniform and isotropic) acoustic fields, omnidirectional sources and microphones, and constant attenuation w.r.t frequency, has been investigated both mathematically and experimentally in a recent research paper. We showed that this model provides a common mathematical framework that unifies several well-known results regarding the statistical properties of reverberation in the space, time and frequency domains.

In this research report, we aim to extend this mathematical analysis to uniform and non-diffuse acoustic fields, and directive sources and microphones. We show that the predictions of the general stochastic model experimentally match the observations, based on both synthetic and real room impulse responses, measured in various acoustic environments.

I. INTRODUCTION

In [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF], we introduced a common mathematical framework for stochastic reverberation models, that aimed to unify several well-known results regarding the statistical properties of reverberation, in the spatial, temporal and spectral domains [START_REF] Cook | Measurement of correlation coefficients in reverberant sound fields[END_REF]- [START_REF]Playing billiards in the concert hall: The mathematical foundations of geometrical room acoustics[END_REF]. This framework was based on the source image principle [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF], [START_REF] Kuttruff | Room Acoustics, Fifth Edition[END_REF], which represents the sound wave reflected by a flat surface as if it was emitted by a so-called source image. In [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF], the positions of the source images were modeled as random, and uniformly distributed according to a Poisson point process.

However, the stochastic reverberation model introduced in [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] was limited to diffuse (i.e. uniform and isotropic 1 ) acoustic fields, omnidirectional sources and microphones, and constant attenuation w.r.t. frequency. In [START_REF] Badeau | General stochastic reverberation model[END_REF], we proposed several extensions of this model, that aim to represent reverberation more realistically, by considering anisotropic and nonuniform acoustic fields, directive sources and microphones, and frequency-varying attenuation coefficients.

In this research report, we aim to investigate the statistical properties of this generalized model. In order to keep the mathematical analysis as simple as possible, we chose to restrict the study to uniform acoustic fields and constant attenuation w.r.t. frequency as in [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF], but we address the generalization to anisotropic acoustic fields and directive sources and microphones as introduced in [START_REF] Badeau | General stochastic reverberation model[END_REF]. We will first provide a mathematical analysis of the statistical properties of the model, regarding its first and second order moments in the spatial, temporal, spectral and time-frequency domains, and its asymptotic Gaussianity. Then we will show that the predicted statistical properties experimentally match the observations, based on both synthetic and real room impulse responses (RIRs), measured in various acoustic environments.

This research report is organized as follows: the uniform stochastic reverberation model will be formally defined in Section II, then its statistical properties will be analyzed mathematically in Section III and experimentally in Section IV. Conclusions will be drawn in Section V, and the proofs of the mathematical results will be presented in Appendices A to G. Throughout the report, we will use the following mathematical notation:

• N: set of whole numbers; • R, C: sets of real and complex numbers, respectively; • R + : set of nonnegative real numbers;

• ı = √ -1: imaginary unit;

• : equal by definition to;

• c =: equal up to an additive constant to; • ∝: proportional to; • x (bold font), z (regular): vector and scalar, respectively;

• z: complex conjugate of z ∈ C; • |.|: absolute value of a scalar or a vector (entrywise); • x m : m-th entry of vector x; • x : transpose of vector x; • S 2 : unit sphere in R 3 (S 2 = {x ∈ R 3 ; x 2 = 1}); • {a, . . . , b}: set including all elements listed from a to b; • L ∞ (V ), where V is a Borel set: Lebesgue space of essentially bounded functions f of support V (i.e. such that f ∞ = ess sup V |f | < +∞); • L p (V ), where V is a Borel set and p ∈ N\{0}: Lebesgue space of measurable functions f of support V , such that f p = ( V |f (x)| p dx) 1 p < +∞;

• . 2 : Euclidean/Hermitian norm of a vector or a function; • J q (.): Jacobian matrix of the multivariate function q(.); • O(.): asymptotically bounded by;

• o(.): asymptotically dominated by;

• ∼: has distribution or is asymptotically equivalent to;

• δ: Dirac delta function;

• E[X]: expected value of a complex random variable X; • φ X (θ) = E[e ıRe (θX) ]: characteristic function of X;

• Covariance of two complex random variables X and Y :

cov[X, Y ] = E[(X -E[X])(Y -E[Y ])];
• var[X] = cov[X, X]: variance of a random variable X;

• Correlation of two complex random variables X and Y :

corr[X, Y ] = cov[X, Y ] var[X]

var[Y ] ;

• P(λ): Poisson distribution of parameter λ > 0:

N ∼ P(λ) ⇔ P (N =n) = e -λ λ n n! ⇔ φ N (θ) = e λ(e ıθ -1) ;

• sinc(x) = sin(x)

x : cardinal sine function;

• 1 A : indicator function of a set A (1 A (x) is 1 if x ∈ A
or 0 if x / ∈ A); • ψ(t) = ψ(-t): conjugate and time-reverse of ψ : R → C; • Convolution of two functions ψ 1 and ψ 2 : R → C:

(ψ 1 * ψ 2 )(t) = u∈R ψ 1 (u)ψ 2 (t -u)du;

• Fourier transform of a function ψ : R → C: ψ(f ) = t∈R ψ(t)e -2ıπf t dt (f ∈ R);

• ⊗: outer product ((ψ 1 ⊗ ψ 2 )(t 1 , t 2 ) = ψ 1 (t 1 )ψ 2 (t 2 ));

• Wigner distribution (a.k.a. Wigner-Ville distribution) of a function γ : R 2 → R:

W γ (t, f ) = τ ∈R γ(t + τ 2 , t -τ
2 )e -2ıπf τ dτ. [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF] Note that the Wigner distribution satisfies several important properties [START_REF] Cohen | Time-frequency distributions-a review[END_REF]:

• even symmetry: ∀t, f ∈ R, W γ (t, -f ) = W γ (t, f );

• real property: if γ is symmetric (∀t 1 , t 2 ∈ R, γ(t 1 , t 2 ) = γ(t 2 , t 1 )), then ∀t, f ∈ R, W γ (t, f ) ∈ R; • projection property: for any functions ψ 1 , ψ 2 : R → R,

t∈R W ψ1⊗ψ2 (t, f )dt = ψ 1 (f ) ψ 2 (f ), (2)
f ∈R W γ (t, f )df = γ(t, t); (3)
• convolution property: if γ(t 1 , t 2 ) = (γ

1 t1 * t2 * γ 2 )(t 1 , t 2 ), W γ (t, f ) = (W γ1 t * W γ2 )(t, f ), (4)
where t * denotes convolution over variable t; • temporal support property: if the support of function γ : R 2 → R is compact, then the temporal support of W γ (t, f ) is also compact. Finally, in order to model the spatial distribution of the source images, we will use the concept of Poisson random measure with independent increments (also refered to as Poisson point process), as we did in [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF]: given a Borel set V ⊂ R 3 of finite volume |V |, we assumed that the number N (V ) of source images contained in V follows a Poisson distribution of rate parameter λ|V |: N (V ) ∼ P(λ|V |) with λ > 0. Formally, given a non-negative, locally integrable function Λ(x) on R p , the Poisson random increment dN (x) ∼ P(Λ(x)dx) corresponds to an infinitesimal volume |V | = dx. Then for any Borel set V ⊂ R p of finite Lebesgue measure, the number N (V ) = V dN (x) of points contained in V follows a Poisson distribution of rate parameter V Λ(x)dx: N (V ) ∼ P( V Λ(x)dx), and for any disjoint Borel sets V 1 and V 2 , N (V 1 ) and N (V 2 ) are independent. Note that in the stochastic reverberation model proposed in [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF], we considered a spatially uniform distribution of the source images in the 3D-space, so that p = 3 and Λ(x) = λ > 0 is constant.

II. UNIFORM STOCHASTIC REVERBERATION MODEL

In this research report, we consider the following stochastic reverberation model (see [START_REF] Badeau | General stochastic reverberation model[END_REF] for a description of the basic geometric principles that underlie this model):

Definition 1 (Uniform stochastic reverberation model). For any sensor i ∈ {1 . . . I}, time t ∈ R, and frequency f ∈ R, let

h i (t) = x∈R 3 y∈R M V i (x, y; y-q(x-x i )) g i x-xi x-xi 2 , . t * s Θ(x,y) x-xi x-xi 2 ,t- x-x i 2 c
e -y α c

x-xi 2 dN (x, y),

(5)

h i (f ) = x∈R 3 y∈R M V i (x, y; y-q(x-x i )) g i x-xi x-xi 2 , f s Θ(x,y) x-xi x-xi 2
,f e -y α+2ıπf x-x i 2 c

x-xi 2 dN (x, y),

where

• h i (t) ∈ R (resp. h i (f ) ∈ C)
is the room impulse response (resp. room frequency response) at sensor i; • g i (u, t) ∈ R (resp. g i (u, f ) ∈ C) is the response of sensor i at direction u ∈ S2 (taking into account both its directivity and orientation) and time t (resp. frequency f ); • s (u, t) ∈ R (resp. s (u, f ) ∈ C) is the response of the source at direction u ∈ S 2 (taking into account its directivity) and time t (resp. frequency f ); • vector x i ∈ R 3 (in meters) is the position of sensor i;

• vector x ∈ R 3 (in meters) represents the possible positions of the source images; • c > 0 is the speed of sound (in meters.hertz);

• the term e -2ıπf x-x i 2 c

x-xi 2 corresponds to the propagation of a monochromatic spherical wave from x to x i ;

• α ∈ R M
+ is a vector of attenuation coefficients (in hertz); • the term e -y α c corresponds to the total attenuation of the sound wave during its propagation from x to x i ; • r ∈ R 3 → q(r) ∈ R M + (in meters) is a 1-homogeneous function (i.e. ∀r ∈ R 3 , ∀ν ∈ R + , q(νr) = νq(r));

• the coordinates (x, y) ∈ R 3 × R M (in meters) are distributed according to a uniform Poisson point process dN (x, y) ∼ P(λdxdy) with λ > 0 (in meters -(3+M ) ); • Θ(x, y) ∈ SO(3) is a random rotation matrix that represents the orientation of the source image of coordinates (x, y). Its distribution is i.i.d. w.r.t. (x, y) and not necessarily uniform on the rotation group SO(3); • V i (x, y; z i ) ∈ {0, 1} is a Boolean that indicates whether the source image of coordinates (x, y) is visible from sensor i or not. Formally, V i (x, y; z i ) is a random field on R 3 × R M , whose probability distribution is parameterized by the vector z i ∈ R M . The joint distribution for all sensors i of the random vector V (x, y; z 1 . . . z I ) [V 1 (x, y; z 1 ) . . . V I (x, y; z I )] i∈{1...I} is i.i.d. w.r.t. (x, y), and it is denoted p(b 1 . . . b I ; z 1 . . . z I ) ∈ [0, 1] where b i ∈ {0, 1} and z i ∈ R M . The marginals for every sensor i are denoted p(b; z) ∈ [0, 1] (they are such that the closer z is to zero, the higher p(1; z)), and the marginals for every pair of sensors (i, j) are denoted p(b i , b j ; z i , z j ) ∈ [0, 1]; • the three random fields dN , Θ and V on R 3 × R M are independent. Then let the attenuation function be defined as ∀r ∈ R 3 , α(r) = q(r) α ≥ 0.

(

) 7
Based on [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF], we define the following prime notation for functions:

• for any function ψ(u, t) defined on S 2 ×R, let ψ (u, t) = ψ(u, t) e α(u)t ; • if function α(.) is constant on S 2 (i.e. ∀u ∈ S 2 , α(u) = α inf ), then for any function ψ(t) defined on R, let ψ (t) = ψ(t) e α inf t . In addition, we assume that:

• all functions defined above are not identically zero;

• r → q(r) is four times continuously differentiable; -∀i ∈ {1 . . . I}, g i (u, 0) = 0 and m s (u, 0) = 0;

• α inf inf u∈S 2 α(u) > 0; • ∀i ∈ {1 . . . I}, function g i (u,
-∀i ∈ {1 . . . I}, g i (u, 0) = 0 and

∂ g i ∂f (u, 0) = 0; -m s (u, 0) = 0 and ∂m s ∂f (u, 0) = 0; where ∀u ∈ S 2 , ∀f ∈ R, ∀(x, y) ∈ R 3 × R M , m s (u, f ) = E s (Θ(x, y)u, f ) ; (8)
• function z → p(1; z) is continuously differentiable and it is not constant;

• the support of function z → p(1; z) is left-bounded, i.e. ∃z inf ∈ R M such that ∀m ∈ {1 . . . M }, z inf m < 0, and ∀z ∈ R M , if ∃m ∈ {1 . . . M } such that z m < z inf m , then p(1; z) = 0; • function (z 1 , z 2 ) → p(1, 1; z 1 , z 2 ) is continuous and
differentiable almost everywhere in R M × R M , and all its partial derivatives belong to

L ∞ (R M × R M ); • (z 1 , z 2 ) → p(1, 1; z 1 , z 2 ) is such that ∀z ∈ R M , p(1, 1; z, z) = p(1; z). (9)
Note that ( 9) is derived from geometrical considerations: when x j → x i , a source image of coordinates (x, y) is visible from both sensors i and j if and only if it is visible from sensor i. This can be expressed as

∀x ∈ R 3 , ∀y ∈ R M , lim xj →xi p(1, 1; y -q(x-x i ), y -q(x-x j )) = p(1; y -q(x-x i )),
which implies [START_REF] Polack | Modifying chambers to play billiards: the foundations of reverberation theory[END_REF]. Moreover, because of (9), function (z 1 , z 2 ) → p(1, 1; z 1 , z 2 ) cannot be smoother than what we assumed in Definition 1 (cf. Lemma 5 in Appendix D) 2 .

Also note that since functions g i (u, t) and s (u, t) are causal ∀u ∈ S 2 , the RIR h i (t) defined in ( 5) is also causal. Besides, all functions g i (u, t), g i (u, t), s (u, t) and s (u, t) are continuous w.r.t. t ∈ R and twice continuously differentiable w.r.t. u ∈ S 2 , and since they all have finite temporal support, all functions g i (u, f ), g i (u, f ), s (u, f ) and s (u, f ) are smooth w.r.t. f ∈ R and twice continuously differentiable w.r.t. u ∈ S 2 . Moreover, since r → q(r) is four times continuously differentiable, function r → α(r) is also four times continuously differentiable.

The stochastic model in Definition 1 is a particular case of the general stochastic reverberation model presented in [START_REF] Badeau | General stochastic reverberation model[END_REF], from which we introduced two simplifications 3 :

• the attenuation coefficients in vector α are constant and do not depend on frequency f ; • the acoustic field is uniform: function r → q(r) does not depend on sensor i. In this report, we will address several particular cases of interest:

Definition 2 (Diffuse acoustic field). Considering the uniform stochastic reverberation model in Definition 1, the acoustic field is diffuse when the two following conditions hold:

• ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant (dimensionless) vector, therefore function α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α; • the random rotation matrices Θ(x, y) are uniformly distributed on SO(3).

Definition 3 (Omnidirectional sensor).

Considering the stochastic reverberation model in Definition 1, sensor i is omnidirectional when

g i (u, t) = g i (t) ∀t ∈ R and g i (u, f ) = g i (f ) ∀f ∈ R do not depend on u ∈ S 2 .
Definition 4 (Omnidirectional source). Considering the stochastic reverberation model in Definition 1, the source is

omnidirectional when s (u, t) = s (t) ∀t ∈ R and s (u, f ) = s (f ) ∀f ∈ R do no not depend on u ∈ S 2 .
Note that in the case of an omnidirectional source, the random rotation matrices Θ(x, y) disappear from Definition 1.

Finally, in order to characterize the correlations between sensors in the time-frequency domain, we will use two kinds of indicators that will be illustrated in the experiments presented in Section IV: Definition 5 (Time-frequency correlation). Considering the stochastic reverberation model in Definition 1, the timefrequency correlation between two sensors i, j ∈ {1 . . . I} is defined ∀t, f ∈ R as

ρ i,j (t, f, x j -x i ) W γi,j (t, f ) W γi,i (t, f )W γj,j (t, f ) (10)
3 Moreover, compared to [START_REF] Badeau | General stochastic reverberation model[END_REF], we simplified the mathematical expressions (without loss of generality), by removing the rotation matrix

Θ i in g i Θ i x-x i x-x i 2 , t and g i Θ i x-x i x-x i 2
, f , and by changing the sign of the rotation matrices Θ(x, y) in s Θ(x, y)

x i -x x i -x 2 , t and s Θ(x, y) x i -x x i -x 2 , f .
whenever the denominator is positive, where W denotes the Wigner distribution defined in (1), and

γ i,j (t 1 , t 2 ) = cov[h i (t 1 ), h j (t 2 )]. (11)
In [START_REF]Playing billiards in the concert hall: The mathematical foundations of geometrical room acoustics[END_REF], we introduced the redundant notation ρ i,j (t, f, x jx i ) (which makes a double use of indexes i and j) in order to insist on the fact that ρ i,j will depend on x i and x j , only via their difference x j -x i , since the acoustic field is uniform (i.e. invariant under any translation). Also note that, since function f → W γi,j (t, f ) is even symmetric, function f → ρ i,j (t, f, x j -x i ) is also even symmetric.

If the time-frequency correlation converges when t → +∞, we also define the following indicator: Definition 6 (Asymptotic correlation function). Considering the stochastic reverberation model in Definition 1, the asymptotic correlation function between two sensors i, j ∈ {1 . . . I} is defined ∀τ ∈ R as

σ i,j (τ, x j -x i ) = f ∈R lim t→+∞ ρ i,j (t, f, x j -x i )e +2ıπf τ df, (12)
whenever the limit exists and the integral converges, where ρ i,j was defined in [START_REF]Playing billiards in the concert hall: The mathematical foundations of geometrical room acoustics[END_REF].

Note that, since function

f → ρ i,j (t, f, x j -x i ) is even symmetric, σ i,j (τ, x j -x i ) ∈ R ∀τ ∈ R.

III. MATHEMATICAL ANALYSIS

In this section, we analyze the statistical properties of the stochastic model introduced in Definition 1, in terms of first order (Section III-A) and second order moments (Section III-B). Then we will be able to relax some restricting assumptions regarding the acoustic field, the source and the microphones, by focusing on late reverberation only, i.e. by assuming that t → +∞. Actually, the mathematical analysis in Section III, as well as the observations in Section IV based on both synthetic and real RIRs, will show that two different asymptotic states can be distinguished: at first, the variations of the attenuation α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] over the direction u ∈ S 2 can be neglected; this period of time will be called early asymptotic state (Section III-C). Then after a while, the directions u that are less attenuated (i.e. such that α(u) > 0 is the lowest) start to dominate the other directions; this period of time will be called late asymptotic state (Section III-D). In both asymptotic states, we will prove that the RIR h i (t) is asymptotically normally distributed.

A. First order moments

Before analyzing the first order moments, we need to define some coefficients that are related to the discrete probability distribution p introduced in Definition 1, and that will play an important role throughout this research report: Lemma 1. Let p(b, z) (where b ∈ {0, 1} and z ∈ R M ) be the probability distribution introduced in Definition 1. Let ∀n ∈ N\{0},

β n = z∈R M p(1; z)e -n z α c dz > 0. ( 13
)
Then we have the following majoration: ∀n ∈ N\{0},

β n = O   e n |z inf | α c n M   . ( 14
)
Lemma 1 is proved in Appendix D. We are now ready to investigate the first order moments of the stochastic reverberation model: Proposition 1 (First order moments). Given the stochastic reverberation model in Definition 1, the room response has the following first order moments: for any sensor i ∈ {1 . . . I},

• Temporal domain: ∀t ≥ T T g + T s , E[h i (t)] = 0. • Spectral domain: let ∀u ∈ S 2 , ∀f ∈ R, m s (u, f ) = E [ s (Θ(x, y)u, f )] (15)
∀(x, y) ∈ R 3 × R M . Then function m s is twice continu- ously differentiable w.r.t. u ∈ S 2 and smooth w.r.t. f ∈ R. Moreover, ∀f ∈ R, E[ h i (f )] = λc 2 β 1 u∈S 2 g i (u, f ) m s (u, f ) ( α (u) + 2ıπf ) 2 du (16)
where β 1 is defined in [START_REF] Badeau | General stochastic reverberation model[END_REF] for n = 1, and function

f → E[ h i (f )] is smooth.
Proposition 1 is proved in Appendix E. Note that this proposition generalizes the results already established in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF]. In particular, h i (t) is centered for t ≥ T (the fact that it is not centered for t ∈ [0, T ] explains why the expected value of the frequency response E[ h i (f )] in ( 16) is not zero).

B. Second order moments

Before analyzing the second order moments, we need to further investigate the properties of the discrete probability distribution p introduced in Definition 1.

Lemma 2. Let p(b i , b j , z i , z j ) (where b i , b j ∈ {0, 1} and z i , z j ∈ R M ) be the joint probability distribution introduced in Definition 1. Let ∀e ∈ R M , β(e) = z∈R M p(1, 1; z -e 2 , z + e 2 )e -2z α c dz ≥ 0. ( 17
)
Then function e → β(e) is even. Moreover, it is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong to

L ∞ (R M ). At e = 0, β(0) = β 2 , where β 2 is defined in (13) for n = 2. Besides, ∀e ∈ R M , β(e) ≤ β 2 e -|e α| c . ( 18
)
In particular, function e → β(e) reaches its maximum at e = 0, and it is not differentiable at e = 0.

Lemma 2 is proved in Appendix D. Note that, in connection with the remark in Footnote 2, the fact that function β(.) is not twice continuously differentiable is the reason for the slower speed of convergence O( 1 √ t ) in late asymptotic state (see the discussion in Section III-E).

In the following lemma, we introduce a function that will play a prominent role to characterize the second order moments of the RIR in the spatial, temporal, spectral and timefrequency domains. Lemma 3. Given the stochastic reverberation model in Definition 1, let β(.) be the function introduced in Lemma 2 in [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept sine technique[END_REF], let β 2 be defined in [START_REF] Badeau | General stochastic reverberation model[END_REF] for n = 2, and let D ≥ 0 denote the distance between two sensors. We also assume that ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant vector, therefore function α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF]

is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α. If D > 0, let b(τ, D) = c 2D β (qcτ ) 1 [-D c ,+ D c ] (τ ) ≥ 0 ∀τ ∈ R, (19) otherwise if D = 0, let b(τ, 0) = β 2 δ(τ ) and b(f, 0) = β 2 ∀f ∈ R. Then function τ → b(τ, D) is even. Moreover,
(τ, D) ∈ L ∞ (] -D c , + D c [). Function f → b(f, D
) is smooth and real-valued. Moreover, we have the two majorations:

∀τ ∈ R, b(τ, D) ≤ cβ2 2D e -α inf |τ | 1 [-D c ,+ D c ] (τ ), (20)
and ∀f ∈ R, | b(f, D)| ≤ β 2 . (21)
In particular, function τ → b(τ, D) reaches its maximum at τ = 0, and it is not differentiable at τ = 0. Besides, when D → 0, we get

b(τ, D) = cβ2 2D + O(1) 1 [-D c ,+ D c ] (τ ) (22)
and b(f, D) = β 2 sinc 2πf D c + O(D). ( 23
)
Lemma 3 is proved in Appendix F. We are now ready to investigate the second order moments of the stochastic reverberation model in the case of diffuse acoustic fields: Proposition 2 (Second order moments in diffuse acoustic fields). Considering the stochastic reverberation model in Definition 1, suppose that the acoustic field is diffuse (cf. Definition 2) and that the source and microphones are omnidirectional (cf. Definitions 3 and 4). For any sensors i, j ∈ {1 . . . I}, let D = x i -x j 2 denote the distance between the two sensors, and let γ i,j (t 1 , t 2 ) be the function defined in [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF]. Finally, let b(τ, D) be the function introduced in Lemma 3. Then the room response has the following second order moments:

• Temporal domain: ∀t 1 + t 2 ≥ 2T + D c with T T g + T s , γ i,j (t 1 , t 2 ) = 4πλc e -α inf (t1+t2) b(., D) * g i * g j * s * s (t 1 -t 2 ), ( 24
) and function (t 1 , t 2 ) → γ i,j (t 1 , t 2 ) is continuous. In particular, if i = j, ∀t ≥ T , var[h i (t)] = 4πλc β 2 e -2 α inf t g i * g i * s * s (0) (25)
and

∀t 1 + t 2 ≥ 2T + D c , the temporal correlation corr[h i (t 1 ), h j (t 2 )]= (b(.,D) * g i * g j * s * s )(t1-t2) β2 ((g i * g i * s * s )(0))((g j * g j * s * s )(0)) (26) only depends on t 1 -t 2 . • Spectral domain: ∀f 1 , f 2 ∈ R, cov[ h i (f 1 ), h j (f 2 )] = 2πλc b( f 1 +f 2 2 ,D) gi(f1) gj (f2) s(f1) s(f2)e -( α inf +ıπ(f 1 -f 2 ))D c α inf +ıπ(f1-f2) , (27) and function (f 1 , f 2 ) → cov[ h i (f 1 ), h j (f 2 )] is smooth. In particular, if i = j, ∀f ∈ R, var[ h i (f )] = 2πλcβ2 | gi(f )| 2 | s(f )| 2 α inf (28) and ∀f 1 , f 2 ∈ R, corr[ h i (f 1 ), h j (f 2 )]= e ı(∠ gi(f1)-∠ gj (f2)+∠ s(f1)-∠ s(f2)) b( f 1 +f 2 2 ,D)e -( α inf +ıπ(f 1 -f 2 ))D c β2 1+ ıπ(f 1 -f 2 ) α inf . ( 29
) • Time-frequency domain: ∀t ≥ T + D 2c , ∀f ∈ R, W γi,j (t, f ) = e -2 α inf t B i,j (f, D) (30)
where

B i,j (f, D) = 4πλc b(f, D) g i (f ) g j (f )| s (f )| 2 , ( 31
)
and function f → B i,j (f, D) is smooth. In particular, if i = j, W γi,i (t, f ) = 4πλc β 2 e -2 α inf t | g i (f )| 2 | s (f )| 2 ≥ 0 (32)
and the time-frequency correlation introduced in Definition 5 only depends on f :

ρ i,j (t, f, x j -x i ) = e ı(∠ g i (f )-∠ g j (f )) b(f,D) β2 . (33)
Finally, if ∀f ∈ R, ∠ g i (f ) = ∠ g j (f ), the asymptotic correlation function introduced in Definition 6 is

σ i,j (τ, x j -x i ) = b(τ,D) β2 ∈ [0, c 2D ], (34)
and function τ → σ i,j (τ, x j -x i ) is continuous in the interior of its support [-D c , D c ],
it reaches its maximum at τ = 0, and it is not differentiable at τ = 0.

Proposition 2 is proved in Appendix F. Note that Lemma 3 shows that when D → 0, the results in Proposition 2 come down to those already established in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF].

C. Early asymptotic state

Proposition 2 holds only in the case of diffuse acoustic fields and omnidirectional source and microphones. If we now focus on late reverberation (i.e. t → +∞), these assumptions can be relaxed: Proposition 3 (Early asymptotic state). Considering the stochastic reverberation model in Definition 1, suppose that ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant vector, therefore function α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α. For any sensors i, j ∈ {1 . . . I}, let γ i,j (t 1 , t 2 ) be the function defined in [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF]. Then the room response has the following statistical properties:

• Temporal domain: let ∀u ∈ S 2 , ∀t ∈ R,

m s * s (u, t) = E s (Θ(x, y)u, .) t * s (Θ(x, y)u, .) (35) ∀(x, y) ∈ R 3 × R M . Then function m s * s (u, t) is con- tinuous w.r.t. t ∈ R and twice continuously differentiable w.r.t. u ∈ S 2 . Moreover, ∀t ≥ T T g + T s , var[h i (t)] = λcβ 2 e -2 α inf t u∈S 2 g i (u, .) * g i (u, .) * m s * s (u, .) (0)du . (36) • Time-frequency domain: let ∀u ∈ S 2 , ∀f ∈ R, m | s | 2 (u, f ) = E s (Θ(x, y)u, f ) 2 (37) ∀(x, y) ∈ R 3 × R M . Then function m | s | 2 is smooth w.r.t. f ∈ R and twice continuously differentiable w.r.t. u ∈ S 2 . Moreover, ∀f ∈ R, when t → +∞, W γi,j (t, f ) = e -2 α inf t B i,j (f, x j -x i ) + O 1 t (38) where ∀f ∈ R, ∀r ∈ R 3 , B i,j (f, r) = λc u∈S 2 β qu r g i (u, f ) g j (u, f )m | s | 2 (u, f )e -2ıπf u r c du, (39)
and function f → B i,j (f, r) is smooth and even symmetric. In particular, if i = j, then ∀t ≥ T ,

W γi,i (t, f ) = λcβ 2 e -2 α inf t u∈S 2 | g i (u, f )| 2 m | s | 2 (u, f )du ≥ 0 (
40) and the time-frequency correlation in Definition 5,

ρ i,j (t, f, x j -x i ) = Bi,j (f,xj -xi) √ Bi,i(f,0)Bj,j (f,0) + O 1 t (41)
only depends on f when t → +∞, and it is asymptotically smooth and even symmetric w.r.t. f . • Asymptotic Gaussianity: when t → +∞, the sequence of random variables

hi(t) √ var[hi(t)]
converges in law to the standard Gaussian distribution.

Proposition 3 is proved in Appendix G. By introducing an additional assumption that holds in various experimental setups (cf. Section IV), the results in Proposition 3 can be simplified:

Corollary 1 (Early asymptotic state). Considering the stochastic reverberation model in Definition 1, suppose that ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant vector, therefore function α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α. Also suppose that the product

g i (u, f ) g j (u, f )m | s | 2 (u, f ) with m | s | 2 defined in (37) can be factorized as g i (u, f ) g j (u, f )m | s | 2 (u, f ) = ξ i,j (u) χ(f ), (42)
where function u → ξ i,j (u) ≥ 0 is twice continuously differentiable, and function f → χ(f ) ≥ 0 is smooth. For any sensors i, j ∈ {1 . . . I}, let D = x i -x j 2 denote the distance between the two sensors and let γ i,j (t 1 , t 2 ) be the function defined in [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF].

Then ∀f ∈ R, ∀t ≥ T T g + T s , W γi,i (t, f ) = λcβ 2 e -2 α inf t χ(f ) u∈S 2 ξ i,i (u)du. (43)
Besides, the asymptotic correlation function σ i,j (τ, x j -x i ) introduced in Definition 6 is nonnegative and it is continuous w.r.t. τ in the interior of its support [-D c , D c ]. If moreover function ξ i,j (.) is even on S 2 (symmetric case), then function τ → σ i,j (τ, x j -x i ) is also even, and its Fourier transform lim t→+∞ ρ i,j (t, f, x j -x i ) is real-valued.

If moreover function ξ i,j (.) is constant on S 2 (isotropic case), then when t → +∞,

ρ i,j (t, f, x j -x i ) = b(f,D) β2 + O( 1 t ) (44)
and σ i,j (τ, x j -x i ) = b(τ,D) β2 ∈ [0, c 2D ], (45)
and function τ → σ i,j (τ, x j -x i ) reaches its maximum at τ = 0, and it is not differentiable at τ = 0. If in addition D → 0, we get

lim t→+∞ ρ i,j (t, f, x j -x i ) = sinc 2πf D c + O(D) (46)
and

σ i,j (τ, x j -x i ) = c 2D + O(1) 1 [-D c , D c ] (τ ). ( 47
)
Corollary 1 is proved in Appendix G.5. Note that the additional assumption in Corollary 1 is necessary to conclude that the asymptotic correlation function σ i,j (τ, x j -x i ) is nonnegative and that its temporal support lies in [-D c , D c ]. Under the more general assumptions of Proposition 3, these properties may not hold.

D. Late asymptotic state

In Section III-C we assumed that the attenuation function α(.) in ( 7) is constant on S 2 . Now we aim to relax this assumption: Definition 7 (Regular attenuation function). Let α : S 2 → R + denote the restriction to S 2 of the function introduced in Definition 1 in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF]. Let α inf = inf u∈S 2 α(u) > 0, and suppose that function u → α(u) reaches its global minimum α inf at a finite set U = {u k } k∈K of distinct points u k ∈ S 2 . We further assume that the Hessian ¨ α k ∈ R of function α at every minimum point u k on the Riemannian manifold S 2 is positive.

Proposition 4 (Late asymptotic state). Considering the stochastic reverberation model in Definition 1, suppose that the attenuation function is regular (Definition 7). Let ∀u ∈ S 2 ,

• J q (u) = J q (u) Iuu T (48)

where J q (x) denotes the Jacobian matrix of function q(.) at x ∈ R 3 . For any sensors i, j ∈ {1 . . . I}, let γ i,j (t 1 , t 2 ) be the function defined in [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF]. Then the room response has the following statistical properties when t → +∞:

• Temporal domain:

var[h i (t)] = πλcβ2 e -2 α inf t t k∈K (g i (u k ,.) * g i (u k ,.) * m s * s (u k ,.))(0) √ ¨ α k + O( 1 t ) (49)
where

∀u k ∈ U, t → m s * s (u k , t) is the continuous function defined in (35). • Time-frequency domain: let ∀f ∈ R, ∀r ∈ R 3 , B i,j (f, r) = πλc k∈K β k (r) g i (u k ,f ) g j (u k ,f )m | s | 2 (u k ,f )e -2ıπf u k r c √ ¨ α k (50) where ∀k ∈ K, ∀r ∈ R 3 , β k (r) = β q (u k ) u T k + • J q (u k ) r (51)
and

∀u k ∈ U, f → m | s | 2 (u k , f ) is the smooth function defined in (37). Then function f → B i,j (f, r) is smooth and even symmetric. Moreover, ∀f ∈ R, W γi,j (t, f ) = e -2 α inf t t B i,j (f, x j -x i ) + O 1 √ t . (52) In particular, if i = j, W γi,i (t, f ) = πλcβ2e -2 α inf t t k∈K | g i (u k ,f )| 2 m | s | 2 (u k ,f ) √ ¨ α k + O( 1 t ) ≥ 0 (53)
and the time-frequency correlation in Definition 5,

ρ i,j (t, f, x j -x i ) = Bi,j (f,xj -xi) √ Bi,i(f,0)Bj,j (f,0) + O 1 √ t (54)
only depends on f when t → +∞, and it is asymptotically smooth and even symmetric w.r.t. f . • Asymptotic Gaussianity: the sequence of random variables

hi(t) √ var[hi(t)]
converges in law to the standard Gaussian distribution.

Proposition 4 is proved in Appendix G. By introducing the same additional assumption as in Corollary 1 that holds in various experimental setups (cf. Section IV), the results in Proposition 4 can be simplified:

Corollary 2 (Late asymptotic state). Considering the stochastic reverberation model in Definition 1, suppose that the attenuation function is regular (cf. Definition 7). Also suppose that the product

g i (u, f ) g j (u, f )m | s | 2 (u, f ) with m | s | 2 defined
in (37) can be factorized as in (42), where function u → ξ i,j (u) ≥ 0 is twice continuously differentiable, and function f → χ(f ) ≥ 0 is smooth. For any sensors i, j ∈ {1 . . . I}, let γ i,j (t 1 , t 2 ) be the function defined in [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF]. Then when t → +∞, ∀f ∈ R,

W γi,i (t, f ) = πλcβ 2 e -2 α inf t t χ(f ) k∈K ξi,i(u k ) √ ¨ α k + O( 1 t ) ,
(55) the time-frequency correlation introduced in Definition 5 is such that

ρ i,j (t, f, x j -x i ) = k∈K a k (x j -x i )e -2ıπf τ k (xj -xi) +O 1 √ t (56) where ∀k ∈ K, τ k (x j -x i ) = u k (xj -xi) c ∈ [-D c , D c ] and a k (r) = β k (r) β2 ξ i,j (u k ) √ ¨ α k l∈K ξ i,i (u l ) √ ¨ α l l∈K ξ j,j (u l ) √ ¨ α l ≥ 0,
where β k (r) is defined in (51). Finally, the temporal support of the asymptotic correlation function introduced in Definition 6 lies in [-D c , D c ]:

σ i,j (τ, x j -x i ) = k∈K a k (x j -x i ) δ (τ -τ k (x j -x i )) . ( 57
)
Corollary 2 is proved in Appendix G.5. Again, note that the additional assumption in Corollary 2 is necessary to conclude that the asymptotic correlation function σ i,j (τ, x j -x i ) is nonnegative and that its temporal support lies in [-D c , D c ]. Under the more general assumptions of Proposition 4, these properties may not hold.

E. Conclusions of the asymptotic analysis

Let us conclude Section III with a brief discussion about the various orders of convergence that appear in Propositions 3 and 4, and in Corollaries 1 and 2. First, in Proposition 3 and Corollary 1 in Section III-C, we observe that:

• var[h i (t)] and W γi,i (t, f ) are in closed-form ∀t ≥ T ; • if i = j, W γi,j (t, f ) and ρ i,j (t, f, x j -x i ) converge to
their asymptotic forms as O( 1 t ). Second, in Proposition 4 and Corollary 2 in Section III-D, we observe that:

• var[h i (t)] and W γi,i (t, f ) converge to their asymptotic forms as O( 1 t ); • if i = j, W γi,j (t, f ) and ρ i,j (t, f, x j -x i ) converge to their asymptotic forms as O( 1 √ t ).

Finally, we can make two additional observations:

• in all cases, the convergence to the asymptotic form is faster in Proposition 3 than in Proposition 4, which justifies the terms early and late asymptotic states; • the convergence speed of O( 1t ) is obtained in late asymptotic state for var[h i (t)] and W γi,i (t, f ), and in early asymptotic state for W γi,j (t, f ) and ρ i,j (t, f, x j -x i ).

This last observation is very important: in the case of real RIRs obtained from real measurements as in Section IV-B, the RIRs are observed on a limited time interval, because the measurement noise ends up dominating the reverberation after a while. It appears that a convergence speed of O( 1 √ t ) is too slow to permit the asymptotic form to be reached within the reverberation time. Therefore we will not be able to observe the asymptotic forms (56) and (57) of ρ i,j (t, f, x j -x i ) and σ i,j (τ, x j -x i ) in late asymptotic state. However, a convergence speed of O( 1t ) is fast enough to permit the asymptotic form to be observed within the reverberation time. This explains why in real conditions, the observed behavior of var[h i (t)] and W γi,i (t, f ) will match the predictions of the late asymptotic state, whereas the behavior of ρ i,j (t, f, x j -x i ) and σ i,j (τ, x j -x i ) will match the predictions of the early asymptotic state.

IV. EXPERIMENTAL RESULTS

In this section, we will study both synthetic and real room impulse responses 4 , and for each of them, we will check whether their asymptotic statistical behavior corresponds to one of the following states:

• diffuse acoustic field (characterized in Proposition 2);

• early asymptotic state of non-diffuse acoustic field (characterized in Proposition 3 and Corollary 1), in the particular symmetric or isotropic cases, or in the general anisotropic case; • late asymptotic state of non-diffuse acoustic field (characterized in Proposition 4 and Corollary 2). In order to perform this classification, we will consider the four following statistical signatures of reverberation5 :

• temporal power profile: var[h i (t)];

• time-frequency power profile: W γi,i (t, f ); • time-frequency correlation: ρ i,j (t, f, x j -x i ) in (10);

• asymptotic correlation function: σ i,j (τ, x j -x i ) in [START_REF] Kuttruff | Room Acoustics, Fifth Edition[END_REF]. In Table I, we summarize some results obtained in Propositions 2, 3 and 4 and in Corollaries 1 and 2. Remember that the results in Corollaries 1 and 2 hold when

g i (u, f ) g j (u, f )m | s | 2 (u, f ) can be factorized as ξ i,j (u) χ(f ), with ξ i,j (u) ≥ 0 and χ(f ) ≥ 0.
This assumption holds in various setups:

1) if the source response does not depend on frequency and if the microphones' responses are nonnegative and do not depend on frequency either, we get

ξ i,j (u) = g i (u) g j (u)m | s | 2 (u) ≥ 0 and χ(f ) = 1 (as illustrated in Section IV-A); 2)
Immediate (∀t ≥ T + D 2c ) var[h i (t)] ∝ e -2 α inf t Wγ i,i (t, f ) ∝ e -2 α inf t χ(f ) ρ i,j = b(f,D) β 2 -→ D→0 sinc( 2πf D c ) σ i,j = b(τ,D) β 2 ∼ D→0 c 2D 1 [-D c , D c ] (τ ) var[h i (t)] ∝ e -2 α inf t Wγ i,i (t, f ) ∝ e -2 α inf t χ(f ) O 1 t (fast) lim t→+∞ ρ i,j smooth, even symmetric; σ i,j ≥ 0, continuous within support ] -D c , D c [ lim t→+∞ ρ i,j real; σ i,j even ρ i,j -→ t→+∞ b(f,D) β 2 σ i,j = b(τ,D) β 2 var[h i (t)] ∝ e -2 α inf t t Wγ i,i (t, f ) ∝ e -2 α inf t t χ(f ) Observed be- fore RT 60 in: Shoebox room, Great Hall (Sections IV-A & IV-B3) Octagon room (Section IV-B2) Classroom (Section IV-B1) All rooms (Section IV) Shoebox room (Section IV-A) O 1 √ t (slow) ρ i,j -→ t→+∞ k∈K a k e -2ıπf τ k σ i,j = k∈K a k δ(τ -τ k ) Observed after RT 60 in: Shoebox room (Section IV-A) TABLE I THE FOUR STATISTICAL SIGNATURES OF REVERBERATION IN VARIOUS EXPERIMENTAL SETUPS and χ(f ) = | g (f )| 2 m | s | 2 (f ) ≥ 0 (as illustrated in Section IV-B1);
3) if the microphones are omnidirectional and have the same response g, and if

m | s | 2 (u, f ) ≥ 0 can be fac- torized as m | s | 2 (u, f ) = ξ s (u) χ s (f ) with ξ s (u) ≥ 0 and χ s (f ) ≥ 0, we get ξ i,j (u) = ξ s (u) ≥ 0 and χ(f ) = | g (f )| 2 χ s (f ) ≥ 0 (as illustrated in Sections IV-B2
and IV-B3). In the following sections, the four statistical signatures of reverberation are estimated from a pair (I = 2) of L observed room impulse responses h

(l) 1 (t) and h (l) 2 (t) for l ∈ {1 . . . L}, that are such that x 2 -x 1 is fixed. For i ∈ {1, 2}, the temporal power profile var[h i (t)] is thus estimated as 1 L L l=1 |h (l) i (t)| 2 . The Wigner dis- tributions W γi,i (t, f ) and W γ1,2 (t, f ) are estimated as W γi,i (t, f ) = 1 L L l=1 |S (l) hi (t, f )| 2 and W γ1,2 (t, f ) = 1 L L l=1 S (l) h1 (t, f )S (l) h2 (t, f ), where S (l) h1 (t, f ) (resp. S (l) h2 (t, f )) is the short time Fourier transform (STFT) of h (l) 1 (t) (resp. h (l) 2 (t)). The distributions W γi,i (t, f ) and W γ1,2 (t, f ) ob- tained in this way are smoothed estimates of W γi,i (t, f ) and W γ1,2 (t, f ) in the time-frequency domain [14]. Then the time- frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) is estimated as ρ 1,2 (t, f, x 2 -x 1 ) = W γ1,2 (t, f ) W γ1,1 (t, f ) W γ2,2 (t, f ) . The asymptotic correlation function σ 1,2 (τ, x 2 -x 1 ) is esti- mated as the inverse discrete Fourier transform (DFT) 6 of 1 |T | t∈T ρ 1,2 (t, f, x 2 -x 1 ) (58)
where

T is a time interval of length |T |, on which ρ 1,2 (t, f, x 2 -x 1 )
is approximately stationary (i.e. does not depend on t) on average.

A. Numerical simulation

We first considered synthetic RIRs generated by the Roomsimove toolbox [START_REF] Vincent | Roomsimove[END_REF], which is a state-of-the-art RIR generator based on the source image principle. Roomsimove is dedicated to parallelepipedic ("shoebox") rooms and applies high-pass filtering above 20 Hz. We used it with the default physical parameters (humidity: 40%, temperature: 20 • C, speed of sound: c = 343 m/s), and we removed the modeling of absorption due to the air 7 . We thus simulated a shoebox room having the same volume as the classroom described in [12, p. 84]: 200 m 3 (the room dimensions are 7.4 m × 9 m × 3 m). The values of the absorption coefficients for the six room surfaces are described in Table II.

TABLE II ABSORPTION COEFFICIENTS FOR THE SIX ROOM SURFACES

Let RT 60 denote the reverberation time, defined as the time the sound pressure level takes to reduce by 60 dB. For this setup, RT 60 = 0.23 s. We considered omnidirectional sources and directional (cardioid) microphones. The distance between the two microphones is 20 cm, and the vector x 2 -x 1 , pointing from the first microphone to the second one, is in the horizontal plane, forming an angle of 50 • from the xaxis and 40 • from the y-axis. We used the sensor orientations described in Table III We thus generated L = 1000 RIRs sampled at 16 kHz, with random source positions and random middle positions of the sensors (both uniformly distributed inside the room volume). We computed all the STFTs with a 128-sample-long Hann window and an overlap of 64 samples in the time domain.

Note that the particular numerical values of the parameters used in this simulation are only provided here for the sake of reproducible research. They were not chosen for their realism, but rather for reducing the computation time and for improving visualization. The observations that we will make below regarding the four signatures of reverberation would qualitatively be the same with different numerical values.

Besides, in the Roomsimove toolbox, the source response does not depend on frequency and the microphones' responses are nonnegative and do not depend on frequency either, thus ∀i, j ∈ {1 . . . I}, g i (u, f )

g j (u, f )m | s | 2 (u, f ) can be factorized as in Corollaries 1 and 2 as ξ i,j (u) χ(f ), with ξ i,j (u) = g i (u) g j (u)m | s | 2 (u) ≥ 0 and χ(f ) = 1.
Thus we expect the measurements on the data to match the results in Table I with χ(f ) = 1 ∀f ∈ R. ] in dB (blue curve), superimposed with a straight red line obtained by linear regression. We can observe that the temporal power profile is slightly bent compared with the straight line, so we are neither in a diffuse acoustic field, nor in the early asymptotic state of a non-diffuse acoustic field in Table I (otherwise, we would have ln(var[h i (t)]) c = -2 α inf t). Fig. 1-(b) represents t var[h i (t)] in dB (blue curve), and a straight red line obtained by linear regression. This time the blue curve is not bent and matches the output of linear regression much better. This corresponds to the behavior predicted by the late asymptotic state of a non-diffuse acoustic field in Table I:

ln(t var[h i (t)]) c = -2 α inf t.
Fig. 2. Time-frequency power profile Wγ i,i (t, f ) in the shoebox room Fig. 2 represents the time-frequency power profile W γi,i (t, f ) in dB. We note that it can be approximately factorized as a function of time multiplied by a fixed spectrum, which is the behavior predicted by the reverberation model in all asymptotic states in Table I. Due to the projection property (3), the function of time is necessarily proportional to the temporal power profile: e -2 α inf t t , which corresponds to the late asymptotic state. Besides, we observe that the fixed spectrum is approximately constant, which is again the behavior predicted by the reverberation model in all asymptotic states, since χ(f ) = 1 ∀f ∈ R. Finally, we can conclude that the time-frequency power profile is the one predicted by the late asymptotic state of a non-diffuse acoustic field (cf. Table I).

Fig. 3 (resp. Fig. 4) represents the real part (resp. the imaginary part) of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ). A very interesting phenomenon can be observed: the timefrequency correlation is not stationary, its "spectrum" evolving from a shape similar to that of a cardinal sine (at low values of t in Fig. 3), to the shape of a sine wave (at high values of t in Fig. 3 and4). Note that the time axis goes up to 0.8 s, far beyond the RT 60 = 0.23 s (we have deliberately modified Roomsimove's code in order to synthesize such long RIRs). Fig. 5 and 6 will help us understand what is at stake here.

The blue curve in Fig. 5-(a) (resp. Fig. 5-(b)) represents the real part (resp. the imaginary part) of the last spectrum of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in Fig. 3 (resp. Fig. 4). The red curves in these figures represent the asymptotic time-frequency correlation predicted by the late I). Indeed, it appears that the minimum of function u → α(u) is reached at a single direction u 0 . Consequently, the real (resp. imaginary) part of lim t→+∞ ρ 1,2 (t, f, x 2 -x 1 ) in ( 56) is a cosine function (resp. a sine function). The good match between the ground truth and the estimation in Fig. 5 shows that the time-frequency correlation in Fig. 3 and4 does converge to the behavior predicted by the late asymptotic state of a non-diffuse acoustic field 8 . 8 Since the acoustic field in this rectangular room is not isotropic, it is not diffuse. The reader might notice that yet, we used a simulated shoebox room in the experimental section of [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] to illustrate the properties of a diffuse acoustic field. Actually, we made the experimental setup isotropic in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] by randomizing the orientation of vector x

(l) 2 -x (l)
1 according to a uniform probability distribution, and by averaging the results. vertical lines), as predicted by the reverberation model in all asymptotic states (cf. Table I). Moreover, within this temporal support (Fig. 7-(b) represents a zoom in on this interval), the correlation function is nonnegative and not even, and it has continuous variations, as predicted by the anisotropic case of the early asymptotic state of a non-diffuse acoustic field in Table I.

We can conclude that up to the RT 60 , the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) represented in Fig. 3 and4 actually behaves as in the early asymptotic state, whereas the late asymptotic state is actually reached much later. This makes a difference with the temporal power profile var[h i (t)] in Fig. 1 and the time-frequency power profile W γi,i (t, f ) in Fig. 2, which match the predictions of the late asymptotic state almost right from the beginning.

This remark confirms our conclusions in Section III-E regarding the asymptotic analysis in Section III. When we will analyze real RIRs in Section IV-B, obtained from measurements in real acoustic environments, we will not be able to estimate the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) beyond the reverberation time, because the measurement noise ends up dominating the reverberation. Therefore we will never be able to observe the late asymptotic state in the timefrequency correlation, which stays stuck in early asymptotic state in its temporal observation interval. This explains why in real conditions, the behavior of the first two reverberation signatures (temporal power profile var[h i (t)] and time-frequency power profile W γi,i (t, f )) will match the predictions of the late asymptotic state, whereas the last two signatures (timefrequency correlation ρ 1,2 (t, f, x 2 -x 1 ) and asymptotic correlation function σ 1,2 (τ, x 2 -x 1 )) will match the predictions of the early asymptotic state.

B. Real RIR measurements

We used the collection of room impulse responses measured in a classroom, the Octagon room, and the Great Hall at the Mile End campus of Queen Mary, University of London in 2008 [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF]. The measurements were created using a sine sweep technique [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept sine technique[END_REF] with a Genelec 8250A loudspeaker (which forms a directive sound source) and two microphones, an omnidirectional DPA 4006 and a B-format Soundfield SPS422B. Each measurement has source and receiver heights of 1.5 m. We used the RIRs measured with the omnidirectional microphone. We resampled the RIRs at 48 kHz, and we truncated them so as to remove both their beginning (formed of early reflections) in order to focus on late reverberation only, and their end, which is dominated by the measurement noise. All STFTs are computed by means of Hann windows of length 2400 samples, with a 95% overlap.

As we will show in the following subsections, the acoustic field in the three rooms is uniform and non-diffuse. Even though we did not report here the results in order to avoid overloading this document, the uniformity of the acoustic field over space was checked by computing the four reverberation signatures (temporal power profile var[h i (t)], timefrequency power profile W γi,i (t, f ), time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ), and asymptotic correlation function σ 1,2 (τ, x 2 -x 1 )) through space averages in different parts of the rooms (we separated the half-left side from the half-right side, and the half-up side from the half-bottom side, along the two horizontal directions). We thus observed a remarkably accurate match of the four signatures in the different parts of the room, which permitted us to conclude that the acoustic field is uniform (at least in the space areas where the measurements were carried out). In the following subsections, the reported signatures are computed through space averages over all available measurements in the rooms.

1) Classroom: In the QMUL Classroom Impulse Response "Omni" dataset, 130 RIRs were measured within a classroom. As described in [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF], the room measures 7.5 × 9 × 3.5 m (236 m 3 ) with reflective surfaces of a linoleum floor, painted plaster walls and ceiling, and a large whiteboard. Measurements were 50 cm apart arranged in 10 rows and 13 columns (over a 9 m × 12 m area) relative to the speaker, with the 7th column directly on axis with the speaker. For this setup, the average measured reverberation time RT 30 is about 1.8 s around 1000 Hz [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF]Fig. 4]. For the correlation measurements (time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) and asymptotic correlation function σ 1,2 (τ, x 2 -x 1 )), we computed space averages over all pairs of microphones placed 50 cm apart, and such that the vector x 2 -x 1 , pointing from the first microphone to the second one, is in the horizontal plane, in the direction of the y-axis.

Fig. 8 represents the temporal power profile var[h i (t)]. As for the simulated shoebox room (cf. Fig. 1), we observe that its behavior is predicted by the late asymptotic state of a nondiffuse acoustic field: var[h i (t)] does not decrease as e -2 α inf t , but rather as e -2 α inf t t (cf. Table I). Fig. 9 represents the time-frequency power profile W γi,i (t, f ). We observe that this power profile can no longer be approximately factorized as a function of time multiplied by a fixed spectrum: the temporal decrease rate depends on frequency, which is not predicted by the reverberation model studied in this research report. Indeed, remember that we simplified the general stochastic reverberation model introduced in [START_REF] Badeau | General stochastic reverberation model[END_REF] by removing the dependency of vector α on frequency f , in order to simplify the mathematical analysis of the model (cf. Sections I and II). Actually, the general model introduced in [START_REF] Badeau | General stochastic reverberation model[END_REF] does account for a frequency-varying decrease rate, but its mathematical analysis is left for future work. Fig. 10 (resp. Fig. 11) represents the real part (resp. the imaginary part) of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ). We notice that the imaginary part is zero on average (the surface in Fig. 11 looks like centered noise), and that the real part is approximately stationary (it does not depend on time), as predicted by the reverberation model in I).

The blue curve in Fig. 12-(a) (resp. Fig. 12-(b)) represents the average over time of the real part (resp. the imaginary part) of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in Fig. 10 (resp. Fig. 11). The red curve in Fig. 12-(a with c = 343 m/s and D = 50 cm, whereas the red curve in Fig. 12-(b) represents the zero prediction. We thus observe that the shape of the timefrequency correlation is well predicted by the isotropic case of the early asymptotic state of a non-diffuse acoustic field when D → 0 (cf. Table I).

Indeed, the results of Corollary 1 hold in this experimental setup, because in this dataset, the microphones are omnidirectional and have the same response g. In addition, if we assume The blue curve in Fig. 13-(a) represents the asymptotic correlation function σ 1,2 (τ, x 2 -x 1 ). We observe that the temporal support of this correlation function lies in the interval [-D c , D c ] (whose boundaries are represented by red vertical lines), as predicted by the reverberation model in all asymptotic states. Moreover, within this temporal support (Fig. 13-(b) represents a zoom in on this interval), the correlation function is nearly constant, as predicted in the isotropic case of the early asymptotic state when D → 0, which confirms what we already observed in Fig. 12. However, the assumption D → 0 does not seem to hold completely, because in Fig. 13-(a), the shape of the asymptotic correlation function

m | s | 2 is also omnidi- rectional, then ∀i, j ∈ {1 . . . I}, g i (u, f ) g j (u, f )m | s | 2 (u, f ) can be factorized as ξ i,j (u) χ(f ), with ξ i,j (u) = 1 and χ(f ) = | g (f )| 2 m | s | 2 (f ) ≥ 0,
σ 1,2 (τ, x 2 -x 1 ) = b(τ,D) β2
is rather peaky at τ = 0. This behavior confirms the prediction of Lemma 3, which states that function τ → b(τ, D) reaches its maximum at τ = 0, and is not differentiable at τ = 0.

We repeated the correlation measurements represented in Fig. 13 with different experimental setups: we tested various distances between microphones, and orientations of the vector x 2 -x 1 , pointing from the first microphone to the second one. We did not include here all the results in order to avoid overloading this document, but all experimental setups led to figures looking like Fig. 13. This confirms that, regarding the time-frequency correlation between sensors, the reverberation in the classroom behaves as in the isotropic case of the early asymptotic state of a non-diffuse acoustic field.

2) Octagon room: In the QMUL Octagon Impulse Response "Omni" dataset, 169 RIRs were measured in the center of the Octagon room (a Victorian building completed in 1888 and originally designed to be a library). As described in [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF], the walls of the room are lined with books, with a wooden floor and plaster ceiling. The room has eight walls each 7.5 m in length and a domed ceiling reaching 21 m over the floor, with an approximate volume of 9500 m 3 . Measurements were 1 m apart arranged in 13 rows and 13 columns (over a 12 m × 12 m area) relative to the speaker, with the 7th column directly on axis with the speaker. For this setup, the average measured reverberation time RT 30 is about 2 s around 1000 Hz [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF]Fig. 4]. For the correlation measurements (time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) and asymptotic correlation function σ 1,2 (τ, x 2 -x 1 )), we computed space averages over all pairs of microphones placed 1 m apart, and such that the vector x 2 -x 1 , pointing from the first microphone to the second one, is in the horizontal plane, in the direction of the x-axis.

Fig. 14 to Fig. 17 respectively represent the temporal power profile var[h i (t)], the time-frequency power profile W γi,i (t, f ), and the real and imaginary parts of the timefrequency correlation ρ 1,2 (t, f, x 2 -x 1 ). All observations made in Section IV-B1 still hold here.

Regarding the asymptotic correlation function σ 1,2 (τ, x 2x 1 ) represented in Fig. 18-(a), the observations are a bit different. First, the temporal support of this correlation function still essentially lies in the interval [-D c , D c ] (whose boundaries are represented by red vertical lines). Besides, Fig. 18-(b) shows that within this interval, the correlation function is still nonnegative and has continuous variations, but it can no longer be considered as approximately constant. Instead, it looks approximately even, which suggests that function u → ξ 1,2 (u) may also be even, as stated in Corollary 1. Therefore, contrary to what we observed in the classroom, the observations in the Octagon room can no longer by explained by the isotropic case, but rather by the symmetric case of the early asymptotic state of a non-diffuse acoustic field (cf. Table I).

Indeed, the results of Corollary 1 hold in this experimental setup, because in this dataset, the microphones are omnidirec- 18-(a) has a particular shape that cannot be modeled as centered noise, whereas Corollary 1 predicts that it should be zero. This discrepancy may be explained by the inaccuracy of the approximation

m | s | 2 (u, f ) ≈ ξ s (u) χ s (f ).
We repeated the correlation measurements represented in 3) Great hall: In the QMUL Great Hall Impulse Response "Omni" dataset, 169 RIRs were measured within the Great Hall (a multipurpose hall that can hold approximately 800 seats). As described in [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF], the hall has a stage and seating areas on the floor and a balcony. The microphones were placed in the seating area on the floor, approximately a 23 m × 16 m area, but the room is significantly bigger as the balcony extends 20 m past the rear wall. Measurements were 1 m apart arranged in 13 rows and 13 columns (over a 12 m × 12 m area) relative to the speaker, with the 7th column directly on axis with the speaker. For this setup, the average measured reverberation time RT 30 is about 2 s around 1000 Hz [16, Fig. 4]. For the correlation measurements (time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) and asymptotic correlation function σ 1,2 (τ, x 2 -x 1 )), we computed space averages over all pairs of microphones placed 1 m apart, and such that the vector x 2 -x 1 , pointing from the first microphone to the second one, is in the horizontal plane, in the direction of the y-axis. Fig. 19 to 22 respectively represent the temporal power profile var[h i (t)], the time-frequency power profile W γi,i (t, f ), and the real and imaginary parts of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ). All observations made in Section IV-B1 still hold here. Note however that the imaginary part of the time-frequency correlation in Fig. 22 is not approximately zero at low frequencies, thus the real-valued asymptotic correlation function in Fig. 23 is not even, as in the case of the shoebox room (cf. Section IV-A, Fig 6). Therefore function ξ 1,2 (u) is not even either, as showed in Corollary 1. This behavior corresponds to the anisotropic case of the early asymptotic state of a non-diffuse acoustic field (cf. Table I).

V. CONCLUSION AND PERSPECTIVES

In this research report, we extended the mathematical analysis of the stochastic model initially proposed in [START_REF] Badeau | Unified stochastic reverberation modeling[END_REF], [START_REF]Common mathematical framework for stochastic reverberation models[END_REF], which was dedicated to the particular case of diffuse acoustic fields, omnidirectional sources and microphones, and constant attenuation w.r.t frequency, to the more general case of uniform The next step in this research project will be the extension of both the mathematical analysis and the experimental validation to the more realistic case of frequency-varying attenuation coefficients, before addressing the most general case of nonuniform acoustic fields. Our purpose is then to develop efficient algorithms for estimating the model parameters, in order to investigate the potential of the general stochastic reverberation model in various signal processing applications. In this appendix, we introduce a common formulation of the two equations ( 5) and ( 6) in Definition 1, which are expressed in the temporal and the spectral domains, respectively. This will permit us to prove general results that hold in both domains.

A.1. Common formulation of temporal and spectral equations Lemma 4. Both (5) and (6) in Definition 1 can be written in the form

h i = x∈R 3 y∈R M Vi(x,y;y-q(x-xi))ϕi(x,y;x-xi)e -y α c x-xi 2 dN (x, y), (59)
ϕ i ∞ = sup R 3 ×R M ×R 3 |ϕ i (x, y; r)| < +∞.
Proof. With r = xx i , ( 5) is derived by applying (59) to 6) is derived by applying (59) to

ϕ i (x, y; r) = g i r r 2 , . t * s Θ(x, y) r r 2 , t -r 2 c , (60)
ϕ i (x, y; r) = g i r r 2 , f s Θ(x, y) r r 2 , f e -2ıπf r 2 c , (61) where ∀(x, y) ∈ R 3 × R M , ∀r ∈ R 3 , |ϕ i (x, y; r)| ≤ sup u∈S 2 ,f ∈R | g i (u, f ) | sup u∈S 2 ,f ∈R | s (u, f ) |.

A.2. First order moments

Proposition 5. Considering the stochastic reverberation model in Definition 1, the expected value of (59) for any sensor i ∈ {1 . . . I} is

E[h i ] = λβ 1 r∈R 3 mϕ i (r)e -α(r) c r 2 dr < +∞ (62)
where function α(.) is defined in (7), β 1 is defined in (13) for n = 1, and ∀x,

y ∈ R 3 × R M , ∀r ∈ R 3 , m ϕi (r) = E[ϕ i (x, y; r)] ∈ L ∞ (R 3 ). (63)
Proof. First, we note that β 1 is finite, as proved in Lemma 1.

We also note that m ϕi ∈ L ∞ (R 3 ) because ϕ i ∞ < +∞, and α inf = inf u∈S 2 α(u) > 0, so that ∀r ∈ R 3 ,

m ϕi (r)e -α(r) c r 2 ≤ m ϕi ∞ e -α inf r c r 2 ,
which proves that the integral in (62) converges.

Let us now prove (62). The expected value of (59) can be written as

E[h i ] = E[I],
where the mathematical expectation is w.r.t. the three random fields ϕ i , V i and dN on R 3 × R M , with I = x∈R 3 y∈R M ψ(x, y)dN (x, y), and ψ(x, y) = Vi(x,y;y-q(x-xi))ϕi(x,y;x-xi)e -y α

c x-xi 2 . ( 64
)
By applying Proposition 4 in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] to Λ(x, y) = λ and I defined above, equation (B4) in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] yields

E[I] = λE x∈R 3 y∈R M ψ(x, y)dxdy , ( 65
)
where the mathematical expectation is now w.r.t. the two random fields ϕ i and V i only. Substituting (64) into (65) yields

E[h i ] = λ E x∈R 3 y∈R M Vi(x,y;y-q(x-xi))ϕi(x,y;x-xi)e -y α c x-xi 2
dxdy .

(66) With the change of variable r = xx i , which is such that dr = dx, since E[V i (x, y; yq(r))] = p(1; yq(r)) and E[ϕ i (x, y; r)] = m ϕi (r) as defined in (63), (66) yields

E[h i ] = λ r∈R 3 y∈R M p(1;y-q(r))mϕ i (r)e -y α c r 2 drdy. (67)
Finally, with the change of variable z = yq(r), which is such that dz = dy, and by substituting [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] and [START_REF] Badeau | General stochastic reverberation model[END_REF] for n = 1 into (67), we get (62).

A.3. Second order moments

Proposition 6. Considering the stochastic reverberation model in Definition 1, the covariance of (59) for any sensors i, j ∈ {1 . . . I} is

cov[h i , h j ] = λ x∈R 3 β(q(x -x i ) -q(x -x j )) mϕ i ,ϕ j (x-xi,x-xj )e - α(x-x i )+ α(x-x j ) c x-xi 2 x-xj 2 dx < +∞ ( 68
)
where function α(.) is defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF], function β(.) is defined in [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept sine technique[END_REF], and

∀(x, y) ∈ R 3 × R M , ∀r i , r j ∈ R 3 , m ϕi,ϕj (r i , r j ) = E ϕ i (x, y; r i )ϕ j (x, y; r j ) ∈ L ∞ (R 3 ×R 3 ). ( 69
)
Proof. First, Lemma 2 shows that function β(.) in ( 17) is upper bounded by β 2 ; however β 2 is finite, as proved in Lemma 1. We also note that

m ϕi,ϕj ∈ L ∞ (R 3 × R 3 ) because ϕ i ∞ < +∞, and α inf = inf u∈S 2 α(u) > 0, so that ∀x ∈ R 3 , β(q(x-xi)-q(x-xj ))mϕ i ,ϕ j (x-xi,x-xj )e - α(x-x i )+ α(x-x j ) c x-xi 2 x-xj 2 ≤ β2 mϕ i ,ϕ j ∞e - α inf ( x-x i 2 + x-x j 2 ) c x-xi 2 x-xj 2 ,
which proves that the integral in (68) converges.

Let us now prove (68). For k ∈ {1, 2}, let

I k = x∈R 3 y∈R M ψ k (x, y)dN (x, y),
with

ψ 1 (x, y) = Vi(x,y;y-q(x-xi))ϕi(x,y;x-xi)e -y α c x-xi 2 , ψ 2 (x, y) = Vj (x,y;y-q(x-xj ))ϕj (x,y;x-xj )e -y α c x-xj 2 . ( 70
) We note that E[I 1 I 2 ] = cov[I 1 , I 2 ] + E[I 1 ]E[I 2 ], (71)
where the mathematical expectation is w.r.t. the random fields ϕ i , ϕ j , V i , V j and dN on R 3 × R M . By applying Proposition 4 in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] to Λ(x, y) = λ and I 1 and I 2 defined above, equation (B4) in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] yields

E[I 1 ] = λE x∈R 3 y∈R M ψ 1 (x, y)dxdy , E[I 2 ] = λE x∈R 3 y∈R M ψ 2 (x, y)dxdy , (72)
where the mathematical expectations in the right members of these equalities are now w.r.t. the random fields ϕ i , ϕ j and V i , V j only. In the same way, by applying Proposition 5 in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] to Λ(x, y) = λ and I 1 and I 2 defined above, equation (B8) 9in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] yields

cov[I 1 , I 2 ] = λE x∈R 3 y∈R M ψ 1 (x, y)ψ 2 (x, y)dxdy .
(73) By substituting (72) and ( 73) into (71), we get

E[I 1 I 2 ] = λE x∈R 3 y∈R M ψ 1 (x, y)ψ 2 (x, y)dxdy +λ 2 E x1∈R 3 y 1 ∈R M ψ 1 (x 1 , y 1 )dx 1 dy 1 x2∈R 3 y 2 ∈R M ψ 2 (x 2 , y 2 )dx 2 dy 2 .
(74) However, since the random increments ψ 1 (x 1 , y 1 ) and

ψ 2 (x 2 , y 2 ) are independent when (x 1 , y 1 ) = (x 2 , y 2 ), we get ∀x 1 , x 2 ∈ R 3 , ∀y 1 , y 2 ∈ R M , λ 2 E ψ 1 (x 1 , y 1 )ψ 2 (x 2 , y 2 ) = λ 2 cov[ψ 1 (x 1 , y 1 ), ψ 2 (x 2 , y 2 )]δ x1,x2 δ y 1 ,y 2 +λ 2 E [ψ 1 (x 1 , y 1 )] E ψ 2 (x 2 , y 2 ) , (75)
where δ x,y denotes the Kronecker delta: δ x,y = 1 if x = y, or δ x,y = 0 if x = y. Since the Lebesgue measure of the support

of δ x1,x2 δ y 1 ,y 2 is zero in (R 3 × R M ) 2 , integrating (75) over (R 3 × R M ) 2 yields λ 2 x1,y 1 ,x2,y 2 E ψ 1 (x 1 , y 1 )ψ 2 (x 2 , y 2 ) dx 1 dy 1 dx 2 dy 2 = λ 2 x1,y 1 ,x2,y 2 E[ψ 1 (x 1 , y 1 )]E[ψ 2 (x 2 , y 2 )]dx 1 dy 1 dx 2 dy 2 = E[I 1 ]E[I 2 ],
(76) where we have used (72). By substituting (74) and ( 76) into (71), we get

cov[I 1 , I 2 ] = E[I 1 I 2 ] -E[I 1 ]E[I 2 ] = λE x∈R 3 y∈R M ψ 1 (x, y)ψ 2 (x, y)dxdy . (77) Since ∀(x, y) ∈ R 3 × R M , E [V i (x, y; y -q(x -x i ))V j (x, y; y -q(x -x j ))] = p(1, 1; y -q(x -x i ), y -q(x -x j ))
and

E ϕ i (x, y; x-x i )ϕ j (x, y; x-x j ) = m ϕi,ϕj (x-x i , x-x j )
as defined in (69), by substituting (70) into (77), we get

cov[I 1 , I 2 ] = λ x∈R 3 y∈R M p(1,1;y-q(x-xi),y-q(x-xj ))mϕ i ,ϕ j (x-xi,x-xj )e -2y α c x-xi 2 x-xj 2 dxdy.
(78) Finally, with the change of variable z = y-

q(x-xi)+q(x-xj ) 2
, which is such that dz = dy, and by substituting ( 7) and ( 17) into (78), we get (68).

A.4. Characteristic function

Proposition 7. Considering the stochastic reverberation model in Definition 1, for any sensors i, j ∈ {1 . . . I}, the characteristic function

φ hi,hj (θ 1 , θ 2 ) = E e ı(Re(θ1hi+θ2hj )) ∀θ 1 , θ 2 ∈ C (79)
of the random vector [h i , h j ] with h i and h j defined in (59) is such that

ln φ hi,hj (θ 1 , θ 2 ) = λ x∈R 3 y∈R M p(1, 1; y -q(x -x i ), y -q(x -x j )) φ ϕi(x-xi),ϕj (x-xj ) θ1e -y α c x-xi 2 , θ2e -y α c x-xj 2 -1 +p(1, 0; y -q(x -x i ), y -q(x -x j )) φ ϕi(x-xi) θ1e -y α c x-xi 2 -1 +p(0, 1; y -q(x -x i ), y -q(x -x j )) φ ϕj (x-xj ) θ2e -y α c x-xj 2 -1 dxdy. ( 80
)
Proof. By substituting (59) into (79) we get

φ hi,hj (θ 1 , θ 2 ) = E[e ıθI ],
where the mathematical expectation is w.r.t. the random fields ϕ i , ϕ j , V i , V j and dN on R 3 × R M , with θ = 1, I = (81) By applying Proposition 4 in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] to Λ(x, y) = λ, and θ, I defined above, equation (B1) in [START_REF]Common mathematical framework for stochastic reverberation models[END_REF] yields

φ hi,hj (θ 1 , θ 2 ) = E e λ x∈R 3 y∈R M (e ıψ(x,y) -1)dxdy = +∞ n=0 1 n! E λ x∈R 3 y∈R M (e ıψ(x,y) -1)dxdy n , ( 82
)
where the mathematical expectation is now w.r.t. the random fields ϕ i , ϕ j , V i , and V j only. Since the random increments ψ(x 1 , y 1 ) and ψ(x 2 , y 2 ) are independent when (x 1 , y 1 ) = (x 2 , y 2 ), we get

E λ x∈R 3 y∈R M (e ıψ(x,y) -1)dxdy n = E λ x∈R 3 y∈R M (e ıψ(x,y) -1)dxdy n (83)
because when developing the product in the left member, all terms but the one in the right member are zero, since they result from integrations on Borel sets whose Lebesgue measure is zero in (R 3 × R M ) n (we use here the same argument as in the proof of Proposition 6). Substituting (83) into (82) yields (84) By substituting (81) into (84), we get:

φ hi,hj (θ 1 , θ 2 ) = +∞ n=0 1 n! E λ
ln φ hi,hj (θ 1 , θ 2 ) = λ x∈R 3 y∈R M E e ıRe θ1
Vi(x,y;y-q(x-xi))ϕi(x,y;x-xi)

x-xi 2 e -y α c e ıRe θ2 Vj (x,y;y-q(x-xj ))ϕj (x,y;x-xj )

x-xj 2 e -y α c -1dxdy.

Then, by considering the conditional expectation given V i (x, y; yq(xx i )) and V j (x, y; yq(xx j )), we get

ln φ hi,hj (θ 1 , θ 2 ) = λ x∈R 3 y∈R M p(1, 1; y -q(x -x i ), y -q(x -x j )) E e ıRe θ1 ϕi(x,y;x-xi) x-xi 2 +θ2 ϕj (x,y;x-xj ) x-xj 2 e -y α c -1 +p(1, 0; y -q(x -x i ), y -q(x -x j )) E e ıRe θ1 ϕi(x,y;x-xi) x-xi 2 e -y α c -1 +p(0, 1; y -q(x -x i ), y -q(x -x j )) E e ıRe θ2 ϕj (x,y;x-xj ) x-xj 2 e -y α c -1 dxdy,
where the mathematical expectation is now w.r.t. the random fields ϕ i and ϕ j only, which finally proves (80).

Note that Proposition 7 provides a straightforward proof for Proposition 1 in [START_REF] Badeau | General stochastic reverberation model[END_REF]. Indeed, equation ( 8) in [START_REF] Badeau | General stochastic reverberation model[END_REF] can be proved by applying Proposition 7 to θ 1 = θ, θ 2 = 0, α = 0, x i = 0, and ∀x,

y ∈ R 3 × R M , ϕ i (x, y; r) = r 2 ψ(r) ∈ L ∞ (R 3 × R M × R 3 ) since ψ ∈ L ∞ (R 3 ) has compact support.

APPENDIX B GEOMETRY WITH TWO MICROPHONES

B.1. Computation of integrals over space

We consider two sensor locations x i , x j ∈ R 3 and D = x i -x j 2 . In Appendices E to G, we will need to compute several integrals of the form:

J ξ = x∈R 3 ξ(x-xi,x-xj ) x-xi 2 x-xj 2 dx, (85)
where

x → ξ (x -x i , x -x j ) ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ).
In the particular case i = j, let r = xx i and

J ξ = r∈R 3 ξ(r) r 2 2 dr, (86)
where

ξ ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ). x j x i x θ r 0 D 2 ρ + vD 2 ρ -vD 2 D 2 θ i θ j
Fig. 24. Geometry with two microphones at x i , x j and a source image at x.

To compute such integrals, we will use the spherical coordinates (θ, ϕ, r), as illustrated in Fig. 24, where θ = 0 corresponds to the direction of vector x j -x i , and the origin of the coordinates is the middle of the line segment [x i , x j ] (so that xi+xj 2

= 0 as represented in Fig. 24). We thus get x = [r sin(θ) cos(ϕ), r sin(θ) sin(ϕ), r cos(θ)] and dx = r 2 dr sin(θ)dθdϕ, with r ∈ R + , θ ∈ [0, π] and ϕ ∈ [0, 2π]. Moreover, as can be deduced from Fig. 24

, • x -x i 2 = r 2 + D 2 4 + rD cos(θ), • x -x j 2 = r 2 + D 2 4 -rD cos(θ),
• θ i = acos r cos(θ)+D/2

r 2 + D 2 4 +rD cos(θ) , • θ j = acos r cos(θ)-D/2 r 2 + D 2 4 -rD cos(θ)
.

By substitution into (85), we get

J ξ = +∞ r=0 π θ=0 2π ϕ=0 ξ θi,ϕ, r 2 + D 2 4 +rD cos(θ),θj ,ϕ, r 2 + D 2 4 -rD cos(θ) r 2 + D 2 4 +rD cos(θ) r 2 + D 2 4 -rD cos(θ)
r 2 dr sin(θ)dθdϕ.

(87) Note that in (87), the two arguments of function ξ have been replaced by their spherical coordinates. In the rest of this research report, we will continue using this notation without any further notice, as long as it is not ambiguous.

Finally, we make a last change of variables, that is also illustrated in Fig. 24:

ρ = x-xi 2 + x-xj 2 2 = r 2 + D 2 4 +rD cos(θ)+ r 2 + D 2 4 -rD cos(θ) 2 , v = x-xi 2-x-xj 2 D = r 2 + D 2 4 +rD cos(θ)-r 2 + D 2 4 -rD cos(θ) D , (88) which is such that ρ ∈ D 2 , +∞ , v ∈ [-1, 1], θ i = acos 2ρv+D
2ρ+vD , θ j = acos 2ρv-D 2ρ-vD , and r 2 dr sin(θ)dθ

r 2 + D 2 4 +rD cos(θ) r 2 + D 2 4 -rD cos(θ) = dρdv. (89)
Indeed, we note that (88) implies that r cos(θ) = ρv and

r 2 = ρ 2 -(1 -v 2 ) D 2 4 , thus x =     ρ 2 -D 2 4 √ 1 -v 2 cos(ϕ) ρ 2 -D 2 4 √ 1 -v 2 sin(ϕ) ρv     , (90)
and dx = ρ 2 -( vD 2 ) 2 dρdvdϕ, which proves (89). Therefore substituting (88) and (89) into (87) yields

J ξ = +∞ ρ= D 2 1 v=-1 2π ϕ=0 ξ 2ρv+D 2ρ+vD , ϕ, ρ + vD 2 , 2ρv-D 2ρ-vD , ϕ, ρ -vD 2 dρdvdϕ.
(91) Again, note that compared with (87), in (91) the angles θ i and θ j have been replaced by their cosine. In the rest of this research report, we will continue using this notation without any further notice, as long as it is not ambiguous since θ i , θ j ∈ [0, π] and the cosine function takes values in

[-1, 1].
Finally, in the particular case i = j (i.e. D = 0), (91) shows that (86) is equivalent to

J ξ = +∞ r=0 1 v=-1 2π ϕ=0 ξ (v, ϕ, r) drdvdϕ. (92)
B.2. Differentiation of function q(.) on S 2

In Appendix G, we will need to express the quantity β q (v, ϕ) vD + ∂q(v,ϕ) ∂v (1 -v 2 )D as a function of u x x 2 ∈ S 2 when ρ → +∞ (which is equivalent to r → +∞), where q (v, ϕ) denotes q(u). To that end, we first note that (90) implies x (x j -x i ) = ρvD, therefore

vD = u (x j -x i ) 1 + O( D 2 ρ 2 ) . (93)
Moreover, ∀x ∈ R 3 , ∂q ∂v (x) = J q (x) ∂x ∂v , where J q (x) denotes the Jacobian matrix of function q(.) at x ∈ R 3 , and (90) implies

∂x ∂v =     -ρ 2 -D 2 4 v √ 1-v 2 cos(ϕ) -ρ 2 -D 2 4 v √ 1-v 2 sin(ϕ) ρ     .
(94) However, by noting that function q(.) is 1-homogeneous, we have J q (x) = J q (u), since J q (x) does not depend on x 2 . Therefore we can write

∂q ∂v (x) = J q (u) ∂x ∂v . (95)
We note that (90) and (94) imply

1 ρ (1 -v 2 )D ∂x ∂v + vDx =   0 0 D   = x j -x i .
Substituting (93) into this last equation, we get

1 ρ (1 -v 2 )D ∂x ∂v = I -uu (x j -x i ) + O( D 3 ρ 2 ). (96)
By multiplying (95) with 1 ρ (1-v 2 )D, and by substituting (96), we get

1 ρ ∂q ∂v (x)(1 -v 2 )D = J q (u) I -uu (x j -x i ) + O( D 3 ρ 2 ).
Since function q(.) is 1-homogeneous, we note that

1 ρ ∂q ∂v (x) = ∂q ∂v (u) + O( D 2 ρ 2
), therefore we end up with

∂q ∂v (u)(1 -v 2 )D = • J q (u)(x j -x i ) + O( D 3 ρ 2 ),
where matrix

• J q (u) is defined in (48). Finally, this last equation and (93) yield

q (v, ϕ) vD + ∂q ∂v (v, ϕ) (1 -v 2 )D = q(u)u + • J q (u) (x j -x i ) + O( D 3 ρ 2 ). (97)
Since function β(.) is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong to L ∞ (R M ) as shown in Lemma 2, then (97) implies

β q (v, ϕ) vD + ∂q ∂v (v, ϕ) (1 -v 2 )D = β q(u)u + • J q (u) (x j -x i ) + O( D 2 ρ 2 ). (98)
APPENDIX C ASYMPTOTIC EXPANSION Proposition 8. Considering the stochastic reverberation model in Definition 1, suppose that the attenuation function α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] is regular (cf. Definition 7). Let ∀x ∈ R,

I(x) = u∈S 2 ψ(u)e -α(u)x du, ( 99
)
where ψ is a continuous function on S 2 . Then when x → +∞,

I(x) = e -α inf x x k∈K 2πψ(u k ) ¨ α k + o(1) . (100)
In addition, if ψ is differentiable almost everywhere on S 2 , and all its partial derivatives belong to L ∞ (S 2 ), then

I(x) = e -α inf x x k∈K 2πψ(u k ) ¨ α k + O 1 √ x . (101)
Finally, if ψ is twice (or more) continuously differentiable, then

I(x) = e -α inf x x k∈K 2πψ(u k ) ¨ α k + O 1 x . (102)
Proof. Since vectors {u k } k∈K are distinct, let {S 2 k } k∈K be a partition of S 2 , such that every u k is in the interior of S 2 k . Then I(x) = k∈K I k (x), where

I k (x) = u∈S 2 k ψ(u)e -α(u)x du. (103)
Since u k is the unique global minimum of the four times continuously differentiable function α(.) in the two-dimensional (2D) Riemannian manifold S 2 k , the following second order expansion holds:

α (u) = α inf +(u-u k ) Äk 2 (u-u k )+O( u-u k 3 2 ), (104)
where the 3 × 3 Hessian matrix Äk has rank 2 (with two positive eigenvalues and one zero eigenvalue), and span( Äk ) = span(u k ) ⊥ is the tangent space of the Riemannian manifold S 2 at point u k . By substituting (104) into (103), we get

I k (x)= u∈S 2 k ψ(u)e -α inf +(u-u k ) Äk 2 (u-u k )+O( u-u k 3 2 ) x du.
With the change of variable v = √ x(u -u k ), the surface element becomes dv = x du, thus

I k (x)= e -α inf x x v∈ √ x(S 2 k -u k ) ψ(u k + v √ x )e -v Äk 2 v+O v 3 2 √ x dv.
(105) If ψ is continuous, when x → +∞, since span(u k ) ⊥ = span( Äk ), we get the asymptotic form

I k (x) = e -α inf x x ψ(u k ) v∈span( Äk ) e -v Äk 2 v dv + o(1) = e -α inf x x ψ(u k ) 2π √ ¨ α k + o(1) , (106)
where ¨ α k > 0 is the Hessian of function α(.) at u k , i.e. the product of the two positive eigenvalues of the Hessian matrix Äk .

In addition, if ψ is differentiable almost everywhere on S 2 , and all its partial derivatives belong to L ∞ (S 2 ), then the next term in the Taylor series expansion of

ψ(u k + v √ x )e -v Äk 2 v+O v 3 2 √ x in (105) is of order 1 √
x . By integrating over span( Äk ), we get

I k (x) = e -α inf x x ψ(u k ) 2π ¨ α k + O( 1 √ x ) . (107)
Finally, if ψ is twice (or more) continuously differentiable, then the term of order 1 √ x in the Taylor series expansion of

ψ(u k + v √ x )e -v Äk 2 v+O v 3 2 √ x
in (105) is odd w.r.t. v, thus its integral over span( Äk ) is zero. The following term is of order 1

x , and since it is even w.r.t. v, its integral over span( Äk ) is not zero in general. Therefore we can write

I k (x) = e -α inf x x ψ(u k ) 2π ¨ α k + O( 1 x ) . (108)
By summing (106) (resp. ( 107) and ( 108)) for all k ∈ K, we finally get (100) (resp. ( 101) and ( 102)).

APPENDIX D PROPERTIES OF PROBABILITY DISTRIBUTION

p Lemma 5. Let p(b, z) (where b ∈ {0, 1} and z ∈ R M ) and p(b i , b j , z i , z j ) (where b i , b j ∈ {0, 1} and z i , z j ∈ R M ) be the probability distributions introduced in Definition 1. Then function (z 1 , z 2 ) → p(1, 1; z 1 , z 2 ) is not differentiable at any point (z, z) such that ∇p(1; z) = 0.
Proof. First, since V i (x, y; z i ) and V j (x, y; z j ) are Booleans, V i (x, y; z i )V j (x, y; z j ) ≤ min (V i (x, y; z i ), V j (x, y; z j )) .

By applying the mathematical expectation to both members of this inequality, we get

p(1, 1; z i , z j ) ≤ min (p(1; z i ), p(1; z j )) . (109)
In other respects, since function z → p(1; z) is continuously differentiable and it is not constant, then there is z ∈ R M such that h ∇p(1; z) = 0. Then ∀t > 0, since we have both [START_REF] Polack | Modifying chambers to play billiards: the foundations of reverberation theory[END_REF] and p(1, 1; z i , z j ) ≤ p(1; z j ), we get

p(1,1;z+th,z-th)-p(1,1;z,z) t ≤ p(1;z-th)-p(1;z) t . Therefore lim sup t→0,t>0 p(1,1;z+th,z-th)-p(1,1;z,z) t ≤ -h ∇p(1; z) = -h 2 2 .
In the same way, ∀t < 0, since we have both (9) and p(1, 1; z i , z j ) ≤ p(1; z i ), we get

p(1,1;z+th,z-th)-p(1,1;z,z) t ≥ p(1;z+th)-p(1;z) t . Therefore lim inf t→0,t<0 p(1,1;z+th,z-th)-p(1,1;z,z) t ≥ h ∇p(1; z) = h 2 2 .
Consequently, the limit of p(1,1;z+th,z-th)-p(1,1;z,z) t when t → 0 and t ∈ R does not exist, thus function

(z 1 , z 2 ) → p(1, 1; z 1 , z 2 ) is not differentiable at (z, z).
Nevertheless, it is still possible to assume that function (z 1 , z 2 ) → p(1, 1; z 1 , z 2 ) is continuous and differentiable almost everywhere in R M ×R M , and that all its partial derivatives belong to L ∞ (R M × R M ), as we did in Definition 1.

Proof of Lemma 1. We have assumed in Definition 1 that the support of z

→ p(1; z) is left-bounded, i.e. ∀m ∈ {1 . . . M }, ∃z inf m < 0 such that ∀z ∈ R M , if ∃m ∈ {1 . . . M } such that z m < z inf m , then p(1; z) = 0. Since ∀z ∈ R M , p(1; z) ∈ [0, 1], (13)
implies the majoration ∀n ∈ N\{0},

β n ≤ z≥z inf e -n z α c dz = c M M m=1 αm e n |z inf | α c n M = O e n |z inf | α c n M
, where symbol ≥ between two vectors is applied entrywise.

Moreover, β n > 0 in (13) because function z → p(1; z) is continuous, nonnegative, and not identically zero.

Proof of Lemma 2. First, function e → β(e) is even because function (z, z ) → p(1, 1; z, z ) is symmetric, since in Definition 1, p is invariant under a permutation of the two sensors i and j.

Moreover, we have assumed in Definition 1 that the support of

z → p(1; z) is left-bounded, i.e. ∀m ∈ {1 . . . M }, ∃z inf m < 0 such that ∀z ∈ R M , if ∃m ∈ {1 . . . M } such that z m < z inf m , then p(1; z) = 0. Consequently, (z, z ) → p(1, 1; z, z ) is also left-bounded: ∀z, z ∈ R M , if ∃m ∈ {1 . . . M } such that z m < z inf m or z m < z inf m , then p(1, 1; z, z ) = 0. Therefore ∀z ∈ R M , if ∃m ∈ {1 . . . M } such that z m < z inf m , then ∀e ∈ R M , p(1, 1; z -e 2 , z + e 2 ) = 0.
z → e -2z α c M m=1 1 [z inf m ,+∞[ (z m ).
Therefore the theorem of differentiability under the integral sign proves that function e → β(e) defined in ( 17) is differentiable almost everywhere in R M , and all its partial derivatives belong to L ∞ (R M ).

Besides, applying ( 17) to e = 0 and substituting (9) yields β(0) = β 2 , with β 2 defined in [START_REF] Badeau | General stochastic reverberation model[END_REF] for n = 2.

With the change of variables z = zi+zj 2

and e = z j -z i , which is such that dzde = dz i dz j , by multiplying both members of (109) with e -2z α c , integrating over z, and substituting [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept sine technique[END_REF] and [START_REF] Badeau | General stochastic reverberation model[END_REF] for n = 2, we get (18). Since z i , z j can take any value in R M , (18) holds ∀e ∈ R M . In particular, function e → β(e) is not differentiable at e = 0 (otherwise, with e = qcτ and τ ∈ R, we would have lim τ →0,τ >0 β(qcτ )-β(0) τ ≤ -β 2 α inf and lim τ →0,τ <0 β(qcτ )-β(0) τ ≥ +β 2 α inf , so the two limits cannot be equal). where function α(.) and β 1 are defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] and [START_REF] Badeau | General stochastic reverberation model[END_REF] for n = 1, and ∀u ∈ S 2 , ∀t ∈ R, ∀(x, y

) ∈ R 3 × R M , m s (u, t) = E [s (Θ(x, y)u, t)] .
Besides, E[h i (t)] = λβ 1 J ξ where J ξ is defined in (86) and

ξ(r) = r 2 g i r r 2 , . t * m s ( r r 2 , t -r 2 c ) e -α(r) c belongs to L 1 (R 3 ) ∩ L ∞ (R 3 ). Thus (92) yields E[h i (t)] = λβ 1 +∞ r=0 1 v=-1 2π ϕ=0 g i (v, ϕ, .) t * m s (v, ϕ, t -r c )e -α(v,ϕ) r c rdrdvdϕ. Let m g i * s (v, ϕ, .) = g i (v, ϕ, .) t * m s (v, ϕ, .) e α(v,ϕ)t = g i (v, ϕ, .) t * m s (v, ϕ, .
) and let T = T g + T s . With the change of variable τ = t -r c , which is such that dτ = dr c , we get ∀t ≥ T ,

E[h i (t)] = λc 2 β 1 1 v=-1 2π ϕ=0 e -α(v,ϕ)t t τ ∈R m g i * s (v, ϕ, τ ) dτ -τ ∈R m g i * s (v, ϕ, τ ) τ dτ dvdϕ = λc 2 β 1 1 v=-1 2π ϕ=0 e -α(v,ϕ)t t m g i * s (v, ϕ, 0) -ı 2π ∂ m g i * s ∂f (v, ϕ, 0) dvdϕ, (110)
where m g i * s (v, ϕ, f ) = g i (v, ϕ, f ) m s (v, ϕ, f ) with m s defined in [START_REF] Polack | La transmission de l'énergie sonore dans les salles (The transmission of sound energy in rooms)[END_REF]. Note that we have proved in Section II that function f → g i (u, f ) is smooth. Moreover, ∀Θ(x, y) ∈ SO(3), ∀u ∈ S 2 , ∀q ∈ N, function f → s (Θ(x, y)u, f ) is q times continuously differentiable and ∂ q s (Θ(x,y)u,f ) ∂f q is dominated by t q s (u, t) 1 . Since SO(3) is a compact set, the differentiability under the integral sign theorem proves that function f → m s (u, f ) is smooth. Therefore functions

m g i * s (v, ϕ, f ) and ∂ m g i * s ∂f (v, ϕ, f ) are well-defined.
Finally, both m g i * s (v, ϕ, 0) = 0 and ∂ m g i * s ∂f (v, ϕ, 0) = 0, since we have assumed in Definition 1 that at least one of the following properties holds ∀u ∈ S 2 at f = 0:

• ∀i ∈ {1 . . . I}, g i (u, 0) = 0 and m s (u, 0) = 0;

• ∀i ∈ {1 . . . I}, g i (u, 0) = 0 and

∂ g i ∂f (u, 0) = 0; • m s (u, 0) = 0 and ∂m s ∂f (u, 0) = 0. Therefore (110) proves that ∀t ≥ T , E[h i (t)] = 0.
Spectral domain: First, ∀Θ(x, y) ∈ SO(3), ∀f ∈ R, function u → s (Θ(x, y)u, f ) is twice continuously differentiable, and all its second order derivatives are bounded on S 2 . Since SO(3) is a compact set, the differentiability under the integral sign theorem proves that function u → m s (u, f ) defined in [START_REF] Vincent | Roomsimove[END_REF] is twice continuously differentiable.

Second, ∀Θ(x, y) ∈ SO(3), ∀u ∈ S 2 , ∀q ∈ N, function f → s (Θ(x, y)u, f ) is q times continuously differentiable and ∂ q s(Θ(x,y)u,f ) ∂f q is dominated by t q s (u, t) 1 . Since SO(3) is a compact set, the differentiability under the integral sign theorem proves that function f → m s (u, f ) is smooth.

By substituting (61) and ( 63) into (62), we get

E[ h i (f )] = λβ 1 r∈R 3 gi r r 2 ,f m s r r 2 ,f e -α(r)+2ıπf r 2 c r 2 dr
where function α(.) and β 1 are defined in ( 7) and ( 13) for n = 1. We note that E[ h i (f )] = λβ 1 J ξ where J ξ is defined in (86) and

ξ(r) = r 2 g i r r 2 , f m s r r 2 , f e -α(r)+2ıπf r 2 c belongs to L 1 (R 3 ) ∩ L ∞ (R 3 ). Thus (92) yields E[ h i (f )] = λβ 1 1 v=-1 2π ϕ=0 g i (v, ϕ, f ) m s (v, ϕ, f ) +∞ r=0 re -( α(v,ϕ)+2ıπf ) r c drdvdϕ. = λc 2 β 1 1 v=-1 2π ϕ=0 gi(v,ϕ,f )m s (v,ϕ,f ) ( α(v,ϕ)+2ıπf ) 2 dvdϕ,
which proves [START_REF] Stewart | Database of omnidirectional and B-format impulse responses[END_REF]. Finally, E[ h i (f )] is obtained as the integral over the compact set S 2 of a function of (u, f ) which is smooth w.r.t. f , and whose derivatives of all orders are bounded; therefore the differentiation under the integral sign theorem proves that function Proof of Proposition 2. Temporal domain: Substituting (60) and ( 69) into (68) yields

f → E[ h i (f )] is smooth.
γ i,j (t 1 , t 2 ) = λ x∈R 3 β(q(x -x i ) -q(x -x j )) e - α(x-x i )+ α(x-x j ) c gi x-x i x-x i 2 ,t1 gj x-x j x-x j 2 ,t2 x-xi 2 x-xj 2 t1 * t2 * m s,s x-xi x-xi 2 , t 1 -x-xi 2 c , x-xj x-xj 2 , t 2 - x-xj 2 c
dx (111) where functions α(.) and β(.) are defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] and [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept sine technique[END_REF], and

∀u 1 , u 2 ∈ S 2 , ∀t 1 , t 2 ∈ R, ∀(x, y) ∈ R 3 × R M , m s,s (u 1 , t 1 , u 2 , t 2 ) = E[s (Θ(x, y)u 1 , t 1 ) s (Θ(x, y)u 2 , t 2 )].
(112) Note that γ i,j (t 1 , t 2 ) = λ J ξ with J ξ defined in (85), and ξ(r 1 , r 2 ) = β(q(r 1 ) -q(r 2 ))e -α(r 1 )+ α(r 2 )

c g i r1 r1 2 , t 1 g j r2 r2 2 , t 2 t1 * t2 * m s,s r1 r1 2 , t 1 -r1 2 c , r2 r2 2 , t 2 -r2 2 c is such that x → ξ (x -x i , x -x j ) ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ).
Thus (91) yields

γ i,j (t 1 , t 2 ) = λ +∞ ρ= D 2 1 v=-1 2π ϕ=0 β q 2ρv+D 2ρ+vD , ϕ (ρ + vD 2 ) -q 2ρv-D 2ρ-vD , ϕ (ρ -vD 2 ) e - α ( 2ρv+D 2ρ+vD ,ϕ ) (ρ+ vD 2 )+ α ( 2ρv-D 2ρ-vD ,ϕ ) (ρ-vD 2 ) c g i 2ρv+D 2ρ+vD , ϕ, t 1 g j 2ρv-D 2ρ-vD , ϕ, t 2 t1 * t2 * m s,s 2ρv+D
2ρ+vD , ϕ, t 1 -

ρ+ vD 2 c , 2ρv-D 2ρ-vD , ϕ, t 2 - ρ-vD 2 c
dρdvdϕ.

(113) If the acoustic field is diffuse (cf. Definition 2) and if the microphones and the source are omnidirectional (cf. Definitions 3 and 4), we get the simplification 19) into (114), we get (24). Moreover, if D > 0, b(., D) ∈ L 1 (R) and has compact support, or if D = 0, b(., 0) = β 2 δ(.). In both cases, since functions g i and s are continuous with compact support, then function b(., D) * g i * g j * s * s is also continuous with compact support. Therefore function (t 1 , t 2 ) → γ i,j (t 1 , t 2 ) is continuous.

γ i,j (t 1 , t ) = 2πλ +∞ ρ= D 2 1 v=-1 β (qvD) e -2 α inf ρ c (g i * s) t 1 - ρ+ vD 2 c (g j * s) t 2 - ρ-vD
In particular, if i = j and t 1 = t 2 = t, we get (25), hence the expression of the temporal correlation in (26).

Spectral domain: Substituting (61) and (69) into (68) yields

cov[ h i (f 1 ), h j (f 2 )] = λ x∈R 3 β(q(x -x i ) -q(x -x j )) g i ( x-xi x-xi 2 , f 1 ) g j ( x-xj x-xj 2 , f 2 )m s, s ( x-xi x-xi 2 , f 1 , x-xj x-xj 2 , f 2 ) e - α(x-x i )+2ıπf 1 x-x i 2 + α(x-x j )-2ıπf 2 x-x j 2 c x-xi 2 x-xj 2 dx
where functions α(.) and β(.) are defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] and [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept sine technique[END_REF], and

∀u 1 , u 2 ∈ S 2 , ∀f 1 , f 2 ∈ R, ∀(x, y) ∈ R 3 × R M , m s, s (u 1 , f 1 , u 2 , f 2 ) = E[ s(Θ(x, y)u 1 , f 1 ) s(Θ(x, y)u 2 , f 2 )].
Note that cov[ h i (f 1 ), h j (f 2 )]=λJ ξ with J ξ defined in (85) and ξ(r 1 , r 2 ) = β(q(r 1 ) -q(r 2 ))

g i r1 r1 2 , f 1 g j r2 r2 2 , f 2 m s, s ( r1 r1 2 , f 1 , r2 r2 2 , f 2 )e -α(r 1 )+2ıπf 1 r 1 2 +α(r 2 )-2ıπf 2 r 2 2 c is such that x → ξ (x -x i , x -x j ) ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ). Thus (91) yields cov[ h i (f 1 ), h j (f 2 )] = λ +∞ ρ= D 2 1 v=-1 2π ϕ=0 β q 2ρv+D 2ρ+vD , ϕ (ρ + vD 2 ) -q 2ρv-D 2ρ-vD , ϕ (ρ -vD 2 ) e -( α ( 2ρv+D 2ρ+vD ,ϕ ) +2ıπf 1) (ρ+ vD 2 )+ ( α ( 2ρv-D 2ρ-vD ,ϕ ) -2ıπf 2 ) (ρ-vD 2 ) c g i 2ρv+D 2ρ+vD , ϕ, f 1 g j 2ρv-D 2ρ-vD , ϕ, f 2 m s, s 2ρv+D 2ρ+vD , ϕ, f 1 , 2ρv-D 2ρ-vD , ϕ, f 2 dρdvdϕ.
If the acoustic field is diffuse (cf. Definition 2) and if the microphones and the source are omnidirectional (cf. Definitions 3 and 4), we get the simplification

cov[ h i (f 1 ), h j (f 2 )] = 2πλ g i (f 1 ) g j (f 2 ) s(f 1 ) s(f 2 ) +∞ ρ= D 2 e -(2 α inf +2ıπ(f 1 -f 2 ))ρ c dρ 1 v=-1 β (qvD) e -2ıπ f 1 +f 2 2 vD c dv. (115)
With the change of variable τ = vD c , which is such that dτ = D c dv, by substituting (19) into (115), we get (27). Moreover, function (f 1 , f 2 ) → cov[ h i (f 1 ), h j (f 2 )] is smooth, as a product of smooth functions.

In particular, if i = j and f 1 = f 2 = f , we get (28), hence the expression of the spectral correlation in (29).

Time-frequency domain: The Wigner distribution of (111) is

W γi,j (t, f ) = λ x∈R 3 β(q(x -x i ) -q(x -x j )) W gi⊗gj x-xi x-xi 2 , x-xj x-xj 2 , ., f t * W ms,s ( x-xi x-xi 2 , x-xj x-xj 2 , t - x-xi 2+ x-xj 2 2c , f ) e - α(x-x i )+ α(x-x j )+2ıπf ( x-x i 2 -x-x j 2 ) c x-xi 2 x-xj 2 dx
where we have used the convolution property (4) of the Wigner distribution, and m s,s is defined in (112). Note that W γi,j (t, f ) = λ J ξ with J ξ defined in (85), and

ξ(r 1 , r 2 ) = β(q(r 1 ) -q(r 2 ))e -α(r 1 )+ α(r 2 )+2ıπf ( r 1 2 -r 2 2 ) c W gi⊗gj r1 r1 2 , r2 r2 2 , ., f t * W ms,s ( r1 r1 2 , r2 r2 2 , t -r1 2+ r2 2 2c , f ) is such that x → ξ (x -x i , x -x j ) ∈ L 1 (R 3 ) ∩ L ∞ (R 3 ). Thus (91) yields W γi,j (t, f ) = λ +∞ ρ= D 2 1 v=-1 2π ϕ=0 β q 2ρv+D 2ρ+vD , ϕ (ρ + vD 2 ) -q 2ρv-D 2ρ-vD , ϕ (ρ -vD 2 ) e - α ( 2ρv+D 2ρ+vD ,ϕ ) (ρ+ vD 2 )+ α ( 2ρv-D 2ρ-vD ,ϕ ) (ρ-vD 2 )+2ıπf vD c W gi⊗gj 2ρv+D 2ρ+vD , ϕ, 2ρv-D 2ρ-vD , ϕ, ., f t * W ms,s ( 2ρv+D 2ρ+vD , ϕ, 2ρv-D 2ρ-vD , ϕ, t -ρ c , f )dρdvdϕ. ( 116
∀t ≥ T + D 2c , W γi,j (t, f ) = 2πλc 2 D e -2 α inf t D c τ =-D c β (qcτ ) e -2ıπf τ dτ t ∈R e 2 α inf t W gi⊗gj (., f ) t * W ms,s (t , f )dt . ( 117
) Finally, by using the projection property (2) of the Wigner distribution and by substituting (19) into (117), we get (30). Moreover, function f → B i,j (f, D) defined in (31) is smooth, as a product of smooth functions.

In particular, if i = j, we get (32), hence the expression of the time-frequency correlation in (33). Then (34) is obtained as the inverse Fourier transform of the right member of (33) when ∀f ∈ R, ∠ g i (f ) = ∠ g j (f ).

Finally, function τ → σ i,j (τ, x j -x i ) in (34) has the same support as function τ → b(τ, D) defined in (19). It is continuous in the interior of this support, it reaches its maximum at τ = 0, and it is not differentiable at τ = 0, because so is function τ → b(τ, D) as shown in Lemma 3. Moreover, (20) yields σ i,j (τ,

x j -x i ) ∈ [0, c 2D ].

APPENDIX G ASYMPTOTIC RESULTS

In this Appendix, we first prove Propositions 3 and 4 in Section III. Then Corollaries 1 and 2 will be proved in Appendix G.5.

G.1. General asymptotic results

First, we compute a few simple asymptotic forms, that will be used in Appendices G.2 and G.3.

Suppose that ρ → +∞. Then we get

2ρv+D 2ρ+vD = v + (1-v 2 )D 2ρ+vD = v + O( D ρ ), 2ρv-D 2ρ-vD = v -(1-v 2 )D 2ρ-vD = v + O( D ρ ), (118)
and since function α(.) is twice continuously differentiable,

α 2ρv+D 2ρ+vD , ϕ (ρ + vD 2 ) + α 2ρv-D 2ρ-vD , ϕ (ρ -vD 2 ) = 2 α (v, ϕ) ρ + O D 2 sup S 2 | ∂ 2 α ∂v 2 | ρ , therefore e - α 2ρv+D 2ρ+vD ,ϕ (ρ+ vD 2 )+ α 2ρv-D 2ρ-vD ,ϕ (ρ- vD 2 ) c = e -2 α(v,ϕ)ρ c 1 + O D 2 sup S 2 | ∂ 2 α ∂v 2 | ρc . (119)
In the same way, since function q(.) is twice continuously differentiable,

q 2ρv+D 2ρ+vD , ϕ (ρ + vD 2 ) -q 2ρv-D 2ρ-vD , ϕ (ρ -vD 2 ) = q (v, ϕ) vD + ∂q(v,ϕ) ∂v (1 -v 2 )D + O D 2 sup S 2 | ∂ 2 q ∂v 2 | ρ .
(120) Since function β(.) is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong to L ∞ (R M ) as shown in Lemma 2, then (120) yields

β q 2ρv+D 2ρ+vD , ϕ (ρ + vD 2 ) -q 2ρv-D 2ρ-vD , ϕ (ρ -vD 2 ) = β q (v, ϕ) vD+ ∂q(v,ϕ) ∂v (1-v 2 )D +O D sup S 2 | ∂ 2 q ∂v 2 | ρ . (121)

G.2. Temporal domain

Since functions g i , g j and s have finite temporal support, when t 1 +t 2 → +∞ with t 1 -t 2 fixed, (113) can be simplified by noting that ρ → +∞. Indeed, substituting the asymptotic forms (118), ( 119) and ( 121) into (113) yields:

γ i,j (t 1 , t 2 ) = λ +∞ ρ= D 2 1 v=-1 2π ϕ=0 e -2 α(v,ϕ)ρ c β q (v, ϕ) vD + ∂q(v,ϕ) ∂v (1 -v 2 )D g i (v, ϕ, t 1 ) g j (v, ϕ, t 2 ) t1 * t2 * m s,s v, ϕ, t 1 - ρ+ vD 2 c , v, ϕ, t 2 - ρ-vD 2 c +O D ρ dvdϕdρ, (122)
where m s,s is defined in (112).

By substituting (93) and ( 98) into (122), we get

γ i,j (t 1 , t 2 ) = λ +∞ ρ= D 2 u∈S 2 e -2 α(u)ρ c β q(u)u + • J q (u) (x j -x i ) g i (u, t 1 ) g j (u, t 2 ) t1 * t2 * m s,s u, t 1 - ρ+ u (x j -x i ) 2 c , u, t 2 - ρ- u (x j -x i ) 2 c +O D ρ dudρ. (123) By applying the change of variable τ = t1+t2 2 -ρ c , which is such that dτ = dρ c , to (123), we get ∀t 1 + t 2 ≥ 2T + D c , γ i,j (t 1 , t 2 ) = λc u∈S 2 e -α(u)(t1+t2) β q(u)u + • J q (u) (x j -x i ) g i (u, .) * g j (u, .) * m s * s (u, t 1 -t 2 - u (xj -xi) c ) + O D c(t1+t2) du, (124)
where m s * s is defined in (35).

Note that ∀Θ(x, y) ∈ SO(3), ∀t ∈ R, function u → s (Θ(x, y)u, .) t * s (Θ(x, y)u, .) is twice continuously differentiable, and all its second order derivatives are bounded on S 2 . Since SO(3) is a compact set, the differentiability under the integral sign theorem proves that function m s * s is twice continuously differentiable w.r.t. u ∈ S 2 . In the same way, ∀Θ(x, y) ∈ SO(3), ∀u ∈ S 2 , function t → s (Θ(x, y)u, .) t * s (Θ(x, y)u, .) is continuous and dominated by T s sup u∈S 2 s (u, .) 2 ∞ . Since SO(3) is a compact set, the continuity under the integral sign theorem proves that function m s * s is continuous w.r.t. t ∈ R.

1) Early asymptotic state: Suppose that ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant vector, therefore function α(.) defined in ( 7) is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α. Then (124) becomes

∀t 1 + t 2 ≥ 2T + D c , γ i,j (t 1 , t 2 ) = λc e -α inf (t1+t2) u∈S 2 β qu (x j -x i ) g i (u, .) * g j (u, .) * m s * s (u, t 1 -t 2 - u (xj -xi) c )du + O( D c(t1+t2) ) .
In particular, if i = j (thus D=0) and t 1 = t 2 = t, we get (36).

2) Late asymptotic state: Suppose that the attenuation function is regular (cf. Definition 7). If we let x = t 1 + t 2 , then (124) shows that γ i,j (t 1 , t 2 )) can be decomposed as

γ i,j (t 1 , t 2 )) = I 1 (x) + I 2 (x) O D c(t1+t2) (125)
where I 2 (x) is the integral defined in (99) with ψ(u) = λc, so that (102) proves that

I 2 (x) = O e -α inf (t 1 +t 2 ) t1+t2 , (126)
and I 1 (x) is the integral defined in (99) with

ψ(u) = λcβ q (u) u T + • J q (u) r g i (u, .) * g j (u, .) * m s * s (u, t 1 -t 2 -u r c ) (127
) where r = x j -x i . Since function β(.) is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong to L ∞ , so is function ψ in (127). If in addition i = j, then r = 0 and function ψ is twice continuously differentiable.

Therefore if r = 0, substituting (101), ( 126) and (127) into (125) shows that when t 1 + t 2 → +∞,

γ i,j (t 1 , t 2 ) = 2πλc e -α inf (t 1 +t 2 ) t1+t2 k∈K β k (r) √ ¨ α k g i (u k , .) * g j (u k , .) * m s * s (u k , t 1 -t 2 -τ k ) + O 1 √ t1+t2
where τ k = u k r c and β k (r) is defined in (51). If i = j (thus r = 0 and τ k = 0) and t 1 = t 2 = t, substituting (102) and ( 126) into (125) proves (49).

G.3. Time-frequency domain

Because of the temporal support property of the Wigner distribution, W gi⊗gj and W ms,s have finite temporal support. Therefore when t → +∞, (116) can be simplified by noting that ρ → +∞. Indeed, by substituting the asymptotic forms (118), ( 119) and ( 121) into (116), we get:

W γi,j (t, f ) = λ +∞ ρ= D 2 1 v=-1 2π ϕ=0 e -2 α(v,ϕ)ρ c β q (v, ϕ) vD + ∂q(v,ϕ) ∂v (1 -v 2 )D e -2ıπf vD c W gi⊗gj (v, ϕ, ., f ) t * W ms,s (v, ϕ, t-ρ c , f ) + O D ρ dvdϕdρ, (128)
where m s,s is defined in (112).

By substituting (93) and ( 98) into (128), we get

W γi,j (t, f ) = λ +∞ ρ= D 2 u∈S 2 e -2 α(u)ρ c β q(u)u + • J q (u) (x j -x i ) e -2ıπf u (x j -x i ) c W gi⊗gj (u, ., f ) t * W ms,s (u, t -ρ c , f ) + O D ρ dudρ.
(129) By applying the change of variable τ = t-ρ c , which is such that dτ = dρ c , to (129), we get ∀t ≥ T + D 2c with T = T g + T s ,

W γi,j (t, f ) = λc u∈S 2 e -2 α(u)t β q(u)u + • J q (u) (x j -x i ) e -2ıπf u (x j -x i ) c g i (u, f ) g j (u, f )m | s | 2 (u, f ) + O D ct du, (130)
where we have used the projection property (2) of the Wigner distribution, and the definition of

m | s | 2 (u, f ) in (37).
Finally, note that ∀Θ(x, y) ∈ SO(3), ∀f ∈ R, function u → s (Θ(x, y)u, f ) 2 is twice continuously differentiable, and all its second order derivatives are bounded on S 2 . Since SO(3) is a compact set, the differentiability under the integral sign theorem proves that function u → m s (u, f )

2 is twice continuously differentiable. Also note that ∀Θ(x, y) ∈ SO(3),

∀u ∈ S 2 , ∀q ∈ N, function f → s (Θ(x, y)u, f ) 2 is q times continuously differentiable and ∂ q | s (Θ(x,y)u,f )| 2
∂f q is dominated by t q s (u, t) 2 1 . Since SO(3) is a compact set, the differentiability under the integral sign theorem proves that function f → m s (u, f ) 2 is smooth.

1) Early asymptotic state: Suppose that ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant vector, therefore function α(.) defined in ( 7) is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α; Then (130) yields (38) and (39). Moreover, B i,j (f, r) is obtained as the integral over the compact set S 2 of a function of (u, f ) which is smooth w.r.t. f , and whose derivatives of all orders are bounded; therefore the differentiation under the integral sign theorem proves that function f → B i,j (f, r) is smooth. The proof of the even symmetry of f → B i,j (f, r) is straightforward. If i = j (thus D = 0), then (130) implies (40), and (38) leads to the expression of the time-frequency correlation in (41). This time-frequency correlation is asymptotically smooth and even symmetric w.r.t. f because so is function f → B i,j (f, x j -x i ).

2) Late asymptotic state: Suppose that the attenuation function is regular (cf. Definition 7). If we let x = 2t, then (130) shows that W γi,j (t, f ) can be decomposed as W γi,j (t, f ) = I 1 (x) + I 2 (x) O D ct (131) where I 2 (x) is the integral defined in (99) with ψ(u) = λc, so that (102) proves that

I 2 (x) = O e -2 α inf t t , (132)
and I 1 (x) is the integral defined in (99) with ψ(u) = λcβ q (u) u T +

• J q (u) r

g i (u, f ) g j (u, f )m | s | 2 (u, f )e -2ıπf u r c ( 133
)
where r = x j -x i . Since function β(.) is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong to L ∞ , so is function ψ in (133). If in addition i = j, then r = 0 and function ψ is twice continuously differentiable. Therefore if r = 0, substituting (101) and ( 132) into (131) proves (52). Moreover, function f → B i,j (f, r) in ( 50) is smooth, as a finite sum of smooth functions. The proof of its even symmetry is straightforward.

If i = j (thus r = 0), substituting (102), ( 132) and ( 133) into (131) proves (53), and (52) leads to the expression of the time-frequency correlation in (54). This time-frequency correlation is asymptotically smooth and even symmetric w.r.t. f because so is function f → B i,j (f, x j -x i ).

G.4. Asymptotic normality

Substituting (60) into (80) applied to θ 1 = θ and θ 2 = 0 yields ln φ hi(t) (θ) = λ x∈R 3 y∈R M p(1; yq(xx i )) where ∀t ∈ R, κ n (t) is the n-th order cumulant of h i (t), and ∀u ∈ S 2 , µ n (u, .) = E g i (u, .) t * s (Θ(x, y)u, .) n is the n-th order moment of g i (u, .) t * s (Θ(x, y)u, .). Note that ∀t ∈ R, ∀Θ(x, y) ∈ SO(3), function u → g i (u, .) The changes of variables r = x-x i and z = y -q(r), which are such that dr = dx and dz = dy, yield where function α(.) is defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] and ∀n ∈ N\{0}, β n is defined in [START_REF] Badeau | General stochastic reverberation model[END_REF]. Note that (137) is of the form J ξ in (86) where 1) Early asymptotic state: Suppose that ∀r ∈ R 3 , q(r) = q r 2 where q ∈ R M + is a constant vector, therefore function α(.) defined in [START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF] is constant on S 2 : ∀u ∈ S 2 , α(u) = α inf , where α inf = q α. Then (138) becomes converges pointwise to that of the standard normal distribution when t → +∞, which proves that it is asymptotically normally distributed.

φ
2) Late asymptotic state: Suppose that the acoustic field is non-diffuse and the attenuation function is regular (cf. Definition 7). Then (138) shows that κ n (t) can be written as the integral I(x) in (99), if we let x = nt and

ψ(u) = λβ n c n-3 T τ =0
e n α(u)τ µn(u,τ ) (t-τ ) n-2 dτ.

Note that ∀t > T , ∀τ ∈ [0, T ], function u → e n α(u)τ µn(u,τ ) (t-τ ) n-2 is continuous since we have proved that function u → µ n (u, τ ) is continuous, and it is dominated by e n ( sup S 2 α(.) ) T C n (t-T ) n-2

. Since [0, T ] is a compact set, the theorem of continuity under the integral sign proves that function u → ψ(u) is continuous on R 3 . Since ψ is continuous, (100) proves that when t → +∞, κ n (t) ∼ 2πλβn , we obtain Therefore the characteristic function of (143) Applying the inverse Fourier transform to the limit when t → +∞ of the right member of (143) yields

which shows that function τ → σ i,j (τ, x j -x i ) is obtained as the product of a continuous function of τ (as shown in Lemma 2) and of the integral on the compact set [0, 2π] of a function of (τ, ϕ) which is bounded by ξ i,j ∞ and continuous w.r.t. τ in ] -D c , D c [, therefore the continuity under the integral sign theorem proves that function τ → σ i,j (τ, x j -x i ) is continuous in ] -D c , D c [. If function ξ i,j (.) is even on S 2 (symmetric case), then (145) shows that function τ → σ i,j (τ, x j -x i ) is also even, thus its Fourier transform lim t→+∞ ρ i,j (t, f, x j -x i ) is real-valued.

If moreover function ξ i,j (.) is constant on S

Finally, applying the change of variable τ = vD c , which is such that dτ = D c dv, and substituting (19) into (146) yields (44). Then (45) is obtained as the inverse Fourier transform of the limit when t → +∞ of the right member of (44). Moreover, (20) yields σ i,j (τ, x j -x i ) ∈ [0, c 2D ]. Besides, function τ → σ i,j (τ, x j -x i ) reaches its maximum at τ = 0, and it is not differentiable at τ = 0, because so is function τ → b(τ, D) as shown in Lemma 3. Finally, (46) is derived by substituting (23) into (44), and (47) is derived by substituting (22) into (45).

Proof of Corollary 2. Substituting (42) into (53) yields (55). Then substituting (50) and (42) into (54) yields (56). Finally, (57) is obtained as the inverse Fourier transform of the limit when t → +∞ of the right member of (56).

  • I: identity matrix; • span(A): column (or range) space of matrix A; • [a, b]: closed interval, including a and b ∈ R; • ]a, b[: open interval, excluding a and b ∈ R;

Fig. 1 .

 1 Fig. 1. Temporal power profile var[h i (t)] in the shoebox room

Fig. 1 -

 1 Fig.1-(a) represents the temporal power profile var[h i (t)] in dB (blue curve), superimposed with a straight red line obtained by linear regression. We can observe that the temporal power profile is slightly bent compared with the straight line, so we are neither in a diffuse acoustic field, nor in the early asymptotic state of a non-diffuse acoustic field in TableI

Fig. 3 .

 3 Fig. 3. Real part of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in the shoebox room

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Last time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in the shoebox room (late asymptotic state)

Fig. 8 .

 8 Fig. 8. Temporal power profile var[h i (t)] in the classroom

Fig. 10 .

 10 Fig. 10. Real part of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in the classroom

  ) represents the cardinal sine function sinc 2πf D c

Fig. 12 .

 12 Fig. 12. Average over time of the time-frequency correlation ρ 1,2 (t, f, x 2x 1 ) in the classroom

Fig. 13 .

 13 Fig. 13. Normalized asymptotic correlation function 2D c σ 1,2 (τ, x 2 -x 1 ) in the classroom

Fig. 14 .

 14 Fig. 14. Temporal power profile var[h i (t)] in the Octagon room

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Real part of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in the Octagon room

Fig. 18 .

 18 Fig. 18. Normalized asymptotic correlation function 2D c σ 1,2 (τ, x 2 -x 1 ) in the Octagon room

Fig. 19 .

 19 Fig. 19. Temporal power profile var[h i (t)] in the Great Hall

Fig. 21 .

 21 Fig. 21. Real part of the time-frequency correlation ρ 1,2 (t, f, x 2 -x 1 ) in the Great Hall

Fig. 23 .

 23 Fig. 23. Normalized asymptotic correlation function 2D c σ 1,2 (τ, x 2 -x 1 ) in the Great Hall

  where ∀(x, y) ∈ R 3 × R M , ∀r ∈ R 3 , |ϕ i (x, y; r)| ≤ min(T g , T s ) sup u∈S 2 ,t∈[0,Tg] |g i (u, t) | sup u∈S 2 ,t∈[0,Ts] |s (u, t) |, and (

x∈R 3 2 Vj

 32 y∈R M ψ(x, y)dN (x, y), and ψ(x, y) = Re θ 1 Vi(x,y;y-q(x-xi))ϕi(x,y;x-xi)x-xi 2 +θ (x,y;y-q(x-xj ))ϕj (x,y;x-xj )x-xj 2 e -y α c .

  x∈R 3 y∈R M (e ıψ(x,y) -1)dxdy n = e λ x∈R 3 y∈R M (E[e ıψ(x,y) ]-1)dxdy .

  Finally, function (e, z) → p(1, 1; ze 2 , z + e 2 )e -2z α c is continuous w.r.t. e, and it is dominated by the integrable function z → e -2z α c M m=1 1 [z inf m ,+∞[ (z m ). Therefore the theorem of continuity under the integral sign proves that function e → β(e) defined in (17) is continuous on R M . In the same way, function (e, z) → p(1, 1; ze 2 , z + e 2 )e -2z α c is differentiable w.r.t. e almost everywhere in R M , and all its partial derivatives w.r.t. e are dominated by a constant multiplied by the integrable function

  FIRST ORDER MOMENTSProof of Proposition 1. Temporal domain: Substituting (60) and (63) into (62) yieldsE[h i (t)] = λβ 1 r∈R 3

  Proof of Lemma 3. First, function τ → b(τ, D) defined in (19) is even because Lemma 2 shows that function β(.) is even. As a consequence, function f → b(f, D) is real-valued.Second, function τ → b(τ, D) is continuous and differentiable almost everywhere in the interior of its support, and∂b ∂τ (τ, D) ∈ L ∞ (] -D c , + D c [),because Lemma 2 shows that function β(.) is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong toL ∞ (R M ). As a consequence, function f → b(f, D) is smooth.Besides, applying (18) in Lemma 2 to e = qcτ leads to β (qcτ ) ≤ β 2 e -α inf |τ | , which proves (20). In particular, function τ → b(τ, D) reaches its maximum at τ = 0, and it is not differentiable at τ = 0 (otherwise we would havelim τ →0,τ >0 b(τ,D)-b(0,D) τ = ∂b ∂τ (0, D) ≤ -β 2 α inf and lim τ →0,τ <0 b(τ,D)-b(0,D) τ = ∂b ∂τ (0, D) ≥ +β 2 α inf ). Moreover, ∀f ∈ R, | b(f, D)| ≤ τ∈R b(τ, D)dτ , thus (20) implies (21). If D → 0, since function β(.) is continuous and differentiable almost everywhere in R M , and all its partial derivatives belong to L ∞ (R M ), we have β (qcτ ) = β 2 +O(τ ) ∀τ ∈ [-D c , D c ]. Therefore (19) yields (22), whose Fourier transform leads to (23).

  the changes of variables t = t1+t2 2 -ρ c and τ = vD c , which are such that dt = dρ c and dτ = D c dv, and by substituting (

  ) If the acoustic field is diffuse (cf. Definition 2) and if the microphones and the source are omnidirectional (cf. Definitions 3 and 4), with the changes of variables t = t -ρ c and τ = vD c , which are such that dt = dρ c and dτ = D c dv, (116) yields

2 , 2 c- 1

 221 x-xi x-xi 2 , t -x-xi dxdy.(134) We consider the following series expansions:ln φ hi(t) (θ) = +∞ n=1 ı n n! θ n κ n (t) φ gi * s (θ, u, t) = +∞ n=0 ı n n! θ n µ n (u, t)(135)

t

  * s (Θ(x, y)u, .) n is continuous, and it is dominated by C n , with C = min(T g , T s ) sup u∈S 2 ,t∈[0,Tg]

  |g i (u, t)| sup u∈S 2 ,t∈[0,Ts] |s(u, t)|.Since SO(3) is a compact set, the theorem of continuity under the integral sign proves that function u → µ n (u, t) is continuous on R 3 . Also note thatsup u∈S 2 ,t∈R |µ n (u, t)| ≤ C n . (136)By substituting (135) into (134) and by identifying the n-th order terms, we get ∀n ∈ N\{0},κ n (t) = λ x∈R 3 y∈R M p(1; yq(xx i )) 2 , t -x-xi 2 cdxdy.

2 , t -r 2 c 3 1 v=- 1 2π

 2231 belongs to L 1 (R 3 ) ∩ L ∞ (R 3 ), because its support does not contain the point r = 0, |µ n | is upper bounded and α(r) ≥ α inf r 2 . Therefore (92) yieldsκ n (t) = λβ n The change of variable τ = t-r c , which is such that dτ = dr c , yields ∀t ≥ T = T g + T s , κ n (t) = λβn c n-ϕ=0 e -n α(v,ϕ)t T τ =0 e n α(v,ϕ)τ µn(v,ϕ,τ ) (t-τ ) n-2 dτ dvdϕ(138)because the support of µ n (v, ϕ, .) lies in [0, T ].

=0 e 2 α 2 . ( 141 )n

 22141 κ n (t) = 4πλβ n c n-3 e -n α inf t T τ =0 e n α inf τ µn(τ ) (t-τ ) n-2 dτ, v, ϕ, τ ) dvdϕ.We note that (136) implies∀τ ∈ R, |µ n (τ )| ≤ C n .(140)In particular, for n = 2, (139) impliesκ 2 (t) = 4πλcβ 2 e -2 α inf t T τ inf τ µ 2 (τ ) dτ.Therefore, ∀n ≥ 3, when t → +∞, =0 e n α inf τ |µn(τ )|dτ c n-3 t n-By substituting (14) and (140) into (141), we thus get ∀n ≥ 3,κ n (t) (κ 2 (t)) |z inf | α c + α inf T n M +1 c n-3 t n-3

nc n- 3 e.τ

 3 -n α inf t t n-1 k∈K T τ =0 e n α(u k )τ µn(u k ,τ )dτ √ ¨ α kIn particular, for n = 2, we getκ 2 (t) ∼ πλcβ 2 e -2 α inf t t k∈K T τ =0 e 2 α(u k )τ µ2(u k ,τ )dτ √ ¨ α k .Therefore, ∀n ≥ 3, when t → +∞,κ n (t) (κ 2 (t)) =0 e n α(u k )τ |µn (u k ,τ )|dτ√ substituting (14) and (136) into (142), we thus get ∀n ≥ 3,

ρ

  i,j (t, f, r) = u∈S 2 β(qu r)ξi,j (u)e -2ıπf u r c du β2 √ u∈S 2 ξi,i(u)du √ u∈S 2 ξj,j (u)du + O( 1 t ).

  σ i,j (τ, r) = u∈S 2 β(qu r)ξi,j (u)δ τ -u r c du β2 √ u∈S 2 ξi,i(u)du √ u∈S 2 ξj,j (u)du ,(144)which proves that τ → σ i,j (τ, r) is nonnegative and that its support is[-D c , D c ]. Moreover, with the change of variables u = √ 1 -v 2 cos(ϕ), √ 1 -v 2 sin(ϕ), v with v ∈ [-1, 1] and ϕ ∈ [0, 2π], (144) can be rewritten asσ i,j (τ, r) = 1 v=-1 2π ϕ=0 β(qvD)ξi,j (v,ϕ)δ(τ -vD c )dvdϕ β2 √ u∈S 2 ξi,i(u)du √ u∈S 2 ξj,j (u)du . Therefore ∀τ ∈ [-D c , D c ], σ i,j (τ, r) = c D β(qcτ ) β2 2π ϕ=0 ξi,j ( τ c D ,ϕ)dϕ√ u∈S 2 ξi,i(u)du √ u∈S 2 ξj,j (u)du ,

1 -

 1 v 2 cos(ϕ), √ 1 -v 2 sin(ϕ), v with v ∈ [-1, 1] and ϕ ∈ [0, 2π] in (143) yields

  ρ i,j (t, f, r) = 1 v=-1 β(qvD)e -2ıπf vD c dv 2β2 + O( 1 t ).

  t) is continuous w.r.t. t ∈ R and twice continuously differentiable w.r.t. u ∈ S 2 . Moreover, ∀u ∈ S 2 , the temporal support of function t → g i (u, t) is included in [0, T g ] with T g > 0; • function s (u, t) is continuous w.r.t. t ∈ R and twice continuously differentiable w.r.t. u ∈ S 2 . Moreover, ∀u ∈ S 2 , the temporal support of function t → s (u, t) is included [0, T

s ] with T s > 0;

• at least one of the three following properties holds ∀u ∈ S 2 at f = 0:

  it is continuous w.r.t. τ and differentiable almost everywhere in the interior of its support, and ∂b ∂τ

  in the isotropic case, if the microphones are omnidirectional and have the same response g, we get ξ i,j (u) = 1

	Convergence	Diffuse acoustic field	Non-diffuse acoustic field, early asymptotic state	Non-diffuse acoustic field,
	speed to	(Proposition 2)	(Proposition 3 & Corollary 1)		late asymptotic state
	asymptotic	Omnidirectional source	Anisotropic	Symmetric	Isotropic	(Proposition 4 & Corollary 2)
	state	Same omnidirectional sensors	(general case)	(ξ i,j even)	(ξ i,j constant)

TABLE III SENSOR

 III . ORIENTATIONS (AZIMUTH, ELEVATION AND ROLL OFFSET IN DEGREES, POSITIVE FOR SLEW LEFT, NOSE UP OR RIGHT WING DOWN)

	Sensor	Azimuth	Elevation	Roll offset
	Sensor 1	15 •	25 •	35 •
	Sensor 2	45 •	55 •	65 •

  .5. Proofs of Corollaries 1 and 2 in Section III Proof of Corollary 1. By substituting (42) into (40), we get (43). Then substituting (39) and (42) into (41) yields

hi(t)

√ var[hi(t)] converges pointwise to that of the standard normal distribution when t → +∞, which proves that it is asymptotically normally distributed.

G

Uniform means invariant under any translation, and isotropic means invariant under any rotation in the three-dimensional (3D) space.

This remark is very important, because the fact that this function is not twice continuously differentiable is the reason for the slower speed of convergence O( 1 √ t) in late asymptotic state (see the discussion in Section III-E).

The Matlab code generating all the figures in Section IV is available at https://perso.telecom-paristech.fr/rbadeau/techreport2019-04-code.zip.

Note that the four signatures of reverberation are based on the second order moments of the RIR. In the following experiments, we will display neither the expected values nor the higher order cumulants of h i (t): the zero mean and asymptotic Gaussianity of RIRs (which are mathematically proved in Propositions 1, 3 and 4) are well-known experimental facts.

In Section IV-B, in order to denoise the estimate of σ 1,2 (τ, x 2 -x 1 ), (58) will be truncated to the frequency band [-5000 Hz, +5000 Hz] before computing the inverse DFT, which is equivalent to smoothing in time domain.

The modeling of absorption due to the air involves a frequency-varying attenuation, which is not accounted for by the model presented in this research report. However, frequency-varying attenuations are handled by the general stochastic reverberation model introduced in[START_REF] Badeau | General stochastic reverberation model[END_REF], and will be analyzed both mathematically and experimentally in future work.

Equation (B8) was proved in[START_REF]Common mathematical framework for stochastic reverberation models[END_REF] in the real case; we use here its extension to the complex case, which is straightforward.

Télécom ParisTech

Institut Mines-Télécom -membre de l'Université Paris Saclay 46, rue Barrault -75634 Paris Cedex 13 -Tél. + 33 (0)1 45 81 77 77 -www.telecom-paristech.fr

Département IDS

© Institut Mines-Télécom -Télécom ParisTech 2019