Sound in Games - Acoustical Leaks in Walls

• Download source code and demo - 166.7 KB

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 in the real world is based directly on our hearing, as we often estimate the position of others by hearing where the sound came from. It is however quite complicated to simulate the real situation, as we use the reflection from our body to determine the location of height where the sound came from, and the difference in time from one ear to the other to determine the x-y position of sound. This approach is called Head Related Transfer Function (HRTF) in Acoustics for people who want to investigate this further.

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 that you need or should understand to fully appreciate this article, are how to deal with:

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

I will however show you some simple equivalent ways of thinking about the problem, if the math is too heavy. You should also be aware that the equations could be constructed using a similar equivalent in either the electrical circuit using resistors, capacitors and inductors or a mechanical configuration using springs viscous damper and masses.

The equation I'm using to develop the formula is quite old fashioned as it deals with Bessel function, which is a general solution to the spherical wave equation among other things, and the article will therefore 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 a good starting point for learning some basic acoustic aspects, as there are lot of situations in acoustics where some of the material is directly transferable.

The basic problem is defined as follows. Given an incoming sound wave from outside (traveling from the left to the right) with an angle of incident based on the normal of the hole, it travels through the hole and comes out the other side. The hole has the diameter of 2a on both sides and the length or depth of the tube is defined as d. How much has the sound level been reduced?

I will expand the original article by including the effect of putting mineral wool inside the hole. This means I can simulate the effect of:

• An empty circular hole
• Putting a light weight material in front of one or two of the holes (like duct tape)
• Putting mineral wool inside the hole that covers the entire length of the duct
• And putting mineral wool inside and duct tape on each side of the openings

With this program, you should be able to understand the effect that the different configurations will have on the end result, or the sound pressure inside. You should also take note of the fact I'm only calculating the normal incident sound wave, in the calculation, but it is easily expanded if you change the angle theta of the incoming wave.

Background

Before I'm going to show you the equations, I'll quickly give you some simple analogies that are used in constructing the solution, or rather how you could think of this problem.

• A room, or a similar closed space could be modeled as a Capacitor or as a Spring
• A pipe, or a leak in a room, could be modeled as an Inductor or as a Mass
• A matrix of thin pipes could be modeled as a Resistor of as a Viscous Damper. For further examination of the acoustical idea, you could see this article.

The solution to the transmitted wave is actually a filer, where the different materials change some of the resonances and add new poles or zeros to the transfer function.

The method used in the program originally stems from solutions of wave equation using Boundary Element Method (BEM), and 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 computer 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 it's 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 perceived loudness of sound are linked to the logarithmic function). The complete general equation looks like this:

It is valid in most settings, although its parameters are idealised, you should always remember that if it doesn't work in the idealised calculation, it won't work in the real world.

A brief explanation for 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 it's 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.e., the air.

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

1. 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.
2. 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.
3. 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.
4. Z3 is the propagation from the right to the left.
5. The two elements of Z4 are the reflected waves from one end to the other. NB: They are only the same values where the radius of the tube doesn't change, and that everything else is also symmetrical.
6. m2 and F2 are the elements that could cover the outgoing hole. The are calculated using similar equations in m1 and F1.
7. Zr2 are the radiation impedance from the outgoing element, and in our case, it's the same as Zr1.

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

1. Specific acoustic impedance. Meaning z = pressure/particle speed.
2. Acoustic impedance. Z = pressure / volume velocity
3. Radiation Impedance. Z<sub>r</sub> = force / particle speed.

They are used in different circumstances, dependent on the equation/problem you are trying to solve. However, one important aspect of combining equations in acoustics is by assuming continuity of particle velocity, meaning the particle velocity can't abruptly change from one place to another. This is vital to understand in many situations where you are trying to combine different types of materials together.

Radiation Impedance

The 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 refer 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 calculation
            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 that's 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<sub>0</sub> _{in the equation below that is the characteristic impedance of air), 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).}

We proceed by finding the absorption coefficient for the material:

We also know from our initial equation for Zg that there is no internal or any other losses in our assumptions, therefore the only way it can get absorbed must be by exiting the other side of the plate. We then write up the solution with replaced values:

The sound reduction could also be calculated:

We insert the values for the characteristic impedance of air and clean up the expression, and arrive at this equation:

The mass law is used quite frequently in predicting the sound insulation at low frequencies of normal housing walls, although the coefficient 42.5 is changed to 47 as a real wall is more complex. It is also considered to be the maximum possible reduction level in the low frequency area, where the sound reduction is normally decided by the share weight of the total wall. It should be noted that for a simple homogenous wall it usually fits rather well.

Normally, it is not valid in the higher frequencies, as a normal wall is also determined by resonant capacities of studs and the limit frequency of the covering plate as well as the distance between the plates. The sound reduction could also be much higher than the mass law predicts above the mass controlled area.

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. It's the simplest model for losses in a medium we have, and Mechel uses this simple model in the low frequency area and some empirical measurements to make it fit at high frequencies. The empirical models have to be used simply because of the complexity of the material configuration in mineral and glass fibre. This loss is in addition of normal geometry losses in a medium, but adds a consistent damping and phase shift per meter. This details of this is not included in the article and I would refer you to this article instead.

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) where 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) where 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 cause 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 in other fields as well. The Bessel function is indeed the solution to the wave equation, so its hardly surprising.

The Gamma function is one of the most universal function in mathematics, and is the backbone in all the functions I use, and it can also be used in statistics among other.

History

This is an updated version that includes the damping effect of a porous material inside the tube, as I initially had just included the possibility of simulation lightweight cover materials to be placed in front or at the back of the tube.

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 its 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.

It should also be noted that it generally underestimates the damping effect at the higher frequencies for an empty hole, as the current implementation does not take into effect the viscous losses or any other effects, which could be significant for an empty hole. With plate covering or with mineral fibre, the calculation fits very well with the expected damping.

References

• "Building acoustics", First edition (2008), Tor Erik Vigran, CRC Press
• "Formulas of Acoustics", 2^nd Edition, F.P. Mechel, Springer
• "Room Acoustics", Fifth edition, Heinrich Kuttruff, Spon press
• "Handbook of mathematical functions", Abramowitz and Stegun (1970), Dover
• "Practical Numerical Methods with C#", Jack Xu,
• "Fundamentals of Acoustics", Fourth edition, Kinsler, Frey, Coppens and Sanders, John Wiley & Sons, Inc.

Online links:

• http://www.phys.unsw.edu.au/jw/z.html