Airborne sound transmission modeling of walls
based on random point process theory
Cédric VAN HOORICKX and Edwin P.B. REYNDERSKU Leuven, Department of Civil Engineering, Belgium, cedric.vanhoorickx@kuleuven.be

Abstract
A novel approach is presented to compute the diffuse sound transmission through modally sparse walls. The
fact that the natural frequencies in the adjacent rooms can be considered to form a sequence of random points on
the frequency axis makes this problem amenable to analysis by random point process theory. The method makes
use of the fact that for any diffuse subsystem the local spacings between the eigenvalues saturate to those of a
Gaussian orthogonal ensemble (GOE) matrix and the statistics of its mode shapes saturate to Gaussian randomﬁelds, with a spatial covariance function depending only on wavelength and distance. This stochastic modelling
approach is computationally very efﬁcient and is illustrated for a calcium silicate block wall and a gypsum blockwall.
Keywords: Random point process theory, Airborne sound transmission, Diffuse ﬁeld

1

INTRODUCTION

In a measurement setup, the airborne sound insulation of a partitioning wall is subject to many sources of
uncertainty. These are often difﬁcult to quantify due to lack of information or model complexity. At highfrequencies, the sound pressure at a speciﬁc location in the room is sensitive to local variations of the acoustic
ﬁeld properties caused by random wave scattering inside and at the boundaries of the source and receiverroom. Usually a diffuse ﬁeld approach is considered, which is a random ﬁeld, composed of a large number ofstatistically independent plane waves with a uniformly distributed spatial phase.
When the local sound pressure in the room is highly sensitive to local variations of the acoustic ﬁeld properties,
it can be modeled as a pure-tone diffuse reverberant ﬁeld, as in Statistical Energy Analysis (SEA). In thisapproach, the system is decomposed into homogeneous subsystems whose response is characterized with asingle random variable, namely its total energy. Next to the mean of the total sound energy over the random
ensemble of rooms, also the variance can be predicted [4], employing a non-parametric model of uncertainty.
In many applications, some of the subsystems of a built-up system are very sensitive to small wave scatterers,
while others are not. In this mid-frequency range, a hybrid method can be used combining deterministic andstochastic energy based methods. This has been achieved in the hybrid FE-SEA approach [10], based on areciprocity relationship that links the mean vibrational energy of a non-parametric random SEA subsystem tothe variance of the nodal forces at its boundary.
This paper presents a novel approach to compute the airborne sound insulation of walls. The same assumptions
are made as for a hybrid FE-SEA approach, namely that (1) the wave ﬁeld of each random subsystem acrossthe random ensemble can be assumed diffuse and (2) the mode shapes saturate to Gaussian ﬁelds. Instead of
directly computing the ensemble mean and variance of the total energy, the approach presented here uses a
Monte Carlo simulation, where every sample represents a realization of the diffuse subsystem. An advantageof this procedure is that parametric uncertainty in the modally sparse subsystems can be easily included byalso varying the uncertain parameters in every Monte Carlo realization. This is useful in e.g. assessing the
reproducibility of a measurement procedure using a virtual round robin test [12]. The approach utilizes therandom point process theory, which relies on the assumption that the natural frequencies can be considered toform a sequence of random points on the frequency axis. This assumption is only valid the natural frequencies
mix with respect to each other across the random ensemble of rooms, which is the case at high frequencies.

1192

2

PREDICTION MODEL FOR AIRBORNE SOUND INSULATION

The problem considered in this paper is the airborne sound transmission through a modally sparse wall inbetween two diffuse rooms. Quantities related to the source room are indicated with subscript 1, quantities
related to the wall with subscript 2, and quantities related to the receiver room with subscript 3.
2.1 System of equations
The transmission loss is computed using the assumed-modes method, approximating the pressure ﬁelds of thesource room p1 and reveiver room p3 and the vibration ﬁeld of the wall u2 into a ﬁnite set of basis functions:n1/3

p1/3 ≈ ∑ φ1/3i (x, y, z)q1/3i (ω) = φ 1/3 q1/3
i=1n2

u2 ≈ ∑ φ2 j (y, z)q2 j (ω) = φ 2 q2

(1)

j=1

Inserting these approximations into Lagrange’s equations of motion and adopting a hysteretic damping modelyields the following linear system of equations:⎤⎧ ⎫ ⎧ ⎫⎡0
D11 K12⎨q1 ⎬ ⎨f1 ⎬⎣ K21 D22 K23 ⎦ q2 = 0(2)⎩ ⎭ ⎩ ⎭
0
K32 D33
q3
0
where

Dkk = −ω 2 Ik + Ω 2k (1 + iηk )(3)with i the imaginary unit and Ω k the diagonal matrix containing the circular undamped eigenfrequencies forsubsystem k corresponding to the mode shapes in φ k . The damping loss factors of the rooms are computed fromtheir reverberation times via η = 4.4π/(ωT ). A point monopole with volume acceleration ap (ω) is assumed,
located at (xp , yp , zp ). Element i of the loading vector f1 then reads
f1i (ω) = −ρa ap (ω)φ1i (xp , yp , zp )(4)

The matrices K12 and K32 are coupling matrices that represent the loading on the room due to the plate movement, and K21 and K23 represent the loading on the plate due to the room pressure. The elements of thecoupling matrices are:
K21, ji =K23, jk =

Ly2

Lz2

0
Ly2

0
Lz2

0

0

φ1i (Lx1 , y, z)φ2 j (y, z)dzdy

K12,i j = −ρa ω 2 K21, ji(5)

φ3k (Lx1 , y, z)φ2 j (y, z)dzdy

K32,k j = −ρa ω 2 K23, jk(6)

2.2 Computation of the sound reduction index
In order to efﬁciently solve equation (2), row reduction is performed on the block matrices to reduce the sizeof the system of equations in equation (2). This yields:
Aq2 = b(7)

where the elements of the matrices A and b are given by:
n3
K23,mk K23,nk
K21,mi K21,ni+ ρω 2 ∑D11,iiD33,kki=1
k=1
n1

Amn = D22,mn + ρω 2 ∑n1

K21,mi f1i
i=1 D11,ii

bm = − ∑

(8)

1193

The energy in the source and receiver room can then be computed as:
2

n2

E1 =

1 n1
∑
ρa i=1j=1|D11,ii |

2

2

n2

f1i − ∑ K12,i j q2 j≈

1 n1 | f1i |2∑ |D |2ρa i=111,ii

n3

E3 = ρa ω 4 ∑k=1∑ K23, jk q2 j

j=1

D33,kk2

The sound reduction index R is ﬁnally computed with the measurement formula:
 
 
 
E1
E3
S2
R = 10 log10− 10 log10+ 10 log10
V1
V3
A3

(9)

(10)

in which V1 and V3 are the volume of the source and receiver room, S2 is the surface area of the partition, andA3 is the absorption of the receiver room.

3

SUBSYSTEM PROPERTIESThe computations require the natural frequencies and mode shapes of the subsystems to be computed. The
rooms are modeled as diffuse subsystems, the wall is assumed to be a deterministic subsystem.
3.1 Deterministic subsystem
If the wall can be modeled as a homogeneous, isotropic, thin plate, its natural frequencies are given by:


 
 
m jy π 2
m jz π 2
D2
ω2 j =
+
m2
Ly2
Lz2

(11)

with D2 the bending stiffness, m2 the mass of the wall per unit area, Ly2 the width, and Lz2 the height of thewall. The integers m jy and n jz represent the number of half wavelengths in the plate dimensions. The modeshapes of the wall are calculated as follows:

 


 

m jy πym jz πz2
φ2 j =  
sin
(12)
sin
= A2 sin k2 jy y sin k2 jz z
Ly2
Lz2
m2 Ly2 Lz2
3.2 Diffuse subsystems
The wave ﬁelds in the rooms are modeled as diffuse ﬁelds. Diffuse ﬁeld models essentially describe rooms
with locally uncertain geometric properties, material properties, and/or boundary conditions, that is, rooms withuncertain wave scatterers. The mode shapes of a room in high-frequency regime can be interpreted as standing
waves that arise from many traveling plane wave components. Adopting a diffuse acoustic ﬁeld model, thepressure mode shape at a given location consists of a summation of independent plane acoustic waves withthe same mean amplitude and uncorrelated phases, coming from all directions with equal probability. It thenfollows from the central limit theorem that the acoustic mode shapes are zero-mean, Gaussian random ﬁelds. A
zero-mean Gaussian random ﬁeld is uniquely determined by its covariance function.
For diffuse reﬂecting boundaries, the mode shapes φsi are statistically homogenous, i.e., the statistics of thepressure mode shape components are independent of their position. The corresponding random wave ﬁeld is adiffuse ﬁeld. For three-dimensional volumes, its covariance function then has the form [1]: As j0 (ksi |x − x |)for 2 points inside the room, away from a deterministic boundary
Csi x, x =(13)

2As j0 (ksi |x − x |) for 2 points on a reﬂecting boundary, away from joints and corners
where j0 (x) = sin(x)/x is the spherical Bessel function of the ﬁrst kind and order zero, ksi := λ2πsi denotes thewavenumber corresponding to the wavelength λsi of mode i for subsystem s, and As is a factor that is independent of position, which can be determined from the mode shape normalization condition. In an acoustic1194

enclosure Ω with volume Vs , the normalization condition reads:
1 2
φ (x)dx = 1
2 si
Ωc

⇔

As =

c2
.
Vs

(14)

The correlation function (13) depends only on the distance between the considered mode shape components,
the wavelength and the total volume. Close to perfectly reﬂecting boundaries, the mean squared sound pressure
equals twice the mean squared sound pressure in the center of the room [13]. Because of this, a factor two
appears in the covariance function (13) for points located on the reﬂecting structural element. In addition, when
the distance between x and x is large compared to the wavelength λsi , it follows from the above expressions
that the mode shapes evaluated at these distinct points are approximately uncorrelated.
3.3 Computation of the coupling matrices
Instead of numerically evaluating the integrals in equations (5) and (6), the statistics of the coupling matrices
are computed directly. When integrating a Gaussian ﬁeld, it is the limit of a linear combination of Gaussian
random variables so it is again Gaussian. The entries of the coupling matrix are therefore Gaussian variables
and are determined by their mean and variance. As the mean of the mode shapes in the rooms equals zero,also the mean of the coupling matrices equals zero. The variance of K21,kl is computed as follows:
Var K21, ji =E φ1i (Lx1 , y, z)φ1i (Lx1 , y , z ) φ2 j (y, z)φ2 j (y , z )dz dzdy dy(15)

The term E [φ1i (Lx1 , y, z)φ1i (Lx1 , y , z )] is the covariance function for two points on a reﬂecting boundary (equation (13)), and this integral becomes:

Var K21, ji = 2A1 A22


 

 
 
 
j0 k1i |x − x | sin k2 jy y sin k2 jy y sin k2 jz z sin k2 jz z dz dzdy dy(16)where the integrals go from 0 to Ly2 for the ones in y and y and from 0to Lz2 for the ones in z and z . The
distance function in the spherical Bessel function is given by |x − x | = (y − y )2 + (z − z )2 . The quadruple
integral can be converted into the following double integral:
2 2
Var K21, ji = 2A1 A22 Ly2
Lz2

1 1

 
 
 

j0 a1iy u2 + γ 2 v2 h a2 jy , u h a2 jz , v dvdu

(17)

0 0

with a1iy = k1i Ly2 , γ = Lz2 /Ly2 , a2 jy = k2 jy Ly2 , a2 jy = k2 jz Lz2 , and:h (a, u) = (1 − u) cos (au) −cos (a)sin (a(1 − u))a(18)

The integral in equation (17) is evaluated numerically. The mean and variance of K23, jk is computed in thesame way.

4

RANDOM POINT PROCESS THEORY

4.1 Introduction
Note that for obtaining the energies in the rooms, sums are needed over expressions with D11,kk or D33,kk in the
denominator, i.e. in Eqs. (8) and (9). At high frequencies, the computation cost of the sound reduction index isdominated by these sums. The natural frequencies can be considered to form a sequence of random points on
the frequency axis, and thus the problem is amenable to analysis by random point process theory [11]. Strictly,
the natural frequency point process is assumed to be stationary, which is only the case if the modal density of1195

the subsystems is constant, and this is not true in this case. On the other hand, only a limited number of modes
signiﬁcantly contribute at any considered frequency. As the modal density is approximately constant over thefrequency range corresponding to these modes, it is reasonable to assume that the point process is stationary.
A distinction is made between type 1 and type 2 sums, with type 1 referring to the sums in Eq. (8) and type2 referring to the ones in Eq. (9). The terms of these sums have D11,kk or D33,kk in the denominator. From thedeﬁnition of the matrices D11,kk and D33,kk , the type 1 sums S1 (ω) and type 2 sums S2 (ω) can be expressed
as:
ck
ck
ck
=∑ 2
(19)
S1 (ω) = ∑S2 (ω) = ∑2
2 + ω 2 (1 + iη)2 )2 + η 2 ω 4
2
2
−ω
(ω
−
ω
k
k −ω + ωk (1 + iη)kkkk
Unless ck is much smaller at ω ≈ ωk than at other frequencies, it follows that the modes mainly contributing to
the summations are the ones in the vicinity of the considered frequency. These sums can also be written as:∞∞

h1 (ω, ω  )ξ (ω  )dω 

S1 (ω) =h2 (ω, ω  )ξ (ω  )dω 

S2 (ω) =

−∞

(20)

−∞

where

ξ (ω  ) = ∑ ck δ (ω  − ωk )

(21)

k

h1 (ω, ω  ) =h2 (ω, ω  ) =1
ω  − ω − iωη/21≈−ω 2 + ω 2 (1 + iη) 2ω (ω  − ω)2 + (ηω/2)2

(22)

1
1
4ω 2 (ω  − ω)2 + (ηω/2)2(23)

1(ω 2 − ω 2 )2 + η 2 ω 4≈

in which the approximations are commonly used [3, 5].
Random point process theory now allows computing the statistics of these sums. These are separately discussed
for the two sum types in the next to subsections.
4.2 Statistics of type 1 sums
Weaver [14] found that, for the generic case where the considered spatial uncertainty does not preserve symmetries, the statistics of the local eigenvalue spacings saturate to those of the Gaussian Orthogonal Ensemble(GOE) matrix from random matrix theory [7]. It can be shown that the GOE assumption leads to the following
statistics (more precisely the ensemble mean and variance) of the type 1 sum S1 (ω) as deﬁned in (20) [5]:


πn
1
2
E
[c
]
α
−
1
+
(24)
k
4ω 3 η
πm
 where n is the modal density, m = ωηn the modal overlap factor and α = E c2k /E [ck ]2 .
As for diffuse subsystems, these sums are performed over many modes, it is assumed, according to the central
limit theory, that the real and imaginary part of the sums have a Gaussian distribution. This is for example
consistent with the assumption of Schroeder [9] that transfer functions written in the form of a modal sum,which are closely related to the type 1 sums1 , are approximately complex Gaussian. A realization of a type 1
sum is therefore obtained from:πniE [ck ] πn12
+ (ξR + iξI )E [ck ] α − 1 +
S1 (ω) = −(25)2ω
4ω 3 ηπm
E [S1 (ω)] = −iE [ck ] πn2ωVar (Re {S1 (ω)}) = Var (Im {S1 (ω)}) =1 The similarity between the type 1 sums and the transfer function in modal form can be easily seen by comparing Eq. (19) and the transfer
function expression e.g. in [5]: the denominator of the sum terms is the same, the numerator consists of the product of two zero-mean Gaussian
factors1196

with ξR and ξI realizations of normally distributed random numbers with mean 0 and standard deviation 1.
The only terms still to be determined are the mean and mean of the square of the factors ck for the sums in (8).
These are calculated from basic probability rules:E [K21,mi K21,ni ] = σK221,m (ω)δmn
E (K21,mi K21,ni )2 = (2δmn + 1)σK221,m (ω)σK221,n (ω)(26)
ρ 2 c2 2
E (K21,mi f1i )2 = σK221,m (ω)σ 2f1i =
σ
(ω)
(27)
E [K21,mi f1i ] = 0
V1 K21,mω in which σK221,m (ω) = Var K21,mi k1i =.
c
4.3 Statistics of type 2 sums
The mean and variance of type 2 sums equal [2, 3]:E [ck ] πn
E [S2 (ω)] =2ηω 3

Var (S2 (ω)) = E [S2 (ω)]

2

1
α −1
+
πm
(πm)2


(28)

The type 2 sums appear in Eq. (9), which represent the total energies in the source and receiver room. As thetotal energy of a subsystem is empirically found to follow a log-normal distribution [6], also the sums appearing
in Eq. (9) are assumed to have a log-normal distribution. A realization of a type 2 sum is therefore obtained
as follows. First the mean μ2 and standard deviation σ2 of the natural logarithm are computed:
⎞
⎛


⎟
⎜
⎟
⎜
E
[S
(ω)]
(ω))
Var
(S
2
2
⎟
σ22 (ω) = ln 1 +μ2 (ω) = ln ⎜
(29)
⎟
⎜
E [S2 (ω)]2
Var (S2 (ω)) ⎠⎝1+
E [S2 (ω)]2after which a realization is obtained fromS2 (ω) = exp (μ2 (ω) + ξ σ2 (ω))(30)

with ξ a realization of a normally distributed random number with mean 0 and standard deviation 1.
For the sums in Eq. (9), the mean and mean of square of the factors ck are given by:  ρ 2 c2
E | f1i |2 = σ 2f1i E ξ1i2 =
V1
⎤
⎡
n2

2



 
ρ 4 c4E | f1i |4 = σ 4f1i E ξ1i4 = 3 2
V1

(31)

n2

2
E ⎣ ∑ K23, jk q2 j ⎦ = ∑ σK223, j (ω) q2 j
j=1

⎡

j=1
n2

4

⎤

E ⎣ ∑ K23, jk q2 j ⎦ =
j=1



2

∑ q2 jj

2

σK223, j (ω)&2%+ 2 ∑ ∑ Re q2i q∗2 j σK223,i (ω)σK223, j (ω)
i

(32)

j

In summary, the ensemble probability distribution of the sound reduction index R is obtained as follows. A
Monte Carlo simulation is performed in which each sample represents a realization of the system. In eachsample, realizations of the sums in Eq. (8) are computed using Eq. (25), after which the system of equations in
Eq. (7) is solved to obtain the vibration ﬁeld of the wall q2 . The energies in the source and receiver room arethen computed from Eq. (9), in which realizations of the sums are computed from Eq. (30). Finally, the soundreduction index is computed with the measurement formula in Eq. (10).

1197

5

APPLICATIONSIn this section, the proposed method is applied to compute the airborne sound insulation of a calcium silicate
block wall and a gypsum block wall. A total of 1000 Monte Carlo realizations are computed. The approachfrom this paper is compared with two other approaches: a deterministic model, where the sound ﬁeld in theemitting and receiver rooms are expanded as a sum of their hard-walled modes, and the hybrid FE-SEA approach as presented in [8]. In all computations, the source and receiver rooms have a volume of 87 m3 and areverberation time of 1.5 s. The wall size is 3 × 3.3m2 .
First, a heavy single-leaf wall is studied. The wall has a thickness of 25 cm and consists of calcium silicate
blocks with mass density ρ = 1800 kg/m3 , Young’s modulus E = 10.8 GPa, and Poisson’s ratio ν = 0.2. The
mean of the predicted sound transmission loss and the corresponding 95 % conﬁdence interval are shown in
Fig. 1a. The dips in the insertion loss correspond to the individual wall modes, indicating that, especially at
low frequencies, the modal density of the wall is generally low (n = 0.0274 modes/Hz). The critical frequencyof the wall is at 102.5 Hz, yet no clear coincidence dip is observed because of the low modal density. Above thecritical frequency, the transmission loss increases with approximately 9 dB per octave. The presented approachcorresponds almost perfectly corresponds with the hybrid FE-SEA approach.
80

60
50

60

R [dB]

R [dB]

40
40

30
20

20
10
0

(a)

63

125

250

500

Frequency [Hz]

1000

0

2000

(b)

63

125

250

500

1000

2000

Frequency [Hz]Figure 1. Mean and 95 % interval of the sound reduction index for (a) the calcium silicate wall and (b) the
gypsum block wall computed using three approaches: deterministic approach (blue dash-dotted line), hybrid
FE-SEA approach (red solid line), and the approach presented in this paper (black dashed line).
The second example is a lightweight single-leaf wall with a thickness of 10 cm. It consists of gypsum blocks
with mass density ρ = 910 kg/m3 , Young’s modulus E = 3.15 GPa, and Poisson’s ratio ν = 0.2. The mean of
the predicted sound transmission loss and the corresponding 95 % conﬁdence interval are shown in Fig. 1b. In
contrast to the heavy wall from the previous section, no clear dips are observed in the harmonic transmission
loss statistics, as the modal density is higher (n = 0.0899 modes/Hz). A clear coincidence dip is observed
around the critical frequency of 334.2 Hz, behind which the transmission loss increases with approximately 9 dB
per octave. Also for this case, a good correspondence with the hybrid FE-SEA approach is obtained, especiallyat high frequencies.

6

CONCLUSIONSIn this paper, a novel approach is presented to compute the airborne sound transmission through a wall locatedbetween two diffuse rooms. The approach is therefore a hybrid approach, coupling the deterministic structuralsubsystem with the random acoustic subsystems. The computation of the airborne sound transmission is performed through a Monte Carlo simulation, in which each sample corresponds to a realization of the random
ensemble. The natural frequencies in the diffuse rooms are considered to form a sequence of random points

1198

on the frequency axis, which makes this problem amenable to analysis by random point process theory. The
approach furthermore uses the fact that for any diffuse subsystem the mode shapes saturate to Gaussian randomﬁelds, with a covariance function depending only on wavelength and distance. Two applications are presented,
a calcium silicate block wall and a gypsum block wall. For both applications, the presented approach shows
good correspondence with the hybrid FE-SEA approach.

ACKNOWLEDGEMENTS
The research presented in this paper has been performed within the frame of the VirBAcous project (project
ID 714591) “Virtual building acoustics: a robust and efﬁcient analysis and optimization framework for noisetransmission reduction” funded by the European Research Council in the form of an ERC Starting Grant. The
ﬁnancial support is gratefully acknowledged.

REFERENCES[1] R. Cook, R. Waterhouse, R. Berendt, S. Edelman, and M. Thompson Jr. Measurement of correlation
coefﬁcients in reverberant sound ﬁelds. Journal of the Acoustical Society of America, 27(6):1072–1077,
1955.
[2] R. Langley and A. Brown. The ensemble statistics of the band-averaged energy of a random system.
Journal of Sound and Vibration, 275(3–5):847–857, 2004.
[3] R. Langley and A. Brown. The ensemble statistics of the energy of a random system subjected to harmonic
excitation. Journal of Sound and Vibration, 275(3–5):823–846, 2004.
[4] R. Langley and V. Cotoni. Response variance prediction in the statistical energy analysis of built-up systems. Journal of the Acoustical Society of America, 115(2):706–718, 2004.
[5] R. Langley and V. Cotoni. The ensemble statistics of the vibrational energy density of a random system
subjected to single point harmonic excitation. Journal of the Acoustical Society of America, 118(5):3064–
3076, 2005.
[6] R. Langley, J. Legault, J. Woodhouse, and E. Reynders. On the applicability of the lognormal distribution
in random dynamical systems. Journal of Sound and Vibration, 332(13):3289–3302, 2013.
[7] M. Mehta. Random Matrices. Elsevier, San Diego, CA, 3rd edition, 2004.
[8] E. Reynders, R. Langley, A. Dijckmans, and G. Vermeir. A hybrid ﬁnite element - statistical energy
analysis approach to robust sound transmission modelling. Journal of Sound and Vibration, 333(19):4621–
4636, 2014.
[9] M. Schroeder. Statistical parameters of the frequency response curves of large rooms. Journal of the AudioEngineering Society, 35(5):299–306, 1987.
[10] P. Shorter and R. Langley. Vibro-acoustic analysis of complex systems. Journal of Sound and Vibration,288(3):669–699, 2005.
[11] R. Stratonovich. Topics in the theory of random noise, volume I. Gordon and Breach, New York, NY, 1963.
[12] C. Van hoorickx and E. Reynders. Uncertainty quantiﬁcation of sound transmission measurement procedures based on the Gaussian Orthogonal Ensemble. In D. Herrin, J. Cuschieri, and G. Ebbitt, editors,
Proceedings of the 47th International Congress and Exposition on Noise Control Engineering, Inter-Noise 2018,
Chicago, USA, August 2018. CD-ROM.
[13] R. Waterhouse. Interference patterns in reverberant sound ﬁelds. Journal of the Acoustical Society of America,
27(2):247–258, 1955.
[14] R. Weaver. On the ensemble variance of reverberation room transmission functions, the effect of spectral
rigidity. Journal of Sound and Vibration, 130(3):487–491, 1989.

1199