Lower bound on frequency validity of energy-stress tensor based diffuse
sound field model
Aidan MEACHAM(1) , Roland BADEAU(2) , Jean-Dominique POLACK(1)(1) Sorbonne Université, UMR CNRS 7190, Institut Jean le Rond ∂ ’Alembert, France(2) LTCI, Télécom ParisTech, Université Paris-Saclay, France(1) aidan.meacham@sorbonne-universite.fr

Abstract
A lower bound on the frequency validity limit is established for an energetic wave equation derived from theenergy-stress tensor, examined in the one-dimensional case [Dujourdy et al, Acta Acustica united with Acustica
103:480-491, 2017]. The method efficiently models diffuse sound fields that dominate reverberation at higher
frequencies and larger distances. Initially noted in the course of an exhaustive search of the solution space of all
valid model parameters, the low-frequency cutoff has implications for the utility of the method in a hybridization
context. In practice, the bound is encountered when determining the absorption and diffusion coefficients by
iteratively approaching the temporal and spatial decay of measured data. As the test frequency decreases, theranges of coefficient combinations that result in less than 10% variation from each decay measure can diverge
until the region where both measures are satisfactory (the intersection of the two domains) disappears. Further
evidence for the bound is provided through comparison with measurements of a long hallway, and stability
concerns in the cases where both coefficients are very small are addressed.
Keywords: Room acoustics, finite difference methods, diffuse field

1

INTRODUCTION

While geometric and wave equation based approaches to room acoustics have dominated theoretical research,in the realm of practical implementations, hybrid methods combining both strategies are increasingly popular.
Since each approach has strengths and weaknesses depending on the frequency range that is being simulated,
it has long been the goal of practicioners to leverage different types of simulations in a cohesive manner tocreate full-band predictions with efficiency and accuracy that would not be possible with a single type alone.
For example, geometric acoustics simulations such as ray tracing (1) or the image source method (2) efficiently
model the early part of an impulse response, even at high frequencies, but cannot account for diffraction ina physically meaningful way, an important part of the low frequency response, and become inefficient when
considering complex enclosure geometry or when directly modeling the late reverberation of a space. On theother hand, wave equation based methods (such as finite element or volume methods in both the time andfrequency domains) excel in the low frequency regime, since they simulate the sound field everywhere in theenclosure, including diffraction and directional effects, with relatively direct implementations of sources andreceivers (3, 4). Extending wave modeling methods to higher frequencies, however, requires increasingly finemeshes that become computationally untenable for real-time applications, even with the advent of modern GPU
implementations (5).
From this brief overview, it is not difficult to predict how hybrid approaches could combine these methods,
utilizing each method in its most effective domain in order to create a full-band model. In practice, however,there are a number of difficulties that arise from combining separate methods in this fashion, not the least ofwhich is the maintenance of machinery (including but not limited to numerical code, meshing tools, automatic
matching of boundary conditions across simulation types, room geometry simplifications, and post-processing
DSP) for each separate part of the simulation. In this paper, we consider how an energetic wave equation based
method mitigates some of these concerns while considering its limitations for future hybridization.

6012

This method, developed by Dujourdy et al. (6), has been found to accurately model the statistical late reverberation, or diffuse field, resulting from an initial perturbation, in both 1- and 2-dimensional contexts (7).
Originating from Ollendorff and Picaut’s work (8, 9), this method aims to take advantage of an additional
degree of freedom given by the inclusion of a diffusion coefficient in addition to the traditional absorption coefficient. Previous work on this topic considered only the 1000 Hz frequency band, but in order to consider
its future fitness for hybrid modeling, we will consider in greater detail the frequency ranges it is capable ofmodeling. Such information provides bounds on the crossover frequencies that would be required in a hybrid
context, implying the necessary ranges that would need to be simulated with other methods in order to createa full-band model, and thus whether or not the diffuse field strategy can be useful in a practical architectural
acoustics context.

2

ENERGETIC WAVE EQUATION METHOD

For this study, we directly used the theory and implementation of the numerical schemes in (6). For the sakeof brevity, we will not reproduce those findings here, but encourage those desiring a full treatment to refer tothe original work.
In brief, the main idea of the energetic wave equation is the development of two conservation laws, one eachfor energy density E and sound intensity J, as defined by Morse and Feshback in (10), rather than the typicalenergy balance based on Gauss’s theorem found in typical finite difference or finite volume formulations. Using
these two laws, it is possible to derive a system of coupled equations relating the energy density, the soundintensity, and the wave-stress symmetric tensor E (11). In the case where one length dominates a space andthe cross-sectional area is relatively constant (as is the case in a long hallway), this system of equations can be
reduced to one dimension by introducing modified absorption and scattering coefficients and integrating (withfurther energy and momentum balance hypotheses) over the minor axis walls.
After dimensional reduction, the system resembles the telegrapher’s equations, and by inspection can be transformed into a linear second-order hyperbolic equation with a single dependent variable (in this case, the energy
density E). Finally, boundary conditions can be derived using the aforementioned hypotheses, and the entire
wave propagation system can be discretized with common finite difference time domain strategies.
The main advantage of this approach is that the propagation of diffuse energy allows for very large spatial discretization. The reason for this is that the modulation frequency of the late energy decay that we are concerned
about is very low. Typical acoustic wave simulations that are concerned with high frequency content have a very
high modulation frequency (since they specifically want to resolve the individual pressure waves), but since weare primarily interested in the stochastic decay of the diffuse high frequency energy, we can assume that the
decay itself is not changing very rapidly, and thus can accommodate relatively large spatial sampling. As an
example, if we are only concerned with the diffuse energy level every tenth of a second, then the spatial sampling rate can be as coarse as 3 meters, with the temporal sampling rate chosen to satisfy the scheme’s stability
conditions.

3

FREQUENCY VALIDITY ANALYSIS

3.1 Room geometry
A similar corridor to the original was used for this study. The hallway had an overall length of 45 meters, witha width of 159 cm and a height of 237.5 cm. In the main narrow portion of the hallway, the ceiling height wasmeasured to a fine metal grating suspended below the actual height of the corridor, which carried lighting andconduit. The hard ceiling was 326 cm, with a decrease every 1.5 meters for metal support beams to 280 cm.
In the recesses, all of which were of uniform length and depth (except for the doorway furthest on the right ofthe plan), the width increased to 239 cm and the height decreased to 220 cm, and ceiling was masonry rather
than the grating mentioned above. In some of the recesses, there were glass display cases or small pieces of6013

Figure 1. Floorplan for the corridor under consideration.
furniture that contributed to the diffusive effects of the recesses themselves. The entire floor was linoleum, andfor the most part, the walls were wooden panel and masonry. There were occasional metal gratings on the flatwall for HVAC, which was audible but not distracting. All doors entering the hallway were closed, and thedoors in the center of the hallway were fully open.
3.2 Measurements
Table 1. T60 s and Spatial Decays
Frequency [Hz]T60 [S]
Spatial Decay [dB/m]

62.5
2.29
-0.63

125
0.35
-0.49

250
0.38
-0.88

500
0.39
-0.94

1000
0.43
-0.66

2000
0.35
-0.54

4000
0.33
-0.66

8000
0.30
-0.84

Impulse responses were collected using a SoundField ST250 microphone and an Outline GRS omnidirectional
speaker, with a MOTU Traveler sound card. As before, the source was positioned 1 meter away from the endof the hall, 1.5 meters above the ground, and centered between the two walls. Beginning 1 meter from thesource, measurements were collected with the microphone’s X-axis aligned along the length of the hallway. A
spacing of 1 meter was used out to 10 meters, which corresponds exactly to the discretization distance in the
numerical simulations. Then, recordings were made every 2 meters until the end of the hallway for a total of
26 sampling locations. Recordings were made using the the swept sine method (12, 13) as implemented in the
Adobe Audition plugin Aurora. The sweep length was 20 seconds, and the source level was adjusted digitally
to maximize the signal-to-noise ratio without clipping as the microphone was moved further and further fromthe source. These gains were recorded in order to recover the true measured energy level for each measurement
location. The sweep responses were then post-processed by convolution with the inverse sweep to recover
impulse responses for each location.
The relevant frequency bands were selected by filtering each impulse response with a standard octave-band
filterbank, and then calculating the desired metrics for each resulting bandlimited response. The T60 (a measureof temporal decay) was calculated by a linear fit to the Energy Decay Curve in dB for each frequency band andat each receiver location. This is sometimes known as Schroeder’s reverse integration (14). Next, each receiverlocation was averaged to arrive at a single T60 for the entire hallway in each band. Finally, spatial decays(similar to strength of sound, or G) were calculated by a linear fit to the sum of energy for each band-limited
response across each receiver location. Both sets of calculated values are presented in Table 1.
3.3 Simulations
In order to determine the frequency bands the model was able to represent, we ran simulations for a sampledsubset of all combinations of absorption and scattering coefficients α and β . In a typical scenario, a scheme toconverge on the desired spatial and temporal characteristics would be used to match a single frequency band,
but here, we want to visualize the entire problem space.

6014

62.5 Hz125 Hz250 Hz500 Hz2000 Hz

4000 Hz

8000 Hz

0.5

Diffraction Coefficient

Regions of Validity
T60

0.4

Spatial Decay
Both

0.3

0.2

0.1

0

1000 Hz

Diffraction Coefficient

0.5

0.4

0.3

0.2

0.1

0

0

0.1

0.2

0.3

Absorption Coefficient

0.4

0

0.1

0.2

0.3

Absorption Coefficient

0.4

0

0.1

0.2

0.3

Absorption Coefficient

0.4

0

0.1

0.2

0.3

0.4

Absorption Coefficient

Figure 2. Agreement between simulated and measured data for the hallway under consideration.
For this study, as was the case before, the spatial sampling step was chosen to be ∆x = 1 meter.
The initial conditions were chosen to be a temporal Gaussian, centered at the same location as the source, inorder to minimize spurious numerical oscillations. Afterward, the same acoustical indices as with the in situ
measurements were extracted from each result.
3.4 Valid domains
Finally, we compared the measured and simulated results for each index. For each frequency band (presented
here in ascending order), we classified which simulations were within 10% of each desired metric, the spatial
and temporal decay rates. To preserve legibility, we have represented in Figure 2 only the region with the most
relevant combinations of coefficients for this hallway. The absorption coefficient α ranges from 0.01 to 0.45,
and the diffusion coefficient β ranges from 0.01 to 0.5. This region corresponds well to the typical ranges
where both T60 and spatial energy decays are valid.
As expected, the two coefficients have different effects on the simulations, given that the validity patterns are
not symmetric, confirming the observation in the original paper. Furthermore, the regions of validity for eachindividual measure appear to be smoothly varying. Perhaps most importantly, for some low frequencies, there
is no region where both indices are valid, implying a lower bound on the frequency that can be represented by
the model for this hallway.
At the lowest frequency band, 62.5 Hz, there are regions of validity outside of the chosen region, but we havechosen to exclude these as either their absorption or diffusion coefficients are unrealistically small or large.
Since the regions of validity appear to smoothly vary, it may be the case that with a particular frequency inbetween this band and 125 Hz, there is a point where the spatial decay and T60 curves “switch” places andbecome valid. We also note the similarity of the 125 Hz band and the 2000 Hz band, since the region ofvalidity is determined by the 10% threshold. In other bands, the regions cross, and therefore can be uniquely

6015

determined to match both measures to arbitrary precision without the need of a threshold. Both of these topics
related to a more continuous perspective on the methodology for determining validity deserve closer examinationin future work.
The case where the individual metric regions cross resulting in two disconnected regions of validity, up to
the lack of symmetry discussed above, can be explained by the observation by Dujourdy et al. (6) that thetwo coefficients can typically be exchanged, where the smaller of the two acts as the absorption coefficient.
While in general this situation is preferable, it implies that the optimization problem, at least in its current
formulation, is not convex, and furthermore, more information about the problem will be required in order touniquely determine a preferred combination of coefficients.
One possible explanation for the lack of a valid region at the lowest frequencies is the mismatch between
long T60 and relatively low spatial decays. At low frequencies, it is common to have more reverberation,and similarly, the energy may be more evenly distributed as a result of diffusion. At higher frequencies, theopposite is the case, where increased absorption of high-frequency energy results in shorter reverberation times
and steeper spatial decay, both of which appear to be more readily represented by this model.
In the case where either coefficient is very small, generally speaking, the simulations become closer and closer
to the lossless case, extending the response times and magnifying small perturbations. As a result, simulations
with some combinations of α and β can veer into instability. This would further constrain the region ofvalidity, but in practice, instability only occurs when approaching the most extreme values of T60 or spatial
decay. Thus, such instability can simply serve as a indicator that careful attention should be paid to the physicalmeasurements that are being matched. It appears that since the region of validity at the lowest frequencies isreasonably well balanced between the two coefficients, for this type of hallway, instability should not present a
major concern in the context of hybridization since the most dangerous simulations can simply be discarded.
3.5 Implications
Since the lowest frequency that contains a convincing overlap is 250 Hz, for this corridor, that is the lowest
bound on the validity of the model. This makes sense, seeing as the late reverberation is more statistical athigher frequencies. For small concert halls, 250 - 500 Hz is a very convenient cutoff frequency for a potential
hybridization since a complementary pressure wave time domain approach would be capable of modeling theregion in real time (5). Thus, at least for this case, it is worthwhile to consider the model as a possible
candidate for accelerating the synthesis of the late tail in a full-band audio simulation.
3.6 Other hallways
In order to confirm these results, the procedure was followed using the recordings from the previous hallway.
In the interest of space, we will not fully reproduce the measurements as we have with the current study, but
the findings were consistent with the conclusions we have drawn. The corridor from the initial work exhibitedregions of validity down to 125 Hz. Given that the two were relatively similar in terms of shape and length,
this is perhaps unsurprising, but a useful confirmation nonetheless.
A third unrelated hallway of a similar length but with fewer diffusing surfaces was also measured and was
only representable (with “crossed” validity regions) by the model in the 125 Hz and 1000 Hz bands. Upon
initial review of these results, it appears that the reason for the lack of validity is primarily due to the reduced
spatial decays across all bands. Because the hallway did not have as many absorptive surfaces, the sound energy
was much more evenly distributed throughout or even modal, and therefore presented difficulty for the model,
which is predicated on the assumption of spatial decay. One possible interpretation is that this incompatibility
between temporal and spatial decays is a manifestation of a truly uniform reverberant field in regions far fromthe source, where the so-called diffuse field theory is valid. Thus, the case where the energy-based model fails
may actually be a criterion for spaces where the Sabine equation holds.
Figures displaying the regions of validity for both of these corridors can be found in the appendix.

6016

4

CONCLUSIONS AND FUTURE WORKIn this paper, we have discussed our findings regarding the lower frequency bound for which an energy-basedlate-reverberation model can accurately represent the stochastic soundfield. The current study focused on a onedimensional model of a long hallway for which the lower bound was 250 Hz. This means that above thiscutoff frequency, the model was able to recreate a diffuse soundfield with acoustical indices consistent withmeasurements made in a real hallway within 10% error.
The investigation of this lower frequency bound brings to mind a number of questions for future research. In
this case, we have examined a single hallway with a given length and certain acoustical indices, but it should bepossible to extend the analysis to one dimensional problems of varying lengths, T60 , and spatial decays, such thatany problem dominated by its length with a sufficient amount of diffusion could be verifiably modeled above
a certain frequency. Furthermore, given the discovery that in some frequency bands there are multiple regionsof validity, it should be possible to examine more optimal methods for converging on desired coefficients for agiven physical measurement in an interative fashion. In such a fashion, it may be possible to better characterize
the contours of the underlying search space in order to accelerate convergence.
In terms of bringing the theoretical model into practical application, it is imperative to begin experimenting witha sonification procedure. While the theoretical aspects of the method are interesting in and of themselves, thereal value of such an approach is the additional accuracy and efficiency that would be enjoyed upon integration
with a full-band hybrid auralization system. To that end, the next step in developing the model is extending thecurrent 1- and 2-dimensional models to 3 dimensions.
Finally, given the large spatial discretization step in this method, it may be the case that the energy modeland a pressure wave approach could use the same grid, reducing the meshing that must be completed before
simulations can be performed. Verifying that meshes (or subsets thereof) could be trivially reused for bothmethods in a hybrid context would be especially relevant in the cases where room geometry is time varying, apossibility that deserves closer inspection.

REFERENCES
1.

Krokstad A, Strom S, and Sørsdal S. Calculating the Acoustical Room Response by the Use of a Ray
Tracing Technique. Journal of Sound and Vibration 1968;8:118–25.

2.

Allen JB and Berkley DA. Image Method for Efficiently Simulating Small-room Acoustics. The Journal of
the Acoustical Society of America 1979;65:943–50.

3.

Savioja L, Rinne TJ, and Takala T. Simulation of Room Acoustics with a 3-D Finite Difference Mesh. In:
Proc. ICMC. 1994.

4.

Botteldooren D. Acoustical Finite-difference Time-domain Simulation in a Quasi-Cartesian Grid. The Journal of the Acoustical Society of America 1994;95:2313–9.

5.

Savioja L. Real-Time 3D Finite-Difference Time-Domain Simulation of Low- and Mid-Frequency Room
Acoustics. In: DAFx. 2010:8.

6.

Dujourdy H, Pialot B, Toulemonde T, and Polack JD. An Energetic Wave Equation for Modelling Diffuse
Sound Fields – Application to Corridors. Acta Acustica united with Acustica 2017;103:480–91.

7.

Dujourdy H, Pialot B, Toulemonde T, and Polack JD. An Energetic Wave Equation for Modelling Diffuse
Sound Fields—Application to Open Offices. Wave Motion. Innovations in Wave Modelling II 2019;87:193–
212.

8.

Ollendorff F. Statistical Room-Acoustics as a Problem of Diffusion (A Proposal). Acta Acustica united
with Acustica 1969;21:236–45.

9.

Picaut J, Simon L, and Polack JD. A Mathematical Model of Diffuse Sound Field Based on a Diffusion
Equation. Acta Acustica united with Acustica 1997;83:614–21.

6017

6018

10.

Morse PM and Feshbach H. Methods of Theoretical Physics. Mc Graw-Hill Book Company, 1953.

11.

Morse PM and Ingard KU. Theoretical Acoustics. Princeton University Press, 1968. 949 pp.

12.

Farina A. Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique.
In: Audio Engineering Society Convention 108. Audio Engineering Society, 2000.

13.

Farina A. Advancements in Impulse Response Measurements by Sine Sweeps. In: Audio Engineering Society Convention 122. Audio Engineering Society, 2007:21.

14.

Schroeder MR. New Method of Measuring Reverberation Time. The Journal of the Acoustical Society ofAmerica 1965;37:409–12.

APPENDIX
62.5 Hz125 Hz250 Hz500 Hz2000 Hz

4000 Hz

8000 HzRegions of Validity

0.6

Diffraction Coefficient

T60
Spatial Decay
Both

0.5
0.4
0.3
0.2
0.1
0

1000 Hz

Diffraction Coefficient

0.6
0.5
0.4
0.3
0.2
0.1
0

0

0.2

0.4

0.6

0

Absorption Coefficient

0.2

0.4

0.6

0

Absorption Coefficient

0.2

0.4

0.6

0

Absorption Coefficient

0.2

0.4

0.6

Absorption Coefficient

Figure 3. Agreement between simulated and measured data for the hallway studied in Dujourdy et al.
62.5 Hz0.25

125 Hz250 Hz500 Hz2000 Hz

4000 Hz

8000 HzRegions of Validity

Diffraction Coefficient

T60

0.2

Spatial Decay
Both

0.15

0.1

0.05

0

1000 Hz

Diffraction Coefficient

0.25

0.2

0.15

0.1

0.05

0

0

0.05

0.1

0.15

Absorption Coefficient

0.2

0

0.05

0.1

0.15

Absorption Coefficient

0.2

0

0.05

0.1

0.15

Absorption Coefficient

0.2

0

0.05

0.1

0.15

0.2

Absorption Coefficient

Figure 4. Agreement between simulated and measured data for a hallway with fewer diffusing surfaces.

6019