Introduction
Simulations of acoustic waves are a topic of intense interest in many circumstances, but they are difficult to achieve with minimal use of computer power. There are several candidates,
but they all have pros and cons, as we shall see. One of the most precise simulation techniques is the Transmission Line Matrix, TLM for short, as it is capable of simulating any situation
given the right preconditions. Some of the best known available techniques in acoustics are:
Ray tracing is normally much used in concert hall acoustics, as the dominant domain (the frequencies it is valid for) often falls within the range where one can apply ray tracing to.
In general this is valid for the higher frequency domain, were a sufficiently large portion can be
The so called Finite Difference Time Domain technique (FDTD for short) is actually the exact
same for acoustics, as the FD-TD method solves the wave equation while the TLM provides a physical model of propagation, this means that the limitations and advantages are exactly
the same for the two. It should be mentioned that the TLM and FDTD can be used both in acoustics and in electrodynamics, as they are both waves that can be
modeled by voltage and
currents in electromechanic waves, and pressure and particle velocity in acoustics.
The FEM method is quite popular algorithm used in many commercial products, and is especially
useful for simulating wave equation at the lower frequency domain, as the method must have points in accordance to the frequency that its simulates. The FEM technique solves
the Helmholtz wave equation, were one needs to insert the frequency one wish to solve for first.
The BEM method is quite popular to generate exact solutions for more limited problems, for instance like the propagation of sound waves through slits and other openings. It very quickly gets complicated if the geometry gets more varied, used in simple cases and is an analytical technique at hart.
The goal of this article is to show just how simple a wave simulation could be done by, and how intuitive the process of designing it actually is. There are also more recent approaches that involves the Lattice Boltzmann Model, and this could be used for the linear wave equation as well as the Navier-Stokes equation.
If you want to know more about the details behind each of these techniques I would recommend to read this PDF document written by one of the Professors I had in numerical acoustics.
Background
The principle behind the TLM method was first described by Huygens, a
Dutch scientist responsible for many developments in both mathematics and physics (he also build the first orrery based on continued fraction), were he reduced the wave propagation to a series of elastic collisions of balls, which incidentally is actually nearly exactly physically right, as an acoustic wave travels with the help of collisions of molecules, that is surrounded by an electromechanical field (the electrons circling around each of the nucleus), and the space between the molecules.
The idea to use Huygens principle to compute the wave equation was first developed by Peter Johns, who also coined the name Transmission Line Matrix (TLM for short), as the analog circuit resembled a series of interconnected points on a matrix grid. The idea is quite simple take an incoming wave at one of the nodes (illustration from Wikipedia):
The first picture (left most), shows an incoming wave at the center node in the diagram, travelling north from the point south of the center point.
The pressure is then scattered (the middle picture), but the total energy remains the same (the sum of all the scatter is still 1), and the wave south is a reflection due to the fact that it travels in one line, and at the node, it sees 3 lines that it can propagate further into, and this leads to the reflection and scatter coefficients calculated by the characteristic impedance the wave sees as it enters the node. In mathematical terms it can be written as follow (from Wikipedia again), for transmission:
And the corresponding reflection coefficient:
The last picture just shows how the pattern repeats over and over. The tree pictures shows the transmission lines of an electrical circuit, while the TLM model in acoustics uses a series of tubes instead, as shown below, where the nodes are shown in gray:
This information is actually all you need to simulate a simple propagation in free space on a computer, but we are interested in the more complicated scenarios, and then the mathematical theory and implications is important to know.
The acoustic tubes needs some additional explanation, as they must have some dimensions too in order to function properly. This is not set explicitly, as one in the TLM model assumes that only plane waves will travel inside the tube from one node to the next, this could actually be described by a matrix of equations:
The angle is usually not included in this model as one assumes that only normal incident plane waves travel inside the tubes. One can also note that the length (or distance) is included in the above matrix, but not in the model. For the model to be valid one needs to choose a small enough distance so that the pressure doesn't change over the tube. The above equation is often used to calculate the behavior of acoustic absorbers, and propagation through walls etc. It also implies that there is a cutoff frequency, were our model is no longer reliable.
In the TLM model it is more interesting to see the propagation matrix in a different form, as one goes from one node to the other. The equation is then dependent on 4 lines of possible travel, so the expression for a wave that propagates from west to east becomes:
The angle described in the above equation is the length l between the nodes and the simulated wavelength and the relation is given below:
This principle could be used to fix a problem for calculating the propagation through thin walls, as they are notorious for creating a near field, and by this one means that some of the sound wave does not propagate from the wall to the room, but just travels perpendicular to the wall, and thereby cancelling the effect. This is most typical at higher frequencies for normal lightweight wall constructions found in houses. I realize that these last three equations are for specially interested people, but it is very useful to know, if you want to use this in the real world. This just goes to show that it doesn't matter how well a computer program you have, as the results depends on the level of knowledge that the user possesses.
The basic TLM propagation model in code
The basic principles have now been established, but this is the CodeProject after all, so we need to take a look at what really needs to be stored for calculating a progression of waves. By a closer review of the model we can establish that the following events might happen in this order for each of the discrete time steps:
- If needed, create the sound pressure for each sound source by setting the incoming pressure
- Calculate scatter from all the incoming sound pressures for each node
- Propagate sound pressure from the incoming to the opposite outgoing node
- Summation of the incoming pressure components gives the current pressure in each node
This clarification is all that is really needed to draw a scheme of the variables needed in each node, were this is connected in each of the directions to a similar node:
Now it is all a matter of counting the number of inputs and outputs for the node, which are two in each of the
four directions, as well as the current pressure. The addverbreations are given the names for the the placement of the connections and if it travels inward or outward. In sum, for a simple propagation model without losses we need 9 matrixes that holds all the variables; like this:
Private Dimensions As Integer = 100
Private _
IncommingEast(Dimensions, Dimensions), _
InncommingNorth(Dimensions, Dimensions), _
IncommingWest(Dimensions, Dimensions), _
IncommingSouth(Dimensions, Dimensions), _
ScatteredEast(Dimensions, Dimensions), _
ScatteredNorth(Dimensions, Dimensions), _
ScatteredWest(Dimensions, Dimensions), _
ScatteredSouth(Dimensions, Dimensions) _
As Decimal
Private CurrentPressure(Dimensions, Dimensions) As Decimal
The algorithm is iterative, as one needs to loop through all the points in the matrix, and I start off by setting the pressure at a given node, and afterwards calculate
the scattered pressures. It is important to calculate all the scattered pulses in a separate for loop before you add up the current pressure.
This actually makes logical sense, as the scattering for all the elements would have to happen before you can say what the next pressure is going to be.
For x As Integer = 0 To Dimensions - 1
For y As Integer = 0 To Dimensions - 1
If Source.Contains(New Point(x, y)) Then
IncommingEast(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + IncommingEast(x, y)
IncommingWest(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + IncommingWest(x, y)
IncommingSouth(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + IncommingSouth(x, y)
InncommingNorth(x, y) = SinAmplitude * Math.Sin(2 * Math.PI * TT / 10) + InncommingNorth(x, y)
End If
ScatteredNorth(x, y) = 0.5 * (+IncommingEast(x, y) - InncommingNorth(x, y) + _
IncommingWest(x, y) + IncommingSouth(x, y))
ScatteredEast(x, y) = 0.5 * (-IncommingEast(x, y) + InncommingNorth(x, y) + _
IncommingWest(x, y) + IncommingSouth(x, y))
ScatteredWest(x, y) = 0.5 * (+IncommingEast(x, y) + InncommingNorth(x, y) - _
IncommingWest(x, y) + IncommingSouth(x, y))
ScatteredSouth(x, y) = 0.5 * (+IncommingEast(x, y) + InncommingNorth(x, y) + _
IncommingWest(x, y) - IncommingSouth(x, y))
Next
Next
The scattered pressure is outgoing, and I need to take into account that each node can have incoming pressures coming from all of the different 4 connection points (North,
West, East and South):
For x As Integer = 0 To Dimensions - 1
For y As Integer = 0 To Dimensions - 1
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
CurrentPressure(x, y) = (0.5 * (ScatteredEast(x, y) + _
ScatteredNorth(x, y) + ScatteredWest(x, y) + ScatteredSouth(x, y)))
Next
Next
Both of these loops happens in each discrete time step, but this code will not work properly without a Perfectly Matched Layer (PML), as we have yet to add end corrections.
The code above will only work if the actual point is surrounded on all sides by air (or the same object in all points). It is now time to introduce boundaries and objects with
reflection factor.
Reflections and terminations of the grid
The reflections have been used and described previously, but in that case it was rather simple as all the acoustic tubes were the same. Reflections and transmissions
are in reality just information that says how much resistance does the wave see when it reaches the node. If the resistance is high it will be mostly reflected,
and if it is low it would mostly transmit.
In acoustics the reflection is determined by the characteristic impedance of the medium, and the difference between
the two coefficients. The reflection factor is defined by the relation:
Where Za is the characteristic impedance of the wall or boundary , while Zt is the characteristic impedance of the free field in the TLM model. The Zt value
is actually not the normal value used in air, and the reason for this is actually the discretisation of the mesh, as it is squared. This has profound implications,
as the speed of propagation (defined as C0 in air with the value of 343 m/s) in 2D is actually c0 / Sqrt(2) (and c0 / Sqrt(3) in 3D). This is just the
Pythagoras of the diagonal travel.
To use this information in the TLM structure we need to calculate the reflection by this method:
This has some interesting effects, for instance; if you want your mesh to have the boundary condition of an infinitely large TLM matrix, you need to match
Za with Zt (equal values, normally called PML) and this you get the following non reflective boundary condition:
Intuitively this actually makes sense, as one would expect that it is not zero, as this would lead to infinitely
absorption, but that the particles actually collide, and some of the energy is transmitted back to the mesh.
The code that makes the reflections from the PML viable is thus given below:
For x As Integer = 0 To Dimensions - 1
For y As Integer = 0 To Dimensions - 1
If x = 0 And y = 0 Then
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncommingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = 0 And y = Dimensions - 1 Then
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncommingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf x = Dimensions - 1 And y = 0 Then
IncommingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = Dimensions - 1 And y = Dimensions - 1 Then
IncommingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncommingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf x = 0 Then
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf y = 0 Then
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf y = Dimensions - 1 Then
IncommingEast(x, y) = ScatteredWest(x + 1, y)
InncommingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncommingWest(x, y) = ScatteredWest(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
ElseIf x = Dimensions - 1 Then
IncommingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
IncommingWest(x, y) = ScatteredEast(x - 1, y)
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
If Wall.Contains(New Point(x + 1, y)) Then
IncommingEast(x, y) = 0.5 * ScatteredEast(x, y)
Else
IncommingEast(x, y) = ScatteredWest(x + 1, y)
End If
If Wall.Contains(New Point(x, y + 1)) Then
InncommingNorth(x, y) = 0.5 * ScatteredNorth(x, y)
Else
InncommingNorth(x, y) = ScatteredSouth(x, y + 1)
End If
If Wall.Contains(New Point(x - 1, y)) Then
IncommingWest(x, y) = 0.5 * ScatteredWest(x, y)
Else
IncommingWest(x, y) = ScatteredEast(x - 1, y)
End If
If Wall.Contains(New Point(x, y - 1)) Then
IncommingSouth(x, y) = 0.5 * ScatteredSouth(x, y)
Else
IncommingSouth(x, y) = ScatteredNorth(x, y - 1)
End If
End If
CurrentPressure(x, y) = (0.5 * (ScatteredEast(x, y) + _
ScatteredNorth(x, y) + ScatteredWest(x, y) + ScatteredSouth(x, y)))
Next
Next
The TLM is a very intuitive model, and one can easily understand what's actually going on without too much effort (at least in the simple cases). If one takes the example here:
InncommingNorth(x, y) = 0.5 * ScatteredNorth(x, y)
I have just set all the walls (not the non-reflective boundary) to the value 0.5, and I just take the energy that was scattered north, and send it back in a loop to the incoming north, only that it is damped by the reflection factor. In the picture of the TLM node it corresponds to short circuit the model in the North end by a resistor. The side walls are set to a non reflective boundary condition (the actual definition of a PML) where the actual value of R is around -0.17.
The code snippets shown in the article is somewhat altered in the actual program, as there are faster ways of checking if it is a free field node. The idea is to set up an array of boolean that is equal in size to the nodes. Here we specify if any node have a wall, or has a neighbor wall in the North, South, East or West. This is done when you add a wall element:
Public Sub SetWalls(ByVal WallPoints As List(Of Point), ByVal ReflectionFactor As Decimal)
Dim R As Decimal = ((1 + ReflectionFactor) - SquareRoot(2) * (1 - ReflectionFactor)) / _
((1 + ReflectionFactor) + SquareRoot(2) * (1 - ReflectionFactor))
For Each p As Point In WallPoints
Walls(p.X, p.Y) = 0.5
WallNeigbor(p.X, p.Y) = True
If p.X <> 0 Then
WallNeigbor(p.X - 1, p.Y) = True
End If
If p.X <> Dimensions - 1 Then
WallNeigbor(p.X + 1, p.Y) = True
End If
If p.Y <> 0 Then
WallNeigbor(p.X, p.Y - 1) = True
End If
If p.Y <> 0 Then
WallNeigbor(p.X, p.Y + 1) = True
End If
Next
End Sub
The code also stores the specified Reflection Factor for each wall elements as it gets added. This is also stored in an array, and can be called everytime a wall element is encountered.
One can now determine if it is a freely propagating field by simply checking if the Wall Neighbor is true:
If Not (WallNeigbor(x, y)) Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
IncomingWest(x, y) = ScatteredEast(x - 1, y)
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
If WallNeigbor(x + 1, y) Then
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
Else
IncomingEast(x, y) = ScatteredWest(x + 1, y)
End If
If WallNeigbor(x, y + 1) Then
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
Else
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
End If
If WallNeigbor(x - 1, y) Then
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
Else
IncomingWest(x, y) = ScatteredEast(x - 1, y)
End If
If WallNeigbor(x, y - 1) Then
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
Else
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
End If
End If
Else
If x = 0 And y = 0 Then
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
IncomingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
IncomingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = 0 And y = Dimensions - 1 Then
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
InncomingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
IncomingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
ElseIf x = Dimensions - 1 And y = 0 Then
IncomingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
IncomingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
ElseIf x = Dimensions - 1 And y = Dimensions - 1 Then
IncomingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
InncomingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
If Walls(x - 1, y) = 0 Then
IncomingWest(x - 1, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x - 1, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
ElseIf x = 0 Then
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
IncomingWest(x, y) = OpenEndedReflection * ScatteredWest(x, y)
ElseIf y = 0 Then
IncomingSouth(x, y) = OpenEndedReflection * ScatteredSouth(x, y)
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
ElseIf y = Dimensions - 1 Then
InncomingNorth(x, y) = OpenEndedReflection * ScatteredNorth(x, y)
If Walls(x + 1, y) = 0 Then
IncomingEast(x, y) = ScatteredWest(x + 1, y)
Else
IncomingEast(x, y) = Walls(x + 1, y) * ScatteredEast(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
ElseIf x = Dimensions - 1 Then
IncomingEast(x, y) = OpenEndedReflection * ScatteredEast(x, y)
If Walls(x, y + 1) = 0 Then
InncomingNorth(x, y) = ScatteredSouth(x, y + 1)
Else
InncomingNorth(x, y) = Walls(x, y + 1) * ScatteredNorth(x, y)
End If
If Walls(x - 1, y) = 0 Then
IncomingWest(x, y) = ScatteredEast(x - 1, y)
Else
IncomingWest(x, y) = Walls(x - 1, y) * ScatteredWest(x, y)
End If
If Walls(x, y - 1) = 0 Then
IncomingSouth(x, y) = ScatteredNorth(x, y - 1)
Else
IncomingSouth(x, y) = Walls(x, y - 1) * ScatteredSouth(x, y)
End If
End If
End If
As it is assumed that most nodes is in the free field this is quite a speed advantage, although the overall code got a lot more complicated by this.
There is also something special about the end nodes, as a wall can also appear at the end of the model, this implies that the node at the top, with a non reflection in the north, can have reflections form a possible wall element in East or West (or both for that matter).
I should mention that this specific way of modeling each of the reflective surface does not allow them to have different absorption and reflection at different frequencies. If one want to use a pulse to find the transfer function between two points, one would have to do this one frequency at the time, and combine the results. There is another possibility, namely to use a FIR filter as the absorption and transmission coefficients, as normal walls and other boundaries is usually a lowpass filter in real life, but this requires some knowledge of DSP (Digital Signal Processing) and one also needs to know what the absorption coefficients are of each material. However this goes far beyond explaining how simple TLM models can be!
A different type of wall
The observant reader have perhaps noticed that the walls used so far can not transmit a wave through the wall, and the energy that is not reflected simply disappears, and this cannot happen in the real world, of course. There is however a problem in implementing such walls in the TLM model if you want to use a uniform mesh, but fortunately we can get around this problem by two methods of tertiary to solve it.
We cannot increase the speed of sound in the TLM model from its original value, but we can however decrease it, and thereby creating a reflection and transmission of lets say two gases. The normal way of implementing this is to use two special branches in each node, that deals with both changes in the speed of sound and attenuation of the sound wave.
The possible changing of the speed of sound also makes simulation of a Non-linear propagation possible, where the speed is changed according to the current pressure at the node. This is not included in the view of course as it is a highly specialized task, but the code in the TLMGrid is fully functional if it is set properly.
We can also set the reflection and transmission dependent on the weight of the two different materials, this would generate a reflection and transmission if the propagating waves sees two different weights while travelling from one node to the next, which also would result in a real reflection and transmission situation in the mesh. However don't try and combine this type of wall with the all absorbing wall, as it would confuse the model completely, and the results would go nuts. A quick introduction from some experts are given in this
PowerPoint presentation (Pdf document).
At last I should mention that the TLM model currently only works with longitudinal waves, and the waves that is produced by a thin wall (with moment of inertia) cannot be combined in this simulation. You would then have to resolve this with either some fancy FIR filter approach or by solving the equations in a FDTD model.
Relating the frequency in the model to the real world
There are many issues in relating the TLM model to the real world, as you have already seen with the speed of sound in the model, which is C0/Sqrt(2) in each direction. One has thus demanded that the real speed C0 appears as the hypotenuse of an equal sided triangle, with the propagation angle of 45 degrees, as seen in the illustration below:
The reason that the 45 degree angle is that this is the worst possible angle for the model to handle correctly, as it has to travel in the x direction and in the y direction to get to the top point. If you set the model to draw the absolute pressure for anything other than 0, to be colored in red, you will see that the model will propagate quite fast in the straight x and y direction, but the amplitudes of the points outside the length of the circular 45 degree radius are extremely small. The wave would then appear as a square around the source point.
There are two similar way of looking at this and to reach a conclusion about how these are related. One way of looking at this is by reasoning that it would take two discrete time pulses for the pressure to reach the point at 45 degrees, as it would have to travel first along the x axis and then along the y axis (or vice a versa) meaning that it would travel 2 Delta L at the speed of Ct.
The standard formula that relates the wavelength or the frequency of a standing wave in a propagation medium (which corresponds to the wave length of the medium) is;
However the sine source in the TLM model needs only to have the time step t as its input (but it could also be adjusted by a constant phase displacement, so you could delay the onset of the Sine wave by a fixed amount. This would allow you to use two sources with different phase to simulate a dipole source ). The generating sine function with a set frequency is thereby given as follows:
To avoid confusion I must emphasize that it is the t that is halfed not the l, as it is written in the equation above it might be confusing. The result is that the frequency is effectively halved
when it is put into the model. Finally, the code for the actual Sine source in the TLM model:
If Sources(x, y) Then
Dim Sine As Decimal = Amplitude * Math.Sin(Math.PI * CurrentTimeStep * _
DeltaLength * SourceFrequency / CT)
IncomingEast(x, y) = Sine + IncomingEast(x, y)
IncomingWest(x, y) = Sine + IncomingWest(x, y)
IncomingSouth(x, y) = Sine + IncomingSouth(x, y)
InncomingNorth(x, y) = Sine + InncomingNorth(x, y)
End If
There is also other ways of generating sources besides the sine function. The most used, and perhaps the most valuable is the Dirac impulse. In theory this is infinitely high point source at an infinitely small time gap, which in mathematical theory makes this a completely flat frequency response at the source, or in other words all frequencies has the same amplitude if a Fourier transform is performed on the source signal. This is extremely easy to code, as it would just have an amplitude when the time step is equal to t=0.
If DeltaSources(x, y) Then
IncomingEast(x, y) = Amplitude + IncomingEast(x, y)
IncomingWest(x, y) = Amplitude + IncomingWest(x, y)
IncomingSouth(x, y) = Amplitude + IncomingSouth(x, y)
InncomingNorth(x, y) = Amplitude + InncomingNorth(x, y)
End If
Finally I mentioned something about the cutoff frequency, where above this limit the TLM is no longer valid. The cutoff is defined as
dl/Lambda = 0.25, but the recommended maximum is dl/Lambda = 0.1, if the frequency you are interested in simulated is above this, you cannot
presume that the results will be accurate. There is however the possibility that you could get a stable output even if you exceed this limit
due to anti aliasing effects.
This results into 10 or more discrete samplings of a Sine wave over a complete period, and the idea is that each of the steps in the Sine could be
modeled good enough so that it will be approximately a straight line. This is by the way not unique to just the TLM model, as both the FDTD
and TLM is dependent on just using the first derivative of the Taylor series,
and if we need to model non linear effects we also have to make adjustments so that the difference between every discrete sample can be
accurately modeled by a straight line.
The cutoff frequency is actually a reason to not use the Dirac delta pulse, as it is guaranteed to include frequencies above the cutoff frequency, and thereby ruin the nice
transfer function effect at the measurement point, and the propagation will also be "corny" in the model, as it would seem to include some random noise in the figure.
For this reason it is more convenient to use a so called Gaussian Pulse as a source instead of the Dirac delta impulse,
and this is defined as:
And this expression will generate the dirac delta pulse when a goes to 0. Remembering that we don't want any frequency components above the TLM cutoff frequency we also need to know
the cutoff frequency (where the amplitude in the frequency spectrum is halved, and this could readily be found using the theory
of Gaussian Filters, and can be found here).
Gaussian sources are however a bit troublesome, as the frequency response degenerates is quite slow. In practice this would imply that the pulse would get extremely long
for very high frequencies, so that a different type of source might be better. A source that uses Butterworth or Chebyshev filter characteristics might be a better way of doing this,
but in this article I will ignore this.
Constructiong the 3D viewer with speed improvements
The design of the 3D viewer is heavily influenced by Sean Sexton's raindrop animation (you can download
his original C# code here). He has however made some simplifications, as he doesn't have to worry about walls and other objects that his wave
is going to travel through. His first implementation did use a delta dirac impulse as a source directly, while his second version used a smoothing filter that is approximately
a simple Gaussian filter. (The Gaussian filter is basically a low pass filter that removes high frequency components.) By using the Gaussian source directly, one get's a smoother surface:
There really isnt much left of his code in my implementation, except for the tips given in this MSDN article,
which basically lays out how to initialize and update a 3D model in WPF the fastest way possible. The basic idea is to initialize all the 3D arrays once, and connect one list of 3D points
to the viewmodel, and when you update the calculations you should update the array that isn't currently linked to the viewmodel, and then swap the two arrays.
Dim RunTimeInMS As Double = 1
Dim LastTimeRendered As Double = 0
Dim IsSimulationRunning As Boolean = False
Dim RunSimulationEvent As New EventHandler(AddressOf RunSimulation)
Private Sub btnStart_Click(sender As System.Object, _
e As System.Windows.RoutedEventArgs) Handles btnStart.Click
If IsSimulationRunning Then
RemoveHandler CompositionTarget.Rendering, RunSimulationEvent
btnStart.Content = "Resume"
IsSimulationRunning = False
Else
AddHandler CompositionTarget.Rendering, RunSimulationEvent
btnStart.Content = "Pause"
IsSimulationRunning = True
End If
End Sub
There is something that is a bit puzzling that is namely that you need the EventHandler in VB, but not in C#. If you just code it like shown below, you would
get a warning that won't go away. I presume that it has something to do with that VB can't give delegates without arguments while C# can.
AddHandler CompositionTarget.Rendering, AddressOf RunSimulation
You will get a warning in VB that does not make much sense, and its solved with the EventHandler delegate.
The rendering can now be set at a specific time interval by coding the sub RunSimulation as follows:
Sub RunSimulation(ByVal sender As Object, ByVal e As RenderingEventArgs)
Dim rargs = DirectCast(e, RenderingEventArgs)
If ((rargs.RenderingTime.TotalMilliseconds - LastTimeRendered) > RunTimeInMS) Then
...
...
...
LastTimeRendered = rargs.RenderingTime.TotalMilliseconds
End If
End Sub
If you however want to just run it as fast as possible you could remove nearly everything in the previous section (This is what I have done in the code as I don't have
any problems with drawing it as fast as possible):
Sub RunSimulation()
...
...
...
End Sub
A DispatcherTimer could also be used, however it might be a bit more lagging, but it gets updated sequentially then too.
The TLM method is also a perfect candidate for parallel processing, since it loops through all the nodes twice in one timestep. To change it to parallel
is actually really easy, all you need to read is this MSDN documentation. Obviously it won't be faster
if you only got one processor in your computer, but I have 6 and it roughly changed the speed to the double of the previous code (others might even be a little faster, lets say with 8 CPU's).
The idea is to change the code in the function called CalculateNextTimeStep in the class TLMGrid. The cleared out function without parallel processing
looks like this:
Public Function CalculateNextTimeStep() As Decimal(,)
For x as Integer = 0 to Dimensions - 1
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
Next
For x as Integer = 0 to Dimensions - 1
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
Next
CurrentTimeStep += 1
Return CurrentPressure
End Function
To implement the code in parallel you firstly add the Imports System.Threading.Tasks you simply change the two for loops to this:
Public Function CalculateNextTimeStep() As Decimal(,)
Parallel.For(0, Dimensions, Sub(x)
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
End Sub)
Parallel.For(0, Dimensions, Sub(x)
For y As Integer = 0 To Dimensions - 1
...
...
...
Next
End Sub)
CurrentTimeStep += 1
Return CurrentPressure
End Function
And that's all there is to it, you should notice that in the original For loop it went up to Dimensions - 1, and in the Parallel.For
it goes until x >= Dimensions to you don't include -1. Sadly I couldn't process the code for the 3D viewer in this way, as other processes have a hold
on the variables used.
There is supposedly one more thing that you could do to improve the performance in WPF, that is to set down the frequency the screen gets updated, called Framerate.
This is set in the Application startup like so:
Imports System.Windows.Media.Animation
Class Application
Protected Overrides Sub OnStartup(e As System.Windows.StartupEventArgs)
MyBase.OnStartup(e)
Timeline.DesiredFrameRateProperty.OverrideMetadata(GetType(Timeline), _
New FrameworkPropertyMetadata() With {.DefaultValue = 60})
End Sub
End Class
By setting it down from 60, to 20 for instance, it should improve the CPU usage. I did not get any improvement so I left it at the default rate of 60, but you might try
and set this down to see what happens. I suspect that this applies to animations from storylines, but here I don't really know exactly what happens and how it works.
With this newfound speed one can now start to experiment with a different type of wave, namely the boat wake (picture from Wikipedia):
This could be modeled using a Cherenkov radiation, as first described to me
by Sergey Kryukov in this
question. The animated picture blow shows the principle:
It is just a series of pulses that moves faster than the wave can propagate. This is also explained very well in this video
from YouTube. This means that a bullet and a boat wave is mathematically quite similar. And the animation produced in the program using a series of Gaussian pulses then shows as follows:
It is amazingly how easy it was to create with the TLM model.
Moving around in the 3D world
There is also the point about how I implemented a very basic changing of the camera look direction, using the mouse. This is done the following way:
Dim StartPos, EndPos As New Point
Dim OriginalLookdirection As New Vector3D
Dim look As Boolean = False
Private Sub mainViewport_MouseDown(sender As System.Object, _
e As System.Windows.Input.MouseButtonEventArgs) Handles Me.MouseDown
StartPos = e.GetPosition(viewport3D1)
OriginalLookdirection = camMain.LookDirection
look = True
End Sub
Private Sub mainViewport_MouseUp(sender As System.Object, _
e As System.Windows.Input.MouseButtonEventArgs) Handles Me.MouseUp
look = False
End Sub
Private Sub mainViewport_MouseMove(sender As System.Object, _
e As System.Windows.Input.MouseEventArgs) Handles Me.MouseMove
If look Then
EndPos = e.GetPosition(viewport3D1)
Dim v As New Vector3D
v = OriginalLookdirection
v.X -= (StartPos.X - EndPos.X) * 0.1
v.Y += (StartPos.Y - EndPos.Y) * 0.1
camMain.LookDirection = v
End If
End Sub
As you can see I have just changed the X and the Y position along with a fixed amount of change in the X and Y of the look direction. This seems to give
a quite natural way of adjusting the view direction.
With some knowledge of vector calculus you can even move the camera position quite easily:
Private Sub mainViewport_KeyDown(sender As System.Object, _
e As System.Windows.Input.KeyEventArgs) Handles Me.KeyDown
If e.Key = Key.Up Then
Dim p As New Point3D
p = camMain.LookDirection
p = Vector3D.Multiply(p, 0.01)
camMain.Position += p
End If
If e.Key = Key.Down Then
Dim p As New Point3D
p = camMain.LookDirection
p = Vector3D.Multiply(p, 0.01)
camMain.Position -= p
End If
If e.Key = Key.Left Then
Dim r As New Vector3D
r = Vector3D.CrossProduct(Normalize(camMain.LookDirection), camMain.UpDirection)
r *= 0.2
camMain.Position -= r
End If
If e.Key = Key.Right Then
Dim r As New Vector3D
r = Vector3D.CrossProduct(Normalize(camMain.LookDirection), camMain.UpDirection)
r *= 0.2
camMain.Position += r
End If
End Sub
Private Function Normalize(ByVal v As Vector3D) As Vector3D
Dim result As New Vector3D
Dim length As Double = Math.Sqrt(v.X ^ 2 + v.Y ^ 2 + v.Z ^ 2)
result = v
result /= length
Return result
End Function
This however messes the easy look direction a bit, as the vector of the look direction is changed. You would have to adjust the length of the vector if you want this to work properly.
Then there is the Zoom.
Private zoomDelta As Vector3D
Private Sub Window_MouseWheel(sender As Object, e As MouseWheelEventArgs)
If e.Delta > 0 Then
camMain.Position = Point3D.Add(camMain.Position, zoomDelta)
Else
camMain.Position = Point3D.Subtract(camMain.Position, zoomDelta)
End If
End Sub
Where the ZoomDelta is dependent on the look direction and is set in the Window_Loaded sub.
Private Sub Window_Loaded(sender As System.Object, e As System.Windows.RoutedEventArgs) Handles MyBase.Loaded
Const ZoomPctEachWheelChange As Double = 0.02
zoomDelta = Vector3D.Multiply(ZoomPctEachWheelChange, camMain.LookDirection)
End Sub
That is a very quick run through of how I implemented the moving around while the simulation is running. Although it is quite laggy while the simulation is running
it functions reasonably well, given that the look direction is in a reasonable vector dimensions. I'm almost sure that there are better way of doing this,
but this works well enough for me at present.
From Vector to Pixels
Using a canvas and a background image is quite good, as I can separate the inputs, source walls placements from the actual calculation algorithm, so that it forms two separate layers.
This is especially good as you don't need to worry about handles for each pixel, but rather store the needed data in an array.
To move from the vector placement on the Canvas to the rasters, it is pretty straight forward for a point, as you could simply scale the point to the width and height,
and multiply with the raster length of your matrix. The code for converting a point on the canvas to a cell in the raster, is pretty straight forward:
Dim PointOnCanvas As New Point
PointOnCanvas = e.GetPosition(cnvTLM)
Dim RasterPointInMatrix As New Point
RasterPointInMatrix.Y = CInt(Math.Truncate(PointOnCanvas.X * Dimensions / cnvTLM.Width))
RasterPointInMatrix.X = CInt(Math.Truncate(PointOnCanvas.Y * Dimensions / cnvTLM.Height))
Oh, the values of X and Y is switched around, as the top left point is 0,0 and not 0, Height. The trouble starts when we define two points, and want them
to generate a continuous line in the matrix (or on a bitmap picture for that matter), one often resorts to Bresenham's line
algorithm. There are others available, but it is widely regarded as the standard implementation.
Given that it is so widely used one should imagine that it would be easy to find a working example around the web, but I was presumably mistaken,
as I couldn't find a ready to use code, that actually worked in all cases. (Some non working examples in VB can
be found here
and here). I didn't bother looking further and implemented a simplified version
of Bresenham's algorithm of my own instead, and it does not work perfectly, as I didn't bother with round off error control, but It should be rather easy to fix that!
Public Function GetBresenham_ishLine(ByVal P1 As Point, ByVal P2 As Point) As List(Of Point)
Dim DeltaX, DeltaY As Integer
DeltaX = Math.Abs(CInt(P1.X) - CInt(P2.X))
DeltaY = Math.Abs(CInt(P1.Y) - CInt(P2.Y))
Dim PositiveX, PositiveY As Double
If P1.X > P2.X Then
PositiveX = -1
Else
PositiveX = 1
End If
If P1.Y > P2.Y Then
PositiveY = -1
Else
PositiveY = 1
End If
Dim LineList As New List(Of Point)
If DeltaX > DeltaY Then
For i As Integer = 0 To DeltaX
Dim TempX, TempY As Integer
TempX = P1.X + PositiveX * i
TempY = CInt(P1.Y + DeltaY / DeltaX * PositiveY * i)
LineList.Add(New Point(TempX, TempY))
Next
Else
For i As Integer = 0 To DeltaY
Dim TempX, TempY As Integer
TempX = CInt(P1.X + DeltaX / DeltaY * PositiveX * i)
TempY = CInt(P1.Y + PositiveY * i)
LineList.Add(New Point(TempX, TempY))
Next
End If
Return LineList
End Function
For my purpose in this program it works good enough, as it is just there to show how the TLM code works. And with this simple thing you can behaviors of different settings.
The problems
This is really all that you need to generate the loss less TLM code, meaning that the only reason that the pressure drops with distance is that the energy is spread over a larger area. This is actual a quite normal way of viewing sound propagation over shorter distances, however with many reflections or at longer distances this might not give a correct view of the reality.
The TLM technique has some properties that are quite preferable to FDTD for parabolic equations, as FDTD are unstable at central difference, and only conditional stable for backward difference, while TLM is unconditional stable. The TLM does however have different problems.
- The wave from a point source does not expand in a uniform pattern, as it is a discrete technique
- Propagation speed is dependent on
1/Sqrt(2)
The code given does not include propagation loss, and to do this one would need two additional matrices, as each node would have a non-reflective tube which consumes energy. One would also have to implement a simple Gaussian filter to improve the results to make any professional use of it.
Other issues is that this is a Time depended algorithm, so to make use of the results, one usually needs to perform an Fourier transform at the receivers, and this is not included in the code as it is.
Faster methods of creating water effect's
The TLM method used here is designed to give an accurate mathematical description of waves, and there are faster options available if you are only interested in creating ripple effects like rain and boat wakes.
The approach is explained pretty well in one of the links provided by Sean, where the goal is just to create something that looks good enough, and is not exact mathematical accurate. The main idea is to use an area averaging technique of the cells around the current one.
History
The main goal of this article is to help you understand how one could write a TLM algorithm reasonably well, and to provide you with some ideas on how to implement a real world physics simulation run on the WPF engine.
The Gaussian pulse still needs some work, as it is now it doesn't correspond to well with changing frequencies, due to an incorrect implementation of the cut off frequency to the filter. There are also the issue of implementing a damping in the system, expressed in acoustics as attenuation by the coefficient Neper. A further development of elements that not only dampens the wave, but also allows for transmission is also quite useful, but requires some recoding in the wall section, as we would need 4 arrays for each of the walls instead of the one we currently have. The speed of sound can also not be changed in a uniform grid as this 2D simulation is using, so it has limited applications other that how to implement the TLM method reasonably fast in VB/C# using WPF.
The 3D TLM method is much more interesting if you want to use the results for
something other than to create raindrop effects, but then again it is not so easy to see what happens, as you need to look at each of the cross-sections. The model is not difficult to expand though, as you would only need to implement Up and Down, with the two corresponding Incoming and OutGoing pressure stubs. It should then be mentioned that there is a better way of constructing the grid than squares.
References
