In this episode we will look at how to make a simple Julia model of one of the simplest systems in acoustics, a rectangular room with rigid walls, assuming adiabatic wave propagation. Even if this system is among the simplest in acoustics it is actually already very complicated. As such, we will focus only on the modes, one part of the problem, without attempting impulse response simulation or other fancy things like that, for now.

Motivation

We use Julia to implement the model by taking advantage of the fact that a closed form solution exists for this problem. So, we will not use a solver. This might feel strange, as in the previous episodes we covered the used of solvers to solve complex equations, with a special focus on FEM. What is the reason to do this?

The reason is that these simple systems constitute very good benchmarks for our solvers. As already discussed in The FEM Pipeline episode, solvers can only provide us with approximate solutions, and the quality of the approximation will depend on how well we set the solver parameters and input data. But if we don’t know the solution, how can we know how good an approximation we got? This is where systems with known solutions come into place: we can directly compare solvers outputs with the known analytical solution, thus being able to assess how much the solver is being accurate and reliable. Once we figure out how to configure our solver thanks to the benchmark model, we can then move to more complex and realistic problems.

We will use what we find in this episode in future episodes.

The Equation

We will very briefly look at the equation. A comprehensive solution and analysis is out of the scope of the episode, but we will look at the most important parts so to clarify what we are doing. If you are an expert, just skip the section. If you instead want to know more, Fundamentals of Acoustics would be a good reference.

The equation governing adiabatic wave propagation in a compressible fluid (such as air) is the so called linear wave equation (which was already introduced in the What is Acoustic Modelling episode):

$$\nabla^{2}p=\frac{1}{c^{2}}\frac{\partial^{2}p}{\partial t^{2}}$$

where $c$ is the speed of sound in air, related to the air thermodynamic properties. At room temperature and standard pressure, assuming the air to be an ideal gas sustaining adiabatic propagation, it is $343$ meters per second.

This equation will describe air in a rectangular room as soon as we define the unknown pressure disturbance field $p$ over all points of a 3D rectangular domain. We can describe length, width and height of the domain with the $3$ numbers $L_{x}$, $L_{y}$ and $L_{z}$. One of the domain corners can be the centre of our reference frame, and we can align the reference frame axes to run along the edges of the domain. We will call the Cartesian coordinates $x$, $y$ and $z$. Then, we make the boundaries of the domain rigid by imposing the normal component of particle velocity to be null at the boundaries. Since particle velocity is proportional to the spatial partial derivatives of the pressure disturbance we will zero those derivatives. In other words, our boundary condition is on the pressure field spatial derivatives at the room walls.

Now, we could place some kind of field source in the room, but for the time being we will only look at the normal modes of the room. In fact, due to geometry, there are special frequencies at which waves bouncing off the walls superpose just right to produce a stationary wave, that is, a wave localised in space, whose amplitude changes in time at the associated special frequency.

This will become clearer as we go along, but for the time being, let’s clarify why we search for stationary waves. The reason is that every steady state vibration of air in the room can be represented by a sum of all of these special modes, each with a different weighting coefficient. This sum is known as modal superposition and we use it to calculate the frequency response of a room. For this reason, and many additional ones, room modes are of crucial importance in room acoustics.

Property of modes is that spatial and temporal properties of the wave are factorised, which means that the pressure disturbance is the product of a function that depends only on space coordinates and a function that depend only on time:

$$p\left(x,y,z,t\right)=S\left(x,y,z\right)T\left(t\right)$$

To search for solutions of this kind we can substitute the equation above into the wave equation which, after some manipulation, brings us to the Helmholtz equation for the spatial part:

$$\nabla^{2}S=-k^{2}S$$

Which has the form of an eigenvalue problem, as we are searching for the all special $S$ and $k$ for which the action of the differential operator $\nabla$ on $S$ is that of simple multiplication by $-k^{2}$. The solution for $S$ has this shape:

$$S\left(x,y,z\right)=\cos\left(\frac{n_{x}\pi}{L_{x}}x\right)\cos\left(\frac{n_{y}\pi}{L_{y}}y\right)\cos\left(\frac{n_{z}\pi}{L_{z}}z\right)$$

Where we normalised the amplitude to $1$ both in value and unit. While, for the resonance frequencies:

$$f_{n_{x},n_{y},n_{z}}=\frac{c}{2}\sqrt{\left(\frac{n_{x}}{L_{x}}\right)^{2}+\left(\frac{n_{y}}{L_{y}}\right)^{2}+\left(\frac{n_{z}}{L_{z}}\right)^{2}}$$

$n_{x}$, $n_{y}$ and $n_{z}$ are the so called mode numbers. They are $3$ integer numbers that identify a single mode. So, one mode could be identified by $0$, $0$, $0$ (a rather not interesting one) or $1$, $0$, $0$, or $4$, $5$, $2$ and so on and so forth. For each choice of the mode numbers we get one spatial part $S$ and one resonance frequency $f$.

It isn’t maybe clear, as we did not do all the maths, but the resonance frequencies are related to the eigenvalues $-k^{2}$ (they turn out to be infinite, but discrete, too, and parameterised by our $3$ integers), and the reason why we have infinite, but discrete, modes parameterised by integer mode numbers is the boundary conditions.

Another equation will hold for $T$, which has a simple solution: an harmonic wave at the associated mode frequency:

$$T\left(t\right)=\cos\left(2\pi f_{n_{x},n_{y},n_{z}}t\right)$$

Depending on $L_{x}$, $L_{y}$ and $L_{z}$, it is possible that different modes (different choices of $n_{x}$, $n_{y}$, $n_{z}$) could have the same frequency, although different $S$. This is known as degeneracy.

The Julia Model

The equations above were implemented in few Julia functions available in this package, which can be installed in your Julia environment. To do so, launch your Julia binary and, in the REPL, press ]. This will make the REPL to go into package mode. Then, issue:

add https://gitlab.com/computational-acoustics/acousticmodels.jl.git

If you want to follow along te rest of the article you will also need something to make plots. Plots is a good choice:

add Plots

This could take some time. Once finished, press CTRL+C to exit package mode. As a note, notice that often Julia code is slow the first time it is run. This happens because precompilation.

Now that we have Plots, let’s compute a bunch of resonance frequencies. This is easier said than done, as there are few caveats. Let’s say we want to compute the first 10 resonance frequencies of the room. Let’s look at the equation for the resonance frequencies again:

$$f_{n_{x},n_{y},n_{z}}=\frac{c}{2}\sqrt{\left(\frac{n_{x}}{L_{x}}\right)^{2}+\left(\frac{n_{y}}{L_{y}}\right)^{2}+\left(\frac{n_{z}}{L_{z}}\right)^{2}}$$

Each frequency depends on $3$ integer numbers. However, it is a bit hard to figure out how to choose these numbers so that we get, say, $10$ frequencies in ascending order of magnitude. Although it is possible to came up with good efficient algorithms, let’s do something maybe a bit computationally expensive, but easy. In the Julia REPL, just include the modelling functions:

using AcousticModels

next, let’s define a 3D grid of integers:

A, B, C = index_grid(100, 100, 100)

In essence, we could think of $n_{x}$, $n_{y}$, $n_{z}$ as 3D coordinates of points in a 3D space, and the arrays A, B and C hold these coordinates. We chosen to include up to 101 values for $n_{x}$, $n_{y}$, $n_{z}$ (the indices start from 0) so to make sure we have $1030301$ of the infinite room modes (in fact, $n_{x}$, $n_{y}$, $n_{z}$ can grow indefinitely). Among those $1030301$ we should really be able to find the first lowest $10$ (or perhaps even the first few hundreds or thousands, but then we will be missing many high frequency modes afterwards, since the higher the mode numbers the more combinations are possible to produce high frequency modes).

By the way, if you press ? in the REPL, the REPL will go into help mode. Each function in the repository has its docstring. For example, searching for index_grid should return this result:

search: index_grid

  index_grid(Nx, Ny, Nz)

  Returns a 3-D grid up to indeces Nx, Ny and Nz along the 3 cartesian directions. The
  grid is represented by three matrices of size Ny by Nx by Nz.

  Example
  ≡≡≡≡≡≡≡≡≡

  julia> Gx, Gy, Gz = index_grid(3, 4, 5)
  ([1 2 3; 1 2 3; 1 2 3; 1 2 3]
  
  [1 2 3; 1 2 3; 1 2 3; 1 2 3]
  
  [1 2 3; 1 2 3; 1 2 3; 1 2 3]
  
  [1 2 3; 1 2 3; 1 2 3; 1 2 3]
  
  [1 2 3; 1 2 3; 1 2 3; 1 2 3], [1 1 1; 2 2 2; 3 3 3; 4 4 4]
  
  [1 1 1; 2 2 2; 3 3 3; 4 4 4]
  
  [1 1 1; 2 2 2; 3 3 3; 4 4 4]
  
  [1 1 1; 2 2 2; 3 3 3; 4 4 4]
  
  [1 1 1; 2 2 2; 3 3 3; 4 4 4], [1 1 1; 1 1 1; 1 1 1; 1 1 1]
  
  [2 2 2; 2 2 2; 2 2 2; 2 2 2]
  
  [3 3 3; 3 3 3; 3 3 3; 3 3 3]
  
  [4 4 4; 4 4 4; 4 4 4; 4 4 4]
  
  [5 5 5; 5 5 5; 5 5 5; 5 5 5])

Next, for each point in the 3D index grid we defined, we calculate the resonance frequency given by the coordinates of that point. For the sake of argument, let’s choose:

• The values of $L_{x}$, $L_{y}$ and $L_{z}$ respectively equal to $5$ $\text{m}$, $4$ $\text{m}$ and $3$ $\text{m}$.
• Let’s assume air in equilibrium at room temperature and ordinary pressure, making for a speed of sound of $343$ $\frac{\text{m}}{\text{s}}$.

Hence, let’s do:

F = mode_frequency.(A, B, C, 5.0, 4.0, 3.0)

And let’s get the first 10 unique modal frequencies:

sort(unique(F))[1:10]

Which in my case returns:

10-element Array{Float64,1}:
  0.0
 34.3
 42.87500000000001
 54.90679033598667
 57.16666666666666
 66.66721666439793
 68.6
 71.45833333333333
 79.26401076641136
 80.89638820738537

$0$ is of course a special case of little interest. The other modes are interesting though, and we see that modes for a room that big start right at the lower end of the audible spectrum, so we quite don’t like to have those in a listening room. Why? to better understand why we will have to consider the modal shapes too, but for the time being let’s visualise one more thing about the frequencies.

We can try to study how the number of modes changes with frequency by tiling the spectrum with many bins and counting how many modes fall within the bins. The best way to choose how big the bins should be is by having bins defined by two edge frequencies which are always at the same ratio between each other. That is, if $f_{n}^{l}$ is the lower frequency of the n-th bin and $f_{n}^{u}$ is the upper frequency of the n-th bin we want:

$$\frac{f_{n}^{u}}{f_{n}^{l}}=\text{constant}$$

This mainly for two reasons:

• Frequency normally varies over many order of magnitude, so it is quite impractical to use uniform width.
• This way matches how humans perceive pitch, as we perceive the distance between pitches as the ratio between them.

The plot below shows the numbers of modes for all frequencies lower than $4$ $\text{kHz}$ by using third octave bands which follow the property above. Third octave bands relate well to perception, but an even better choice would have been ERB bins. We focus on the modes lower than $4$ $\text{kHz}$ as the higher the frequency the higher the amount of modal frequency we are missing by computing modes on a $100$ by $100$ by $100$ grid. In short, for higher frequencies, this plot will not be accurate unless we raise the size of the grid, but that would raise computational and memory requirements.

This plot, presented on logarithmic scales, tells very interesting things. We notice that the number of modal frequencies per bin gets very high very fast. Hence, the higher frequency bins are overpopulated with modes, which means that there are many more modes per unit frequency the higher the frequency. So, we conclude that at low frequency the mode frequencies are far and sparse (up to $50$ $\text{Hz}$ we only have one mode per third octave band), while they are many and closely packed at high frequency (the third octave band close to $4$ $\text{kHz}$ has almost $150000$ modes in it). This is why room modes are a problem mainly at low frequency. Why? This is another question that the shape of the modes can answer, so let’s have a look at them!

Let’s make a couple of plots of the first modal shape, associated with the frequency of $34.3$ $\text{Hz}$. The plots are shown below.

The plots represent the mode shape at the first modal frequency, $34.3$ $\text{Hz}$, on an horizontal plane at $1.80$ $\text{m}$ from the floor. In essence, if we were a $1.80$ $\text{m}$ tall person, the plots tell us the what the value of the mode is at our head height if we were to walk to any x-y spot on the floor. Notice that the mode is not a pressure wave yet: if we put a source in the room the air, at steady state, will vibrate on a modal superposition, as we said. As we will see in the next episodes, however, at low frequency single modes dominate the superposition. Hence, we can think of the plot above as the pressure field we would have in the room if it was driven by a point source at $34.3$ $\text{Hz}$ (a part for a scaling factor controlled by the sound source).

The animation shows us the wave mode action in slow motion (representing it at $34.3$ $\text{Hz}$ would just make it hard to understand, as the animation would be really fast). We can clearly see that the wave does not travel, but that positive and negative peaks keep on exchanging periodically as time goes on. In fact, the whole spacial pressure distribution is modulated by a cosine wave at $34.3$ $\text{Hz}$, as we seen above.

We can see that there is a line right in the middle of the floor where we would not hear anything, as in there the mode is $0$, no matter how loud we would crank up our sound source. This is the problem of modes in rooms. To be able to listen to music properly, we would like our room to not alter the sound coming from the speakers. However, we see now that depending on where we are in the room there are places where, if the frequency is just right, we do not hear any sound at all!

For different modes we will have different spots of deafness. These are called nodal surfaces (as they are 2D surfaces in a 3D space). Likewise, there are gonna be spots of maximal sound pressure, which are the blue and red areas above. Close to the east and west walls of the room we would hear our $34.3$ $\text{Hz}$ at its loudest, for example.

To better understand nodal surfaces, let’s make a plot of a complex high frequency mode.

The mode plotted above is the $185.3$ $\text{Hz}$ mode, with $n_{x}$, $n_{y}$ and $n_{z}$ chosen as $2$, $3$, and $2$ respectively. It was plotted in the volume of the room as isosurfaces: each surface is drawn at a constant value of the mode. It is possible to see that, in essence, a mode is made of zones of high pressure (either positive or negative) separated by surfaces of silence. By the way, note how there are exactly $2$, $3$ and $2$ zones of high pressure per axis (you should count half zones as $0.5$). This is not a coincidence: this is what the mode numbers make to happen.

Now, think about the frequency histogram again. If we are at low frequency, the mode frequencies are far away and the mode shape are very different from one another. So, there will be certain frequencies, in the low frequency range, for which one particular spot is in a nodal surface for one mode and in a zone of high pressure for another. This means that, for a rigid room, the low frequency response is bound to be uneven.

At high frequency instead the modes are all close together. So, if one spot happens to be at a node surface for one mode, there will be a mode very close for which this is not true. What happens is that the great density of modes at high frequency smooths the response, even if the room walls are ideally rigid.

Conclusion

In this episode will only skimmed at the surface of the problem of modes in enclosed spaces, mainly focusing on the information that will be relevant to us. However, few important points were uncovered:

• Modes are stationary waves that exist in enclosed spaces.
• Any steady state vibration of air can be expressed as a modal superposition.
• Modal frequencies are sparse at low frequency and dense at high frequency.
• Modes cause the low frequency response of enclosures to be uneven.

Now we will be able to move onto modelling the same system with FEM, but that should be also done step by step:

• Getting accustomed to Elmer first.
• Getting accustomed to Elmer Helmholtz solver second.
• Solving a rigid room acoustic problem with Elmer.

These will be covered in the next episodes.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.

• Introductory
• Physics
• Modal