## Introduction

This is an article for acoustical or people who need to deal with sound in their program. Leaks through openings are a very common problem, especially in building acoustics, and therefore in realistic game situations. One should note that a lot of information on the room and source localization is based directly on our hearing. The article is therefore most useful to developers that are interested in game auralization, or if you want to incorporate some acoustical effects in your own program.

The mathematical aspects the you need or should understand to fully appreciate this article, are how to deal with:

- Complex numbers normal arithmetic operations, as well as trigonometry
- Basic calculus involving integration and derivation
- ODE and their solutions

Much of the rules of calculating these equations are quite old, meaning that they include some heavy mathematics that is not commonly treated in online articles or blogs, at least not in direct relation to this problem. It should also be said that it is also an a good starting point for learning some basic acoustic aspects, as there are lot of situations in acoustics were the material is directly transferable.

The basic problem is defined as follows. Given an incoming wave, traveling from the left to the right, it travels through the opening and comes out the other side, how much is the sound level been reduced?

## Background

The method used in the program stems from solutions of wave equation using Boundary Element Method (BEM). I will explain briefly the basics of the derivation of the equation. As a side note it is also worth knowing that acoustical simulations are relatively closely connected to problems regarding visualisation, especially in the acoustical high frequency area. The first ray-tracing program was developed by acousticians, and later adopted by visual designers. (You could read more about the programming concerning acoustical development of raytracing here).

Some explanations of what the greek letter tau actually represents in this case, before I describe the entire formula. It represents the constant, that the pressure in N/m^2 is multiplied with in order to get the resulting pressure on the other side of the wall, and its usually given as 10 *log10(tau), as we are very found of logarithms in acoustics (The reason for the obsession about logarithms are that our senses of the loudness of sound are linked to the logarithmic function). The complete general equation looks like this:

It is valid in most settings, although it's parameters The different elements in the formula are given below:

Some additional explanations of the equation are in order to properly understand what's going on here. and the hole is just a complex filter that would reduce or enhance some effects at different frequencies in the output., It is generally much easier to understand an acoustic equation of the sort above when its broken into smaller steps. There are several way to do this but the simplest one I could think of, is to construct an electric analogy circuit:

Now our equation is broken into simple steps that could be explained separately. Our source of sound is a AC voltage source, that is marked 2Pe, this would include the resistance in the surrounding material i.a. the air.

The other are explained in the order from the left to the right:

`Zr1`

, this is the radiation impedance that the opening has. It means that the incoming wave with a force F, would be translated into the movement u by a factor of the radiation impedance.`m1`

is the mass of the potential object that is placed in front of the opening. This is where the duct tape would be simulated.`F1`

is the resilience of the material`m1`

, which can be ignored for thin plate materials like duct tape. If the opening is not covered`m1`

and`F1`

will not be included in the equation.`Z3`

is the propagation from the right to the left.- The two elements of
`Z4`

are the reflected waves from one end to the other.**NB:**They are only the same values were the radius of the tube doesn't change. `m2`

and`F2`

are the elements that could cover the outgoing hole. The are calculated using similar equations in`m1`

and`F1`

.- Zr2 are the radiation impedance from the outgoing element, and in our case its the same as Zr1.

To examine the different parts I will first go through the three different kinds of acoustic impedance:

- Specific acoustic impedance. Meaning
`z = pressure/particle speed`

. - Acoustic impedance.
`Z = pressure / volume velocity`

- Radiation Impedance.
`Z`

._{r}= force / particle speed

### Radiation impedance

First thing is to define radiation impedance which is: `Z = F/u,`

where F is the force, u is the particle velocity and Z is, of course, the radiation impedance. It can be a rather difficult task to evelop this equation setp by step using integration, so I have omitted these steps, and reffer instead to the book reference given at the end of the article. A full derivation of the equation is for instance given in the book *"Fundamentals of Acoustics*" (p. 185-186).

Here the R1 term is defined as the resistance function, an X1 as the reactance function. The resistance function is defined using the solutions found in Bessels equation:

The X1 term is found using Struves function as seen below:

Radiation impedance is defined as the ratio of the force a radiator exerts on a medium, to the velocity of the radiator.

```
''' <summary>
''' Finds the radiation impedance of a circular piston
''' </summary>
''' <param name="radius">The radius of the piston</param>
''' <param name="frequency">The frequency in question</param>
''' <returns>A complex radiation impedance</returns>
''' <remarks></remarks>
Public Function RigidPistonRadiator(ByVal radius As Double, ByVal frequency As Double) As Complex
'If the user, for some silly reason puts in 0 or a negative number as the radius
If radius <= 0 Then
Return (New Complex(0, 0))
End If
Dim r1, x1 As Double
'Medium constants
Dim c0 As Double = 340
Dim p0 As Double = 1.21
'Wavenumber
Dim k As Double
k = 2 * System.Math.PI * frequency / c0
'The resistance function
r1 = 1 - (2 * Acoustics.SpecialFunctions.BesselJ(2 * k * radius, 1) / (2 * k * radius))
'The reactance function
x1 = 2 * Acoustics.SpecialFunctions.Struve(2 * k * radius, 1) / (2 * k * radius)
'Stor the results from calcualtion
Dim Result As New Complex(r1, x1)
'Multiply with the caracteristic impedance with area
Result = System.Math.PI * radius ^ 2 * c0 * p0 * Result
Return Result
End Function
```

### Mass impedance

I actually made a simplification in modeling the damping effect of simple duct tape or similar materials. For a complete mathematical treatment, you need to include the stiffness of the material also, but thats usually an irrelevant property in thin lightweight materials intended for use in this program. You can actually derive the mass law in acoustics from these equations. It simply states that a thin wall surrounded by air (Z_{0} in the equation below), the damping at a given frequency, is dependent on the frequency times the mass (and by mass I mean the total mass per square meter).

The mass law was initially developed by measuring the sound insulation of normal housing walls. It is also considered to be the maximum possible reduction level of any frequency. Normally it is not valid in the higher frequencies, as they are determined by resonant capacities of studs and the limit frequency of the covering plate as well as the distance between the plates.

### Characteristic impedance and propagation function in a medium

The reason that characteristic impedance is important so you can see that it is useful, as it links the particle velocity with pressure by a constant.

It is however a special case of the solution to the one dimensional wave equation, and not valid for all conditions. It basically describes the resistance in a pressure wave multiplied with the particle velocity.

The real solution for spherical wave is actually `Zk = rho C *(jkr/(1+jkr))`

and only when` kr >> 1`

will the characteristic impedance be` rho c.`

If one assumes the airflow resistivity is dependent on only on a constant times the particle velocity one gets an one dimensional model (called Rayleigh-model) of loss in a medium. Its the simplest model for losses in a medium we have, and some modified version of it forms the Delany-Bazley porous loss model. This loss is in addition of normal geometry losses in a medium, but adds a consistent damping per meter. This is not included in the article, but should definitely be on your todo list.

## Special mathematical functions used in the program

I have three special functions that are needed to use this program, and I have placed these inside a module. The source of the code is also given. They do originally come from the book "Numerical methods for C#" by Jack Xu (with thanks to the author for allowing me to publish this part). However, the Gamma function is changed based on the references that are given, and the Struve function is my own implementation based on the Wikipedia link given in the code comment.

```
Namespace Acoustics
Module SpecialFunctions
''' <summary>
''' http://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind_:_J.CE.B1
''' </summary>
''' <param name="x">Jv(X) were x is the number and v is the order</param>
''' <param name="v">Bessel function order</param>
''' <returns></returns>
''' <remarks></remarks>
Public Function BesselJ(ByVal x As Double, ByVal v As Double) As Double
Dim sum As Double = 0.0R
Dim term As Double = 0.0R
Dim i As Integer = 0
Do
term = System.Math.Pow(-1, i) * System.Math.Pow(0.5 * x, 2 * i + v) / Gamma(i + 1) / Gamma(i + v + 1)
sum += term
i += 1
Loop While System.Math.Abs(term) > 0.000000000001R
Return sum
End Function
''' <summary>
''' http://en.wikipedia.org/wiki/Lanczos_approximation
''' </summary>
''' <param name="x">number x dependent on gamma function</param>
''' <returns>Gamma of the number x</returns>
''' <remarks></remarks>
Public Function Gamma(ByVal x As Double) As Double
Dim g As Integer = 7
Dim p As Double() = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 0.0000099843695780195716, 0.00000015056327351493116}
Dim y As Double
If x < 0.5 Then
Return Math.PI / (Math.Sin(Math.PI * x) * Gamma(1 - x))
End If
x -= 1
y = p(0)
For i As Integer = 1 To g + 1
y += p(i) / (x + i)
Next
Dim z As Double = x + (g + 0.5)
Return System.Math.Sqrt(2 * System.Math.PI) * System.Math.Pow(z, x + 0.5) * System.Math.Exp(-z) * y
End Function
''' <summary>
''' http://en.wikipedia.org/wiki/Struve_function#Power_series_expansion
''' </summary>
''' <param name="x">The Hv(x) were x is the number and v is the order</param>
''' <param name="v">Struve function order</param>
''' <returns></returns>
''' <remarks></remarks>
Public Function Struve(ByVal x As Double, ByVal v As Integer) As Double
Dim result, term, sum As Double
'Check if the first term is 0, if so then everything
' would be zero
If (0.5 * x) ^ (v + 1) = 0 Then
Return result
End If
Dim i As Integer = 0
Do
term = ((-1) ^ i) * (0.5 * x) ^ (2 * i)
term = term / (Gamma(i + (3 / 2)) * Gamma(i + v + (3 / 2)))
'Avoid adding NonNumbers as they will couse termination
If Not Double.IsNaN(term) Or Not Double.IsNegativeInfinity(term) Then
sum += term
End If
i += 1
Loop While System.Math.Abs(term) > 0.000000000001R ' i < 100
If Double.IsNegativeInfinity(result * sum) Or Double.IsNaN(result * sum) Then
Return 0
Else
Return result * sum
End If
End Function
End Module
End Namespace
```

The three simple functions are used quite frequently in several different acoustic problems, and will often form the cornerstone of many numerical solutions. The Bessel function is indeed the solution to the wave equation, so its hardly surprising. The Gamma function is also used indirectly in many methods, also outside the realm of acoustics.

## History

You should know that the program is rather simple in nature, and does not include the behavior of the reduction when you put porous material inside the tube. The results are also normalized to 1 m^2, meaning it is created for comparison from different hole sizes. You should remove the term 1/area in the program to get it's real reduction value.

The program could be used in the following way. assume that you have a noise source on the outside of the wall (like a road) that produces 40 dB in each frequency band outside your wall (not far from a heavily traffic road), This would mean that the noise inside your room would be perceived as 40 - the reduction at any given frequency (it's a little more complicated than that, as it is also determined by the reduction of the wall, and of the reverberation time inside the room).

Jack Xu should also be thanked for allowing me to include his modified source code from the book "Numerical methods for C#" in this program.

## References

*"Building acoustics"*, First edition (2008), Tor Erik Vigran, CRC Press*"Formulas of Acoustics", 2nd Edition, F.P. Mechel, Springer**"Handbook of mathematical functions"*, Abramowitz and Stengun (1970), Dover*"Practical Numerical Methods with C#"*, Jack Xu,*"Fundementals of Acoustics"*, Fourth edition, Kinsler, Frey, Coppens and Sanders, John Wiley & Sons, Inc.

Online links: