Air absorption of sound as a digital filter – Part 1: Theory

There is a very important future of sounds coming from a distance besides their obvious decrease in loudness: they sound muffled. Thunderstrike is the best example: when striking close to us, we hear an ear-piercing, high-pitched crack, but when afar, we hear a low-pitched, thumping, growl. This is an effect of sound absorption by the atmoshphere. The air acts as a low-pass filter, and this filtering is the “stronger”, the bigger the distance and higher the frequency. It is also dependent on the air’s temperature, humidity and pressure.

This future is one of the properties crucial for the brain to estimate the distance to the sound source (the other being level of the sound). Since it’s fundamental to how we perceive sound , it cannot be ignored in a realistic auralization. Sound coming from a distance that is rich in high frequencies , sounds unnatural. This is a well known fact and must be accounted for if one wants to achieve audio immersion. Most common and also simplest way to do that is to apply a LP filter on the sound source, with cutoff frequency changing according to the distance. More time-consuming way would be to record the sound from the distance it is meant to be played. Although this gives probably the most realistic effect (because it also includes attenuation over ground and reflections from it), it is also very unpractical. What if there is a sample that will be played from various distances in regard to a player (i.e. gunshots) ? Or the weather (in the game) changes? One could obviously record all the samples in all possible conditions…but this is nowhere near universal.

We can instead hit somewhere in between, which is to create a filter that will be basing on the mathematical model of the sound absorption by air. Such a model exists, moreover, it is an ISO standard. Given all the equations needed, we can get the set of absorption coefficients that we will use as our filter. This approach brings a very big benefit: since all it requires is an absorption model and a distance the sound traveled, it can be reused everywhere where custom DSP is possible, own sound engine being the best example.

So, let’s get to it. In the rest of the post i will talk about the model and filter design in this particular case. This will be rather some basic DSP. No coding though, actual implementation will be covered in the next post.
Be wary! Altough i’ll try to explain most of the stuff, signal processing is not something one can learn in a day. It would be impossible to explain all the technicalities, how and why they work, in this post. If You have some basics in filter design, there will be nothing new for You. If You don’t, then I strongly suggest checking out the links on the bottom, first.

1. The model
A set of equations that i’m going to use can be found here: LINK.Those are equations that can be found in ISO 9613-1 “Attenuation of sound during propagation outdoors”. I won’t rewrite them here, just briefly go through them. The most important formula for us is an analytical way to calculate the absorption coefficient, a, measured in dB/m. As we can see, it’s a function of sound frequency, temperature, pressure and two frequencies, called “relaxation frequencies” of molecular oxygen and nitrogen (don’t ask, don’t know) that are also functions of temperature, pressure and humidity. So we need only three things to calculate an absorption coefficient for a given frequency: temperature, humidity and pressure. To calculate the attenuation over the distance, we simply multiply the absorption coeffient by that distance: A (r) = a*r [dB]. The result is the decrease in sound’s pressure level, in dB.
Here are some plots of change in level (attenuation*-1) for frequency range 0-2kHz and different distances (0dB is sound (pressure) level of the source):
att1att2att_0C_10%att_0C_90%

2. Designing the filter

Ok, so now that we are able to calculate air absorption coefficient for any frequency we desire, let’s move to designing our digital filter. I will be using a method called “Frequency sampling” (details below) and our filter will be of FIR (Finite Impulse Response) type. This is the simplest, yet, a very popular way to design digital filters mainly because of how easy it is to implement. Basically, we are creating samples of the desired frequency response and then, using inverse discrete fourier transform (IDFT) on those samples we obtain filter’s coefficients ( it’s impulse response). FIR filter comes with several benefits: they require no feedback and therefore are always stable, introduce smaller rounding errors and can be “coded” to be very efficient using FFT (more on that later). They also have a very neat future: they can be easily designed to have a linear phase. Linear phase means that every frequency component passing through a filter will be shifted (delayed) by the same amount – consequently there will be no phase distortions. Whether such phase distortions are audible is a very controversial topic, with many claiming that they aren’t audible in most cases. Nonetheless, the conclusion from many conducted experiments is that linear phase should be desired, where possible, especially in the audio domain. Since it is so easy to achieve in case of FIR filters, why not take the opportunity? (For those interested about the audibilty of phase distortions, I posted few links that elaborate on the subject).

Ok, after that little explanation of my choices, let me now lay the steps of our algorithm of designing our filter and then walk you through the details on each one.

1. Get the distance that the sound wave travelled.
2. Choose the length of the filter. This must be a power of 2 (to take advantage of very fast FFT).
3. Calculate air absorption coefficient for every required frequency in the [0 , sampling_rate/2 ] range.
4. Multiply array of absorption coefficients by the distance. This, multiplied by -1 is our ideal filter amplitude response.
5. Convert the response from the decibel scale to linear one.
6. Mirror the array with respect to the N/2 sample. This is required if we want the result of the IDFT to be real-valued, not complex.
7. Perform a N-point IDFT (using FFT) on the array. The result is an array of filter’s coefficients/impulse response.
8. Shift the impulse response by N/2.
9. Optionally, multiply by a window function.

Ad 1. We will be creating a new filter everytime the distance changes.

Ad 2. The longer the impulse response, the more our filter will resemble the ideal we design (it’s actuall frequency response will be closer to the designed one) …obviously at the cost of computational power. Since we will be taking an advantage of FFT (Fast Fourier Transform) in it’s fastest implementation (Radix-2) to perform IDFT on our designed frequency response array, it must be of length N = 2^p, for example N = 256 or 512.

Ad 3. Domain of a product of Discrete Fourier Transform is a vector of complex spectrum values that correspond to frequencies spaced in interval equal to fs / N, where fs is the sampling rate and N is the length of signal. For example: if we performed DFT on a signal of length N = 6 and with fs = 10, the result would be a vector of values corresponding to: 0, 10/6, 2*10/6, 3*10/6, 4*10/6, 5*10/6 [Hz]. We can see that the Nyquist frequency (fs/2, in this case 3*10/6 = 5Hz) is the (N/2 + 1)’th sample.
In order to obtain a real-valued impulse response, we need to make our vector symmetrical. Therefore we are sampling values only in the [0, Nyquist] range, that is for first N/2+1 samples -> we calculate absorption coefficients only for f = 0, fs/N, 2*fs/N…fs/2 Hz.

Ad 4-5. Since the attenuation over the distance r is equal to absorption_coefficient*r , we need to multiply every coefficient by the current distance in meters. If we then multiply by -1, we actually get our ideal filter’s amplitude response (in the 0…fs/2 range). Here is why: Filter amplitude/magnitude response (often denoted by A(f) ) is an information of how much the given frequency component will be changed in amplitude. It is an absolute value of a product of DFT done on the filter’s impulse response (magnitude of a complex number, hence “filter’s magnitude plot”). It is a function of frequency with values >= 0. Those values basically are the ratio between magnitude of a frequency component (sinusoid) at the output, to it’s magnitude at the input. Therefore a value of 1 means no change in amplitude, 0-1 means attenuation and a value > 1 means that the filter will amplify this frequency.
In a logarithmic (sometimes called decibels) scale, this is given by:

L(f) = 20*log10( A(f) ) [dB]

where log10 is a base-10 logarithm. L = 0 dB means no change in amplitude (20*log10(1) = 0), L = -Inf means total attenuation (logarithm goes to -infinity for x->0) and L > 0dB means amplification. How is this all related? If we multiply our vector of absorption coefficients from step 3 by the distance r, we get a result in dB, which is a vector of attenuations for each frequency -> Att(f) = 20*log10( p(0) / p(r) ), where p(r) is the sound pressure at distance r and p(0) is the initial pressure. Since amplitude response is a ratio of output/input, and our input is p(0) and output is p(r), our L(f) should be 20*log10( p(r) / p(0) ) which is equal to -20*log10( p(0) / p(r) ) = -Att(f)! Now, we only need to convert Att(f) to linear scale so we can perform IDFT on it. Conversion back from decibel scale is given by: A(f) = 10^( – Att(f)/20 )

Ad 6. This is a very important property of the Discrete Fourier Transform and the essence of the method we use to design our filter. If an impulse response caluclated with IDFT is not to be complex we must assure that: A(n) = A(N-n), where n = 1…N/2 if zero-indexed. Obviously we don’t want the impulse response to be complex, because then the filtered audio signal would be complex too. And how would that sound? :)

Ad 7. Since length of our vector is the power of 2, we can use a Radix-2 FFT (or rather Inverse-FFT) algorithm to calculate our filter’s impulse response.

Ad 8. If the impulse response is symmetrical: h(n) = h(N-n) or assymetrical: h(n) = -h(N-n) the filter has a linear phase (phase is a linear function of frequency). However, in our case it’s not quite the case: you may noticed that after shifting by N/2 the impulse response is symmetrical except the rightmost sample. That is the downside of this method of designing filter. We have two choices: we can either accept the fact that for lower order of the filter (order = length – 1) the phase will be a “little” non-linear or we can discard the rightmost sample which will guarantee linear phase although the resulting frequency response won’t be interpolated through the points we designed. However, for higher filter lengths, like 128 or more, phase will be very close to linear if we choose to leave it as it is anyway and that still can be improved by multiplying by a window (proof on the graphs below)

Ad 9. That is an optional step , but highly recommended since multiplying our impulse response by a smooth window function (i.e. Blackman window) comes with big benefits. Main one is better interpolation of frequency response between the points we sampled, since by default it is interpolated with a sinc function (since that is the fourier transform of a rectangular window and that’s what every discrete signal is multiplayed with, being of finite length). Another benefit is the more linear phase. To visualize those benefits i created two filters of different lengths using our algorithm with impulse response multiplied by a blackman window and not (click to enlarge):

N = 64:
N = 64

N = 128:
N = 128

Both benefits, frequency response more similar to our design and linear phase are clearly visible.

Now, we have all that is required to…filter any type of sound data with our air-absorption-filter :) Filtering itself with an example of implementation, in the Unity engine, in the next post!

——
LINKS:

Audibility of phase distortion:
http://www.silcom.com/~aludwig/Phase_audibility.htm
http://www.audioholics.com/room-acoustics/human-hearing-phase-distortion-audibility-part-2

DFT / FFT:
A very good paper on DFT for everyone: http://www.analog.com/static/imported-files/tech_docs/dsp_book_Ch8.pdf
http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf

Frequency sampling method:
http://www.bulletin.zu.edu.ly/issue_n15_3/Contents/E_04.pdf

Wiki has a great article on window functions:
http://en.wikipedia.org/wiki/Window_function