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

General stochastic reverberation model Modèle stochastique général de réverbération

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

In a recent research paper, we proposed 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, spectral and temporal domains. This model was dedicated to diffuse (i.e. isotropic and uniform) acoustic fields, omnidirectional sources and microphones, and constant attenuation coefficients w.r.t. the frequency. In this technical report, we introduce several extensions of this model, that aim to model reverberation more realistically, by considering anisotropic and non-uniform acoustic fields, directive sources and microphones, and frequency-varying attenuation coefficients.

I. INTRODUCTION

In [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF], we proposed 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, spectral and temporal 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 stochastic model 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 located outside the room. In [START_REF] Badeau | 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. More precisely, given a Borel set V ⊂ R 3 of finite volume |V |, we assume that the number N (V ) of source images contained in V follows a Poisson distribution of rate parameter λ|V |: N (V ) ∼ P(λ|V |).

Mathematically, this is formalized through the concept of Poisson random measures with independent increments: 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).

In the stochastic reverberation model proposed in [START_REF] Badeau | 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.

In other respects, we assumed that the source images emit spherical acoustic waves, which undergo an exponential attenuation along their trajectories, that is due to the multiple reflections on the room surfaces and to propagation in the air. We considered a diffuse (i.e. uniform and isotropic) acoustic field, so we further assumed that this attenuation is isotropic. We also assumed that the attenuation is constant w.r.t. the frequency, so that it only depends on the length of the trajectory. Finally, we considered omnidirectional sources and microphones. These assumptions resulted in the following stochastic model: for any sensor i ∈ {1 . . . I}, time t ∈ R, and frequency f ∈ R,

h i (t) = x∈R 3 (gi * s) t- x-x i 2 c e -α x-x i 2 c x-xi 2 dN (x) ĥi (f ) = ĝi (f ) ŝ (f ) x∈R 3 e -(α+2ıπf ) x-x i 2 c x-xi 2 dN (x) (1) where
• h i (t) ∈ R (resp. ĥi (f ) ∈ C) is the room impulse response (resp. room frequency response) at sensor i;

• g i (t) ∈ R (resp. ĝi (f ) ∈ C) is the impulse response (resp.
frequency response) of the omnidirectional sensor i; • s (t) ∈ R (resp. ŝ (f ) ∈ C) is the impulse response (resp. frequency response) of the omnidirectional source; • vector x i ∈ R 3 (in meters) is the position of sensor i;

• vector x ∈ R 3 (in meters) represents the positions of the source images. As mentioned above, it is distributed according to a uniform and isotropic Poisson point process dN (x) ∼ P(λdx) with λ > 0 (in meters -3 ); • c > 0 is the speed of sound (in meters.hertz);

• α > 0 is the constant attenuation coefficient (in hertz);

1 • 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 . Note that (1) is a re-parametrization from [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF], where g i (t) and s(t) were merged in a single function g(t) = (g i * s)(t + T ) e αt , with T > 0.

This model was thus limited to diffuse acoustic fields, omnidirectional sources and microphones, and constant attenuation coefficients w.r.t. the frequency. In the following sections, we will propose successive extensions of this model, that aim to model reverberation more realistically, by considering anisotropic and non-uniform acoustic fields, directive sources and microphones, and frequency-varying attenuation coefficients. For this purpose, we will make several experimental observations 1 regarding some statistical properties of the source images (Sections II to V), that we do not aim to prove mathematically (this could be the topic of other papers), but rather to describe mathematically with our new stochastic reverberation model (that will be mathematically formulated in Section VI). In future papers, we will show that this model is able to predict various statistical properties of reverberation (regarding the first and second order moments in the spatial, spectral, temporal and time-frequency domains, and the asymptotic Gaussianity), that we experimentally observed in a wide variety of both synthetic and measured room impulse responses. In the following sections, we will use the same mathematical notation as in [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF].

II. DIFFUSE REVERBERATION MODEL FOR OMNIDIRECTIONAL SOURCES AND SENSORS

A. Simulation results: ergodic room

We consider a bi-dimensional room 2 whose geometry is depicted in Fig. 1. The straight blue lines represent the room walls, the small black circles represent the sensors' positions, the red cross represents the original source position, and the black straight line originating from the source represents the source orientation.

We assume that the reflections on the walls of the sound wave emitted from the source are specular, and they will be modeled by using the source image principle: from within the room, the sound wave reflected from a given wall is the same as the one that would be emitted from a source image located outside the room. Fig. 2 represents one such source image, which is obtained by mirroring the original source (and its orientation) w.r.t. the top-right wall, along with its visibility region (green surface), which is included in a cone emerging from the source image position, and whose boundaries match the edges of the reflecting wall. Here we note that the first (i.e. top-left) sensor can see this source image, whereas the second (i.e. bottom-right) one cannot.

This source image is what we call a first generation source image (obtained by mirroring the original source, whose 1 The Matlab code generating all the figures in this report is available at https://perso.telecom-paristech.fr/rbadeau/techreport2019-code.zip. 2 All the experiments in this technical report illustrate the source image principle in 2D-space for convenience, but of course the reverberation model will be defined in the 3D-space. generation index is zero by convention). By iterating this mirroring process on every wall, we can successively obtain infinitely many source image generations. Fig. 3 represents the two constellations of all source images up to the 30th generation, that are visible from the first sensor (red crosses) or from the second one (magenta circles). We observe that some source images are visible from both sensors, whereas others are visible only from one sensor (an accurate count shows that about 32 % of the source images visible from one sensor are also visible from the other one). We also observe that the spatial distribution of the source images looks random, and it is approximately uniform and isotropic. For this reason, this reverberation room is called ergodic, which means that computing averages over time is equivalent to computing averages over space. Indeed, on the one hand, space is sampled uniformly and isotropically by the source images, thus averaging over space is equivalent to averaging over all source images, and on the other hand, the propagation time of sound waves is proportional to the distance of a source image to the sensor, thus every instant corresponds to a circle centered at each sensor, so that averaging over time is also equivalent to averaging over all source images by varying the radius of the circle. of source images as in Fig. 3, but a vertical axis is added, which represents the generation index of every source image. We observe that now both constellations form a cone, that we will refer to as a reflection cone, since the generation index corresponds to the total number of reflections undergone by the sound wave. In the same way, Fig. 6 represents another kind of reflection cone: the vertical axis does no longer represent the generations of the source images, but rather their log-magnitudes, which are computed as follows: the log-magnitude of the original source is 0, and every reflection on every wall decreases the log-magnitude by the attenuation coefficient of the wall. As a numerical example, we used the following values of attenuation coefficients for the four walls: {1, 2, 3, 4} (the walls being sorted counterclockwise, starting from the bottom wall). As previously, the first column corresponds to the first sensor, the second column to the second sensor, the first row represents the cones seen from an oblique direction, the second row represents the cones seen from directly overhead, and the third row represents a slice of the cones along the horizontal axis. We notice that the properties of this reverberation cone are the same as in Fig. 4 and 5: the different attenuation coefficients of the walls do not affect its shape.

Finally, Fig. 7 represents the same kind of reflection cone as in Fig. 6, but this time the attenuation coefficients of the four walls depend on the angle of incidence of the sound wave on the wall (as a numerical example, the coefficients {1, 2, 3, 4} were multiplied by the cosine of the angle of incidence of the sound wave w.r.t. the normal of the wall). We notice that the properties of this reverberation cone are the same as in Fig. 6: the dependence of the attenuation coefficients on the angle of incidence does not affect its shape either.

B. Reverberation model for omnidirectional sources & sensors

We can now introduce a stochastic reverberation model based on these empirical observations. Using the same notation as in [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF], for any sensor i ∈ {1 . . . I} and frequency f ∈ R,

ĥi (f ) = ĝi (f ) ŝ (f ) x∈R 3 y∈R M V i (x, y; y -q x -x i 2 ) e -y α(f )+2ıπf x-x i 2 c x -x i 2 dN (x, y) (2) where:
• For all m = 1 . . . M , the reflections undergone on the m-th wall by the sound wave emitted by the source image at position x are characterized by a unique frequency response e -ym αm (f )

c

, where αm (f ) > 0 is the frequencyvarying attenuation coefficient of the m-th wall (in hertz), and y m ≥ 0 (in meters) is related to the total number of reflections on the m-th wall. Therefore the accumulated frequency response of all reflections on all walls is e -y α(f ) c , where

α(f ) = [α 1 (f ) . . . α M (f )] ∈ R M + and y = [y 1 . . . y M ] ∈ R M + .
• Vectors (x, y) (in meters) are jointly distributed according to a uniform and isotropic Poisson point process dN (x, y) ∼ P(λdxdy) with λ > 0 (in meters -(3+M ) ). • In practice however, not all possible pairs (x, y) are visible from sensor i: in the same way as the source images visible from sensor i form an isotropic reflection cone (as observed in Fig. 4 and5), the number of reflections on every wall is approximately proportional to xx i 2 , hence the condition y m ≈ q m xx i 2 where q m > 0 (dimensionless quantity) is related to the rate of reflection on the m-th wall, which depends on the dimensions of the wall and on the room geometry. • We thus introduce a random Boolean V i (x, y; yq xx i 2 ), where q = [q 1 . . . q M ] ∈ R M + , which determines whether the source image parameterized by (x, y) is visible from sensor i. Its distribution (that we do not need to specify precisely) is parameterized by z i = yq xx i 2 : the closer z i is to zero, the higher the probability of visibility from sensor i (in order to fit the reflection cone). The joint distribution for all sensors i of the random vector [V 1 (x, y; z 1 ) . . . V I (x, y; z I )] i∈{1...I} is assumed i.i.d. w.r.t. (x, y), and it will be denoted p(b 1 . . . b I ; z 1 . . . z I ) where b i ∈ {0, 1} and z i ∈ R M (note that the Boolean variables V i and V j for very close sensors i and j are expected to be mutually dependent). The marginals for one sensor are the same ∀i and will be denoted p(b; z), and the marginals for a pair of sensors (i, j) are the same ∀i, j and will be denoted p(b i , b j ; z i , z j ). We can make several remarks about model (2):

• Contrary to the original model (1), model ( 2) is written in the frequency domain only (not in the time domain), because the dependency of α(f ) on frequency f prevents the inverse Fourier transform of e -y α(f ) c to be written in closed form in the general case. Nevertheless, as we will show in future papers, investigating the statistical properties of model [START_REF] Cook | Measurement of correlation coefficients in reverberant sound fields[END_REF] in the temporal and in the timefrequency domains will be still feasible.

• Considering the above definitions, Proposition 1 in the Appendix proves that the spatial distribution of the source images visible from every sensor i is uniform and isotropic of same rate parameter ∀i, as observed in Fig. 3. • As observed in Fig. 6-(a) to (d), the points of coordinates

x, -y α(f ) c (at fixed frequency f ) form a cone, since the Boolean V i (x, y; yq xx i 2 ) selects values of y close to q xx i 2 . Moreover, since the probability distribution of V i (x, y; yq xx i 2 ) only depends on yq xx i 2 , the thickness of this cone is uniform, as observed in Fig. 6-(e) and (f).

• As observed in Fig. 7, this model is also applicable even when the attenuation coefficient of every wall actually depends on the angle of incidence of the sound wave on the wall. • The attenuation due to sound propagation in the air can be modeled in exactly the same way as the attenuation due to the successive reflections on any wall; so it is readily taken into account by the model, if we just assume that one index (e.g. m = 1) corresponds to the propagation in the air rather than to the reflections on a wall of the room.

• Even though the source image principle only takes specular reflections into account in its geometric formulation, the stochastic model ( 2) allows for source images randomly distributed in space, with random visibility regions; therefore it readily takes diffuse reflections into account. In brief, compared to the original model (1), model (2) brings the following novelties:

• the attenuation coefficients are now frequency-varying;

• the reflection cone is no longer parameterized as a perfect isotropic cone (cf. the term e -α x-x i 2 c in (1)), but more realistically as a noisy isotropic cone (cf. the term e -y α(f ) c in [START_REF] Cook | Measurement of correlation coefficients in reverberant sound fields[END_REF], where y is random);

• the visibility regions of the source images from different parts of the room are taken into account.

III. UNIFORM AND ISOTROPIC MODEL FOR DIRECTIVE SOURCES AND SENSORS

A. Simulation results: orientations of the source images

We now intend to take the source directivity into account. Fig. 2 shows that the orientations of the source images are obtained by iteratively mirroring that of the original source. Fig. 8 represents the resulting orientations of the source images in a vicinity of the original room. It can be noticed that, similarly to the positions of the source images, their orientations also look random. Suppose now that the original source has a directivity, so that its response can be written ŝ (u, f ) instead of ŝ(f ), where u is a unit direction vector. Then at any sensor i, the response of a source image of parameters (x, y) can be written ŝ Θ(x, y) xi-x xi-x 2 , f , where Θ(x, y) is a random rotation matrix, i.i.d. w.r.t. (x, y) and uniformly distributed on the rotation group SO(3).

B. Reverberation model for directive sources and sensors

In order to take the source and sensor directivities into account, model ( 2) is modified as follows: for any sensor i ∈ {1 . . . I} and frequency f ∈ R,

ĥi (f ) = x∈R 3 y∈R M V i (x, y; y -q x -x i 2 ) ĝi Θ i x-xi x-xi 2 , f ŝ Θ(x, y) xi-x xi-x 2 , f e -y α(f )+2ıπf x-x i 2 c x-xi 2 dN (x, y) (3)
where

• ĝi (u, f ) ∈ C is the response of sensor i at direction u and frequency f (taking into account its directivity); • ŝ (u, f ) ∈ C is the response of the source at direction u and frequency f (taking into account its directivity); 3) is a deterministic rotation matrix that represents the orientation of sensor i; • Θ(x, y) ∈ SO(3) is a random rotation matrix that represents the orientation of the source image of parameters (x, y). Its distribution is i.i.d. w.r.t. (x, y) and uniform on the rotation group SO(3).

• Θ i ∈ SO(

IV. UNIFORM AND ANISOTROPIC MODEL A. Simulation results: rectangular room

We now consider the rectangular room depicted in Fig. 10 (using the same symbols as in Fig. 1).

Fig. 11 represents the constellation of all source images up to the 30th generation. Note that in this very particular geometric setup, all source images are visible from the whole room 3 . Compared with Fig. 3, the spatial distribution of the source images no longer looks random, but it is still uniform. 11, but a vertical axis was added, which represents the generation index of every source image. We observe that the constellation still forms a reflection cone, but this cone is no longer isotropic. Fig. 12-(c) represents a slice of the cone cut along the horizontal axis, which shows that the thickness of the cone is still approximately uniform. Fig. 12-(b) represents the same reflection cone as in Fig. 12-(a), as a surface plot obtained by means of Delaunay triangulation. Fig. 12-(d), represents the same reflection cone as in Fig. 12-(b), but seen from directly overhead, which confirms that the shape of the cone is not isotropic. The left column in Fig. 13 repeats the experiment in Fig. 6: the vertical axis does no longer represent the generations of the source images, but rather their log-magnitudes, which are computed as follows: the log-magnitude of the original source is 0, and every reflection on a wall decreases the logmagnitude by the attenuation coefficient of this wall. As a numerical example, we used the same values of attenuation coefficients for the four walls: {1, 2, 3, 4} (the walls being sorted counter-clockwise, starting from the bottom wall). As previously, the first row represents the cone seen from an oblique direction, the second row represents the cone seen from directly overhead, and the third row represents a slice of the cone along the horizontal axis. We notice that the properties of this reverberation cone are the same as in Fig. 12: the different attenuation coefficients of the walls do not affect its shape.

The right column in Fig. 13 repeats the experiment in Fig. 7: now the attenuation coefficients of the four walls depend on the angle of incidence of the sound wave on the wall (as a numerical example, the coefficients {1, 2, 3, 4} were multiplied by the cosine of the angle of incidence w.r.t. the normal of the wall). Contrary to what was observed in Fig. 7, we notice that the shape of the reverberation cone is modified by the dependence of the attenuation coefficients on the angle of incidence: the cone is more isotropic. But more importantly, it is still a cone.

Finally, Fig. 14 represents the orientations of the source images in a vicinity of the original room. It can be noticed that, similarly to the positions of the source images, their orientations no longer look random: their angles take 4 different values only, that are equidistributed.

B. Reverberation model with anisotropic reflection cone

In order to take the anisotropy observed in Section IV-A into account, model ( 3) is modified by simply replacing the isotropic reflection cone, which was parameterized as y ≈ q x-x i 2 , by a more general, possibly anisotropic, reflection cone, parameterized as y ≈ q(xx i ): for any sensor i ∈ {1 . . . I} and frequency f ∈ R,

ĥi (f ) = x∈R 3 y∈R M V i (x, y; y -q(x -x i )) ĝi Θ i x-xi x-xi 2 , f ŝ Θ(x, y) xi-x xi-x 2 , f e -y α(f )+2ıπf x-x i 2 c x-xi 2 dN (x, y) (4)
where

• q(x) ∈ R M + is a 1-homogeneous function;

• Θ(x, y) is a random rotation matrix, i.i.d. w.r.t. (x, y), whose probability distribution is not necessarily uniform on the rotation group SO(3). Compared to the original model, this new model permits to represent a uniform and anisotropic acoustic field.

V. NON-UNIFORM AND ANISOTROPIC MODEL

A. Simulation results: church geometry

We consider a bi-dimensional room whose geometry is shaped as that of church, and depicted in Fig. 15. Fig. 16 represents the two constellations of all source images up to the 30th generation, that are visible from the first (i.e. top-left) sensor (red crosses) or from the second (i.e. bottom-right) one (magenta circles). As in Section II-A, we observe that some source images are visible from both sensors, whereas others are visible only from one sensor (an accurate count shows that about 7 % of the source images visible from one sensor are also visible from the other one). 16, but a vertical axis was added, which represents the generation index of every source image. Again we observe that both constellations form a reflection cone. In the same way, Fig. 17-(c) (resp. Fig. 17-(d)) represents a slice of the cone in Fig. 17-(a) (resp. Fig. 17-(b)) cut along the horizontal axis, which shows that the thickness of the cone is approximately uniform. Fig. 19 represents the second kind of reflection cones: the vertical axis represents the log-magnitudes of the source images, which are computed as follows: the log-magnitude of the original source is 0, and every reflection on a wall de-creases the log-magnitude by the attenuation coefficient of this wall. As a numerical example, we used the following values of attenuation coefficients for the 14 walls: {1, . . . , 14} (the walls being sorted counter-clockwise, starting from the bottom wall). As previously, the first column corresponds to the first sensor, the second column to the second sensor, the first row represents the cones seen from an oblique direction, the second row represents the cones seen from directly overhead, and the third row represents a slice of the cones along the horizontal axis. We notice that the properties of these reverberation cones are the same as in Fig. 17 Fig. 20 represents the same kind of reflection cones as in Fig. 19, but now the attenuation coefficients of the 14 walls depend on the angle of incidence of the sound wave on the wall (as a numerical example, the coefficients {1, . . . , 14} were multiplied by the cosine of the angle of incidence w.r.t. the normal of the wall). As in Section IV-A, we notice that the shapes of the reverberation cones are modified by the dependence of the attenuation coefficients on the angle of incidence: they are more isotropic. But more importantly, they are still cones.

Finally, Fig. 21-(a) (resp. 21-(b)) displays the histogram of the orientations of all source images up to the 30th generation, that are visible from sensor 1 (resp. sensor 2). These histograms look a bit noisy, but clearly show that the distribution of orientations of the source images is the same for both sensors, and that it is not uniform (the angles take 8 different values only).

B. Reverberation model with sensor-dependent reflection cone

In order to take this dependence of the reflection cone on the sensor location into account, model ( 4) is modified by simply making the 1-homogeneous function q depend on x i : for any

ĥi (f ) = x∈R 3 y∈R M V i (x, y; y -q(x -x i , x i )) ĝi Θ i x-xi x-xi 2 , f ŝ Θ(x, y) xi-x xi-x 2 , f e -y α(f )+2ıπf x-x i 2 c x-xi 2 dN (x, y)
where function q(x, x i ) is 1-homogeneous w.r.t. its first parameter. Compared with the original model, this new stochastic model permits to represent non-uniform acoustic fields: different constellations of source images may form different reflection cones in different parts of the room (however, note that the random variables V i (x, y; z), Θ(x, y) and dN (x, y) are still i.i.d. for all (x, y); they are thus uniformly distributed over space).

VI. GENERAL STOCHASTIC REVERBERATION MODEL

To sum up, we can now formally define the most general stochastic reverberation model, that was introduced in Sec-tion V-B: for any sensor i ∈ {1 . . . I} and frequency f ∈ R,

ĥi (f ) = x∈R 3 y∈R M V i (x, y; y -q i (x -x i )) ĝi Θ i x-xi x-xi 2 , f ŝ Θ(x, y) xi-x xi-x 2 , f e -y α(f )+2ıπf x-x i 2 c x-xi 2 dN (x, y) (5)
where

• ĥi (f ) ∈ C is the room frequency response at sensor i as a function of frequency f . • ĝi (u, f ) ∈ C is the response of sensor i at direction u and frequency f (taking into account its directivity). • ŝ (u, f ) ∈ C is the response of the source at direction u and frequency f (taking into account its directivity). • vector x i ∈ R 3 (in meters) is the position of sensor i.

• vector x ∈ R 3 (in meters) represents the positions of the source images.

α(f ) ∈ R M
+ is a vector of frequency-varying attenuation coefficients (in hertz). . The marginals for one sensor are the same ∀i and will be denoted p(b; z) (they are such that the closer z is to zero, the higher p(1; z)), and the marginals for a pair of sensors (i, j) are the same ∀i, j and will be denoted p(b i , b j ; z i , z j ).

• q i (x) = q(x, x i ) ∈ R M + is a 1-homogeneous function. • vector (x, y) ∈ R 3 × R M (in
This general model encompasses all the particular cases of interest presented in the previous sections:

• if the acoustic field is uniform, q i (x) = q(x) ∀i;

• if the acoustic field is isotropic, q i (x) = q i x 2 and rotation matrices Θ(x, y) are uniformly distributed; • therefore if the acoustic field is diffuse, q i (x) = q x 2 ;

• if sensor i is omnidirectional, ĝi (u, f ) does no longer depend on u; • if the source is omnidirectional, ŝ (u, f ) does no longer depend on u and rotation Θ(x, y) disappears.

Note that in practice, functions ĝi (u, f ), ŝ (u, f ), α(f ), and the parametric probability distribution p, have to satisfy some technical conditions for the model ( 5) to be mathematically well-defined and physically meaningful. These conditions depend on the physical assumptions about the acoustic field (either isotropic or anisotropic, either uniform or non-uniform), the source and the sensors (either directive or omnidirectional). They will be mathematically formulated in future work.

VII. CONCLUSION AND PERSPECTIVES

In this technical report, we proposed several extensions of the stochastic reverberation model introduced in [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF], that aim to model reverberation more realistically, by considering anisotropic and non-uniform acoustic fields, directive sources and microphones, and frequency-varying attenuation coefficients. In future work, we will show that this model is able to predict various statistical properties of reverberation (regarding the first and second order moments in the spatial, spectral, temporal and time-frequency domains, and the asymptotic Gaussianity), that we experimentally observed in a wide variety of both synthetic and measured room impulse responses. APPENDIX Proposition 1. Let dN (x, y) denote a uniform Poisson random measure on R 3 × R M with independent increments dN (x, y) ∼ P(λdxdy) where λ > 0, and ∀(x, y) ∈ R 3 × R M , let V (x, y; z) be independent Boolean random variables, whose probability distribution p(b; z) with b ∈ {0, 1} is parameterized by z ∈ R M , and such that z → p(1; z) ∈ L 1 (R M ). Finally, let q : R 3 → R M + be any function. Then

dN (x) = y∈R M
V (x, y; yq(x))dN (x, y)

defines a uniform Poisson random measure on R 3 , of rate parameter λ = λ z∈R M p(1; z)dz.

Proof of Proposition 1. In Proposition 4 in [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF], we showed that a random measure is distributed as a uniform Poisson point process with independent increments dN (x) ∼ P(λ dx) if and only if any stochastic integral of the form

I = x∈R 3 ψ(x) dN (x), (7)
where ψ : R 3 → R is any essentially bounded function with compact support, has a characteristic function of the form φ I (θ) = e λ x∈R 3 (e ıθψ(x) -1)dx [START_REF] Polack | Modifying chambers to play billiards: the foundations of reverberation theory[END_REF] with λ > 0. Thus, in order to prove that (6) defines a uniform Poisson random measure on R 3 , we can sketch the following mathematical derivations:

ln (φ I (θ)) = ln E[e ıθI ] (9)
=ln E e ıθ x∈R 3 ψ(x) y∈R M V (x,y;y-q(x))dN (x,y) (10)

= x∈R 3 y∈R M
ln E e ıθψ(x)V (x,y;y-q(x))dN (x,y) (11)

=

x y ln 1+p(1; y-q(x)) E e ıθψ(x)dN (x,y) -1 (12) = x∈R 3 y∈R M p(1; y-q(x)) E e ıθψ(x)dN (x,y) -1 (13) = x∈R 3 y∈R M p(1; y-q(x)) e λdxdy(e ıθψ(x) -1) -1 (14) =λ x∈R 3 y∈R M p(1; y-q(x)) e ıθψ(x) -1 dxdy (15)

=λ x∈R 3
e ıθψ(x) -1 dx (16) where

• [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF] is obtained from (9) by substituting ( 7) and ( 6);

• [START_REF] Kuttruff | Room Acoustics, Fifth Edition[END_REF] is obtained from [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF] by considering the independence for all (x, y) of the random variables V (x, y; yq(x)) and dN (x, y); • (12) is obtained from [START_REF] Kuttruff | Room Acoustics, Fifth Edition[END_REF] by considering the conditional expectation of e ıθψ(x)V (x,y;y-q(x))dN (x,y) given the Boolean random variable V (x, y; y-q(x)); • (13) is obtained from (12) by using the first-order expansion of the logarithm; • (14) is obtained from (13) by using the closed-form expression of the characteristic function of the Poisson distribution; • (15) is obtained from (14) by using the first-order expansion of the exponential; • (16) is obtained from (15) by substituting z = yq(x) and λ = λ z∈R M p(1; z)dz. These derivations do not really form a rigorous mathematical proof of ( 8), but they provide an intuitive understanding of how it is obtained. Actually, (8) can be formally proved by following the construction of the Lebesgue integral, as we did in [START_REF] Badeau | Common mathematical framework for stochastic reverberation models[END_REF]Appendix B]. We conclude that dN (x) defines a uniform Poisson random measure on R 3 with independent increments dN (x) ∼ P(λ dx).

Fig. 1 .Fig. 2 .

 12 Fig. 1. Geometry of the ergodic room

Fig. 4 -

 4 (a) and (b) represent the exact same constellations

Fig. 3 .

 3 Fig. 3. Two constellations of source images up to 30th generation

Fig. 4 -

 4 (c) (resp. Fig. 4-(d)) represents a slice of the cone in Fig. 4-(a) (resp. Fig. 4-(b)) cut along the horizontal axis, which shows that the thickness of the cone is approximately uniform.

Fig. 4 .

 4 Fig. 4. Constellations of source images with their generations

Fig. 5 .

 5 Fig. 5. Reflection cones of the generations of the source images

Fig. 6 .Fig. 7 .

 67 Fig.[START_REF] Schroeder | Statistical parameters of the frequency response curves of large rooms[END_REF]. Reflection cones of the log-magnitudes of the source images (attenuations do not depend on the angle of incidence on the walls)

Fig. 8 .

 8 Fig. 8. Orientations of the source images

Fig. 9 .

 9 Fig. 9. Histograms of the orientations of the source images

Fig. 10 .

 10 Fig. 10. Geometry of the rectangular room

Fig. 11 .

 11 Fig. 11. Constellation of source images up to 30th generation

Fig. 12 .

 12 Fig. 12. Reflection cone of the generations of the source images

Fig. 13 .

 13 Fig. 13. Reflection cone of the log-magnitudes of the source images

Fig. 15 .

 15 Fig. 15. Geometry of the church

Fig. 16 .

 16 Fig. 16. Two constellations of source images up to 30th generation

Fig. 17 .

 17 Fig. 17. Constellations of source images with their generations

Fig. 18 -

 18 Fig. 18-(a) (resp. Fig. 18-(b)), represents the same reflection cone as in Fig. 17-(a) (resp. Fig. 17-(b)), as a surface plot obtained by means of Delaunay triangulation. Fig. 18-(c) (resp. Fig. 18-(d)) represents the same reflection cone as Fig. 18-(a) (resp. Fig. 18-(b)), but seen from directly overhead. Contrary to what we observed in Section II-A, the two reflection cones are neither isotropic nor equal: the reflection cone now depends on the sensor location.

Fig. 18 .

 18 Fig. 18-(a) (resp. Fig. 18-(b)), represents the same reflection cone as in Fig. 17-(a) (resp. Fig. 17-(b)), as a surface plot obtained by means of Delaunay triangulation. Fig. 18-(c) (resp. Fig. 18-(d)) represents the same reflection cone as Fig. 18-(a) (resp. Fig. 18-(b)), but seen from directly overhead. Contrary to what we observed in Section II-A, the two reflection cones are neither isotropic nor equal: the reflection cone now depends on the sensor location.

  and 18: the different attenuation coefficients of the walls do not clearly affect their shapes.

Fig. 19 .

 19 Fig. 19. Reflection cones of the log-magnitudes of the source images (attenuations do not depend on the angle of incidence on the walls)

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Reflection cones of the log-magnitudes of the source images (attenuations depend on the angle of incidence on the walls)

  • Θ(x, y) ∈ SO(3) is a random rotation matrix that represents the orientation of the source image of parameters (x, y). Its distribution is i.i.d. w.r.t. (x, y) and not necessarily uniform on the rotation group SO(3).• c > 0 is the speed of sound (in meters.hertz).

•

In practice, the source images still have visibility regions, but the particular geometry of the room is such that every source image in Fig.11actually results from the superposition of a finite number of source images of same generation and same orientation, located at the same position, whose visibility regions form a partition of the rectangular room.