White is zero pressure, toward red is positive, towards blue is negative.

(click title to see video)

So I made a deal with Gry; I will challenge Odeon with any complex geometry I wish, and if it instantly works as promised I will write her a testimonial. After lunch I felt creative, and so I designed some strange looking building, which can be at best characterised as some twisted hybrid between a museum, concert hall and the fuselage of a 747. The world should be thankful that I am not an architect.

Back to our experiment, I’ve really invested a lot of energy (not easy after a conference lunch) in making the geometry complex and confusing to human beings and machines alike. The model was finished, and guess what, at a click of a button Odeon imported it. A cup of tea and a few clicks later, the model was up and running with wall materials and multiple sources. So, as I am a man of my word, here is my testimonial: Odeon handles complex geometries like no other acoustic modelling software I have encountered; it is quick, simple to use and most importantly troublefree. Very well done.

Jon.

But I do have a good excuse this time, as I’ve been preparing for the EAA’s conference “Forum Acusticum 2011″ that will take place next week in Aalborg, Denmark.

I’ll be delivering a paper on a simple method to solve frequency dependent problems in numerical time domain methods within one of the sessions on computational acoustics.

My presentation is scheduled to be on Friday 1-Jul at 8:20.

If you’re around – come say hi!

The rest of the conference programme is available here.

With a little help from my home-brewed FDTD model and the (very good and very free) visualisation platform ParaView, I have created some 3D animations of sound wave propagation. In this demonstration you can see a wideband acoustic monopole (Gaussian pulse) propagating in a small (6x6x6 meters) room. In the middle(ish) of the room I’ve put a small panel-shaped fully-reflecting phase-reversing (R=-1) obstacle  so it is evident how diffraction occurs. The boundaries of the room are somewhat reflective (Normalised impedance = 9).

A more thorough review to follow right after the conference.

The tutorial is very basic and intended for students, scientists and engineers with background in acoustics and mathematics but without any experience in numerical time domain methods. The text is written in an informal language, so please excuse me. As such, I won’t reference anything directly in the text; however, more formal and elaborate discussions on the topic can be found in many scientific papers and some books (see Bilbao and some chapters in Smith III).

The FDTD method can be applied to many partial differential equations. For room acoustics it  is usually used to solve either Euler’s equations of continuity (mass and momentum), or the second order wave equation. The main idea is to sample our problem in space and time and solve it iteratively. Spatial samples are called nodes or junctions (depending on who you ask…) and are interconnected in different schemes or topologies (depending on who you ask…) to form a grid or mesh (depending on who you ask…). These schemes refer to the way nodes are interconnected on the grid, forming differently shaped cells, such as cubic, interpolated-cubic or octahedral, amongst others. If you want an overview read this post and its corresponding paper.

In this post we shall derive the Standard Leapfrog Scheme (hereinafter SLF) for the 2nd order acoustic wave equation. This scheme is often referred to as the Kirchhoff Digital Waveguide Mesh (or K-DWM). The reasons for this choice are as follows:

1. The 2nd order wave equation results in a non-staggered solution, which is very easy to understand.
2. The SLF is probably the most common scheme in room acoustics (although far from being the most computationally accurate).
3. The SLF scheme is very computationally efficient, so we can easily implement it in Matlab.
4. The SLF scheme is very easy to derive.

Don’t get me wrong. The SLF scheme is far from perfect, and nowadays I hardly use it. Nevertheless, it is an ideal starting point for understanding FDTD.

The 2nd Order Wave Equation

(The following derivation is mainly based on Bilbao).

Let’s have a quick look at our governing PDE:

Also bear in mind that , which is our sought quantity – the pressure at a specific point in space at a specific time.

If we can numerically approximate the partial derivatives in the wave equation, and discretise space and time, we can come up with an algorithm that will allow an iterative solution of our governing PDE. Let’s do just that.

Finite Difference Approximation

Recall that by the fundamental theorem of calculus, we can approximate the derivative of a function at a point by means of finite differences. A 2nd order finite difference is defined as:

h-squared is the distance between the two points on our function, so for the time derivative it corresponds to the temporal period (delta-t), and for the space derivatives it corresponds to the spatial period (delta-x, delta-y and delta-z respectively).

Applying that to the four partial derivatives in the wave equation (one for time, three for space) we can show that:

These will soon become very useful…

Derivation of a Discrete Update Equation

Now, before we proceed there is one more ‘preparation’ to make. We need to normalise the temporal and spatial discretisation steps to integer quantities. We do that by introducing four discrete indices corresponding to their continuous counterparts: n for discrete time, and i,j and k for discrete space, such that:

Now, we can rewrite the approximations to the derivatives using the above discrete notation and substitute back into the original wave equation,. yielding:

Next, we have to consider two more things. First, if we assume uniform discretisation in all spatial directions (i.e. delta x = delta y = delta z) we can simplify things a whole lot, so let’s do that:

Also, we want to consider the relationship between change in time, change in space and velocity of propagation. In a normal physical situation, the change in space would equal the change in time, times the velocity of propagation; however if we try to do that here some nasty things would happen. This is because our algorithm is recursive, and a growing solution will render it unstable. As such, we adhere to the Courant-Friedrichs-Lewy criterion, which for our scheme states that:

Where D is the dimensionality of the problem.

So by plugging that into our difference equation we get:

Solving for , the equation can be rearranged to show that:

Which is the update equation for the SLF scheme, typically shown in scientific publications on this topic.

The equation can be further simplified, if we set the Courant variable to the scheme’s upper stability limit, which for a 3D SLF is . Then we get:

Which in words mean that the pressure on any node on the grid one sample ahead in time, equals a weighted sum of its surrounding pressures minus its own pressure one sample back in time… Simple, right?

This equation can help us to solve for sound propagation in the medium, which is assumed to be air. At the limits of the domain the update equation becomes slightly different, as it should account for the applied boundary conditions. But let’s not be bothered about that for now.

In ‘part 2′ of the tutorial, we’ll walk through the implementation of this in Matlab. Stay tuned!

Another very detailed model (shown in the image above), is kindly provided by “Ten24“. It is a full head and semi-torso laser-scanned to a high resolution poly mesh – very rare quality for a freebie.

GP-GPUs have been around for quite some time, and are a very cost effective solution for executing parralel algorithms. For example, a high end consumer graphics card (around 400£) has about 240 cores and 1.5GB of memory. There are also specialised products dedicated to scientific computing, such as nVidia’s Tesla systems. The nice thing about nVidia’s architecture, termed CUDA, is that you don’t need to be a graphics programmer in order to handle it.

Truly, CUDA is simple to learn but difficult to push to the limit, however even simple algorithms run very well and getting to know the platform is definitely worth your time (given that your algorithms are suitable for parralel computation). I guess that the most difficult thing is getting used to the GPU programming paradigm; the rest is just small technical details.

I am not a computer scientist, so if any of you readers are, please forgive my choice of words. In traditional CPU programming we tend to think of processor activity as a ‘costly’ operation wheras memory traffic is treated as ‘cheap’, so programmers tend to store a lot of information in short term memory rather. In GPU programming, the opposite happens. Memory access is slow but computation is fast, so we prefer to re-compute rather than pre-compute. That’s the main idea, and if you can put your head around it, you’ve already made a big step towards a good GPGPU algorithm.

This raises some new issues, and calls for rexamination of methods that were ‘abandoned’ due to their high computational expense. One of these is the pseudo-spectral time-domain (PSTD) method (see Liu or Spa et al.) which strongly rely on calculating the partial derivatives in the wave equation by means of Fourier transforms – an ideal operation for GPGPU parralelisation. My co-supervisor, Professor Jamie Angus has written about this in a recent AES paper.

Back to the original scope of this post,  ‘traditional’ FDTD methods also benefit from GPGPU technology. In a recent study I conducted with my supervisor Dr. Bruno Fazenda, we have shown that up to x30 improvement in calculation speed can be achieved. Other studies (for example see Savioja, or Southern et al. ) have shown similar improvements. As I mentioned earlier, this is definitely ‘worth your while’.

WaveCloud is the title of my ‘home-brewed’ GPU-enabled FDTD model, that makes use of the recently published compact explicit schemes. It consists of a core engine written in plain C with CUDA extensions, and a wrapper/GUI written in C#. Here’s a screenshot of ‘WaveCloud Desktop’:

In case you were wondering, there are no commercial plans for this software and it is meant to remain within the University of Salford. However, I am now working on a MATLAB toolbox utilising the same engines, that will most likely be released for research purposes in the future (another reason to follow this blog…).

A unique feature of this toolbox is the binaural processors, currently supporting Lindemann’s model; and hopefully also Gaik and Breebaart soon (from what I hear).

Note that you also need the LTFAT toolbox to run these functions.

Here’s an example a Gaussian pulse being scattered from a round object:

(yes. I know. The object is phase reversing…)

The rest of the simulations can be seen at Jon’s Acoustics Tube.

Enjoy,

Jon.