Click here to Skip to main content
15,878,959 members
Articles / Programming Languages / C#

Propagation of Sound in a Room

Rate me:
Please Sign up or sign in to vote.
5.00/5 (11 votes)
10 Mar 2023CPOL9 min read 15.6K   277   23   10
Calculate sound transmission between a point source and a receiver point within an enclosed space
This article will go through the steps and limitations of calculating the sound transmission between two points in a room by using modal analysis. This will use Fourier transformation in order to do the calculations based on the wave equation and fully reflective walls. In addition, it will also give two different models for adding dampening to the room, one uses reverberation time and the other will use the absorption coefficient.

Image 1


Sound transmission between points within a room is one of the most important aspects in practical acoustical problems. In many aspects, it is the cornerstone of what we do, since so many things are dependent on the sound pressure transmission from a source to its listeners. The relation calculated by the geometrical damping and the damping at the boundaries (or wall, ceiling and floor) given a sound source with a reference sound pressure that spreads the sound equally in all directions. The diagram shown is the damping of a flat response of 0 at the source and therefore the negative pressure at the receiver.

The general procedure for calculating the sound pressure level of a receiver have two main approaches: geometrical acoustics which is only valid for high frequency and Fourier acoustics generally only used for the low frequency part due to the number of nodes needed in the calculation increase drastically for increasing frequency.

Wave equation

We start off with the classical time and space depended wave equation:

$ \frac{d^2 \, f(t,x)}{dx^2} = \frac{1}{c^2} \frac{d^2 \, f(t,x)}{dt^2} \tag{1} \label{eq:Wave}$

This equation describes all the details of the sound transmission between all possible sources and receivers, but it is very time consuming to do the calculations. There is a simple way assuming that you do not need all details of the transmission. We perform Fourier (or Laplace) transformation on the differential equation itself and get Helmholtz equation instead:

$ \nabla^2 f(x) + k^2 \cdot f(x) = 0 \tag{2} \label{eq:Helmholtz}$

This is what acousticians call a steady state solution to the wave equation, since you integrate over all time to eliminate the time dependency in the equation. For a mathematician, this can lead to another problem: Gibbs phenomena. The Fourier transformation of the wave equation actually eliminates some disjoint initial conditions and solutions that it cannot correct. Actually, the history of Fourier transform begins in 1750 with Euler when he considered the initial condition which is the Fourier transformation of the curve. An English translation of Euler's original work is given here: PDF version of an English translation and there are more translations of more of Euler's work that are also given here. In the article in question, he wrote the initial condition as:

$ \alpha \sin(\pi s/a) + \beta \sin(2 \pi s/a) + \gamma \sin(3 \pi s/a) \tag{3} \label{eq:EulerIC}$

Which is exactly a summation of the Fourier coefficients for as also indicated by Daniel Bernoulli for suspended strings, which he published in 1753. Euler had some reservations as he didn't think that the Fourier transformation was general enough, and he is actually correct in his analysis. The Weierstrauss theorem will actually not fix this. In fact, the only reason that Fourier transformation is valid, is an article by Carleson that finally set the question on firm ground. The formulated proof is known as Carleson's theorem states that Fourier transform converges almost everywhere. The issue is that there must be a finite number of disjoint steps in the function, then it will converge. And it will converge almost everywhere except for the disjoint step. This can be illustrated by a step function in blue and the Fourier transformation approximation in red shown below:

Image 2

What you do to reproduce the error is to give a step function as the initial condition, take the Fourier transformation of it, and at last the inverse Fourier transformation.

The reason for the discrepancy in the figure is simply that the Fourier transform is infinitely smooth, meaning that infinitely many derivatives exist in every point. While the step functions it tries to emulate is missing the derivative at the step point. The reality is that the Fourier transform can accurately replicate any function in \(L^2(f(x)) \), which is a fancy way of saying that on average it will be correct, but that does not mean that the approximation from the Fourier transform are equal in value at every point.

Solving Helmholtz equation

We now assume that all the boundaries of the room are infinitely hard. This will mean that the particle velocity at the wall must be zero, which in return will mean that the sound pressure will be the maximum at the same point. Now we are searching for a solution that has its maximum value at the boundaries, so this would be achieved by harmonics of the cos function which is also a solution to the homogenous wave equation.

$ p(x,y,z) = C \cdot \cos\left(\frac{M \pi x}{V_x}\right) \cdot \cos\left(\frac{L \pi y}{V_y}\right) \cdot \cos\left(\frac{K \pi z}{V_z}\right) \tag{4} \label{eq:HelmholtzBC}$

Assuming that we solve Helmholtz equation in each direction using a power series method. We generally want to add a source to the differential equation:

$ \nabla^2 f(x) + k^2 \cdot f(x) = - i \cdot \omega \cdot \rho_0 q(r_0) \tag{5} \label{eq:HelmholtzInhomogenous}$

This is a inhomogeneous Helmholtz equation that can be solved using a power series method. The source distribution function. This can be expressed as individual frequency amplitudes:

$q(r_0) \cdot e^{i \omega t} = \sum_{N=0}^{\infty} Q_N \Psi_N(r) e^{i \omega t} \tag{6} \label{eq:InsertedEquation}$

One also needs the property of the form function to follow certain laws. Mainly that the orthogonality property will give:

$ V = \iiint\limits_{V} \Psi^2_N (r) \, dV $

This will require some shape property and can be illustrated by each of the energy integrals from 1,2 and 3 nodes contributing to the sound level.

$ w_{axial} = \frac{1}{L_{x}} \int_0^{L_x}{\cos^2{ \left( x \frac{L \cdot \pi}{L_x} \right) \, dx}} = \frac{1}{2}$
$ w_{tangential} = \frac{1}{L_{x} L_y} \int_0^{L_x} \int_0^{L_y} \cos^2{\left( x \frac{L \cdot \pi}{L_x} \right)} \cos^2{\left( y \frac{M \cdot \pi}{L_x} \right)} \, dx dy = \frac{1}{4}$
$ w_{oblique} = \frac{1}{L_x L_y L_z} \int_0^{L_x} \int_0^{L_y} \int_0^{L_z} \cos^2{\left( x \frac{L \cdot \pi}{L_x} \right)} \cos^2{\left( y \frac{M \cdot \pi}{L_y} \right)} \cos^2{\left( z \frac{K \cdot \pi}{L_z} \right)} \, dx dy dz = \frac{1}{8}$

This means that the axial mode will have double the energy than an tangential node, and that a tangential node will have twice the energy of the oblique node. However the same integrals will appear with the orthogonal requirement, but there we demand that each of these three integrals yields the same result. Therefore, each time either L, M or K is zero, we will have to compensate with 2 for each mode that is zero and 1 for each mode that are not zero.

Applying the power series to every element in the differential equation gives you the following resulting equation:

$ \sum_{N=0}^\infty A_N \nabla \Psi_N(r) + \sum_{N=0}^\infty A_N k^2 \Psi_N(r) = -i \omega \rho_0 \sum_{N=0}^\infty Q_N \Psi_N(r) \tag{7} \label{eq:TotalHH}$

Rearranging the terms:

$ \sum_{N=0}^\infty A_N \Psi_N(r) \left[ k^2 - k_N^2 \right] = -i \omega \rho_0 \sum_{N=0}^\infty Q_N \Psi_N(r) \tag{8} \label{eq:TotalClean}$

Setting up the final solution with regards to the pressure at the receiving point:

$p(x,y,z) = \rho_0 c_0^2 Q \sum_{L=0}^\infty \sum_{M = 0}^\infty \sum_{K = 0}^\infty \frac{\omega \cdot \Psi_{L,M,K}(x,y,z) \cdot \Psi_{L,M,K}(x_0,y_0,z_0)}{V_{L,M,K} \cdot \left( \omega^2 - \omega^2_{L,M,K} \right)} \tag{9} \label{eq:Final}$

Include damping

In room acoustic the damping of the room from soft surfaces gives a substantial effect on the sound level in a room.

$m \cdot \frac{d^2 x}{dt^2} + k (1 + i \cdot \eta) \cdot x = 0 \tag{10} \label{eq:PendulumWave}$

Assuming that the solution can be written with the complex \( s = \sigma + i \tau \) on the form \( e^{i s \cdot t} \). Inserting this expression into the previous differential equation will give the following

$s = \sqrt{\frac{k}{m} (1 + i \eta)} \approx \omega_0 \left( 1 + i \frac{\eta}{2} \right) \tag{11}\label{eq:ComplexSolutionWave}$

We now assume that the energy of any sound waves will be negligible after it has decayed by 60 dB. This is a historical picked value as speech typically have a level at 60 dBA so this was assumed to be the level which you could no longer hear the sound \( \eqref{eq:ComplexSolutionWave} \) . This was suggested in the 1900's before the advent of electrical equipment that actually made it possible to measure the actual sound level, so this logic was assumed to be correct. Now it is just a convenient way of expressing the damping value of sound and vibrations in a confined space, and the resulting damping value is dB per time unit. In practice signal to noise ration determines the usable range, and this is seldom more than 20 - 30 dB in room acoustics. Its used in a general way for acousticians to describe damping in both rooms and in plates, wall etc.

$\eta = \frac{\log \left(10^{6} \right) }{T_{60} \cdot \omega} = \frac{\log\left(10^{6} \right)}{T_{60} \cdot 2 \pi f} = \frac{2.2}{T_{60} \cdot f} \tag{12} \label{eq:Damping}$

All that is left is to gather this information to the reverberation time or the absorption factor of the surfaces in the room. \( \eqref{eq:Damping} \)

$ \omega_{L,M,K} = \omega_{L,M,K} \cdot \sqrt{1+ i \cdot \eta } = \omega_{L,M,K} \cdot \sqrt{1 + i \frac{4.4 \cdot \pi}{\omega_{L,M,K} \cdot T_{60}}} \tag{13} \label{eq:TotalComplete}$

We now insert this damping into the previous equation:

$p(x,y,z) = \rho_0 c_0^2 Q \sum_{L=0}^\infty \sum_{M = 0}^\infty \sum_{K = 0}^\infty \frac{\omega \cdot \Psi_{L,M,K}(x,y,z) \cdot \Psi_{L,M,K}(x_0,y_0,z_0)}{V_{L,M,K} \cdot \left( \omega^2 - \omega^2_{L,M,K} - \frac{i \cdot \omega_{L,M,K} \cdot 4.4 \cdot \pi}{T_{60}} \right)} \tag{14} \label{eq:TotalCompleteReverberation}$

There is another way of estimating the damping term by simply inserting the experimentally found variable of damping \( \eta = \frac{\alpha}{8} \) by using experience from measurements and together with some approximation techniques.

$p(x,y,z) = \rho_0 c_0^2 Q \sum_{L=0}^\infty \sum_{M = 0}^\infty \sum_{K = 0}^\infty \frac{\omega \cdot \Psi_{L,M,K}(x,y,z) \cdot \Psi_{L,M,K}(x_0,y_0,z_0)}{V_{L,M,K} \cdot \left( \omega^2 - \omega^2_{L,M,K} - i 2 \omega_{L,M,K} \left( \epsilon_L \frac{\alpha_{x}}{8 L_x} + \epsilon_M \frac{\alpha_{y}}{8 L_y} + \epsilon_K \frac{\alpha_{z}}{8 L_z} \right) \right)} \tag{15} \label{eq:TotalCompleteWithAlpha}$

This concludes this article. Plans for extending the solutions for boundaries with impedances is on its way, unless you can beat me to it.



  • 23rd February, 2023: Initial version



This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Written By
Chief Technology Officer
Norway Norway
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

QuestionTough Room Pin
Member 1488449217-Mar-23 12:27
Member 1488449217-Mar-23 12:27 
AnswerRe: Tough Room Pin
Kenneth Haugland18-Mar-23 6:00
mvaKenneth Haugland18-Mar-23 6:00 
GeneralRe: Tough Room Pin
Member 1488449220-Mar-23 11:27
Member 1488449220-Mar-23 11:27 
GeneralRe: Tough Room Pin
Kenneth Haugland21-Mar-23 6:58
mvaKenneth Haugland21-Mar-23 6:58 
GeneralRe: Tough Room Pin
Member 1488449222-Mar-23 8:00
Member 1488449222-Mar-23 8:00 
SuggestionToo dirty Pin
Sergey Alexandrovich Kryukov23-Feb-23 17:11
mvaSergey Alexandrovich Kryukov23-Feb-23 17:11 
GeneralRe: Too dirty Pin
Jeffamn24-Mar-23 3:55
Jeffamn24-Mar-23 3:55 
AnswerRe: Too dirty Pin
Sergey Alexandrovich Kryukov24-Mar-23 4:32
mvaSergey Alexandrovich Kryukov24-Mar-23 4:32 
Questionstrange program Pin
Member 1531459123-Feb-23 8:50
Member 1531459123-Feb-23 8:50 
AnswerRe: strange program Pin
Kenneth Haugland23-Feb-23 8:56
mvaKenneth Haugland23-Feb-23 8:56 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.