This is the first in a series of three blog posts about the basic use of Acoular. It explains some fundamental concepts and walks through a simple example.
Acoular is a Python library that processes multichannel data (up to a few hundred channels) from acoustic measurements with a microphone array. The focus of the processing is on the construction of a map of acoustic sources. This is somewhat similar to taking an acoustic photograph of some sound sources.
To use Acoular, we first have to import Acoular into our notebook.
import acoular
In this example, we want to analyze time histories from 64 microphones stored in the file "three_sources.h5". The file contains data from a measurement (actually a simulated measurement in this special case) of scene with three different sound sources.
This file is in HDF5 format, which is an open all purpose numerical data container file format. Besides the time histories it contains information about the sampling rate. To learn about the internal structure of the file, have a look at the file! You will need an HDF5 file viewer (e.g. https://www.hdfgroup.org/downloads/hdfview/). For now, we can skip this step and simply load the file into our notebook. The file can be found in the same directory as the notebook itself.
ts = acoular.TimeSamples( name="three_sources.h5" )
The file is not directly opened using Python commands, but instead we use the TimeSamples class from Acoular. This class manages the data in the file in an intelligent way which even allows to use huge data files that do not fit into the memory.
If we inspect the ts object, we see not the data itself, but just the type and the location in memory.
ts
We can use this object now to answer some questions about the data:
- How many channels,
- and how many samples do we have?
- What is the sampling frequency?
print( ts.numchannels, ts.numsamples, ts.sample_freq )
The signal processing can either take place in the time domain or in the frequency domain. To work in the frequency domain, the time history data must be transformed into power spectra. More specifically, we need the cross spectral matrix (CSM) which contains the pairwise cross spectra of all possible combinations of two channels. The cross spectral matrix is computed using Welch's method. For this, an FFT block size and a window type have to be chosen.
ps = acoular.PowerSpectra( time_data=ts, block_size=128, window="Hanning" )
ps.fftfreq()
We see that after the FFT we do have spectra with 400 Hz frequency spacing. For most applications this would be too coarse. We use it here to get faster processing. If we wish to have finer spacing, we need larger block sizes.
Up to now we have defined how the processing should be done, but we did not compute the actual CSM. Nearly all expensive computations in Acoular are only done on demand using a lazy evaluation paradigm. We can trigger the computation of the cross spectral matrix by just asking for it. In this example we do want to print this already large matrix. Instead we only print its shape.
ps.csm.shape
This matrix actually has the dimensions 65 (number of frequencies) times 64 by 64 (number of microphone channels).
To continue with our task to generate an acoustic photograph, we need the microphone positions. In Acoular, one option is to read them from an XML file. The file looks like this:
<?xml version="1.0" encoding="utf-8"?>
<MicArray name="array_64">
<pos Name="Point 1" x="0.152" y="0.1141" z="0"/>
<pos Name="Point 2" x="0.134" y="0.1021" z="0"/>
...
<pos Name="Point 64" x="0.0218" y="-0.0329" z="0"/>
</MicArray>
(most of the lines are omitted here)
A MicGeom object handles the file:
mg = acoular.MicGeom( from_file="array_64.xml" )
mg.mpos
The mg object now contains all information about the microphone positions. Now let us this to plot the microphone geometry. For plotting we use the matplotlib package.
%matplotlib notebook
import matplotlib.pylab as plt
plt.plot( mg.mpos[0], mg.mpos[1], 'o' )
plt.axis( 'equal' );
This gives a nice impression of the microphone arrangement which consists of 7 inertwined planar logarithmic spirals (hard to see the spirals, I admit).
To map the sound, we need a mapping grid. Here we construct a simple regular and rectangular grid. Note that we have to decide about the size, spacing and the distance (z coordinate) from the array. Size and spacing define the number of points in the grid. The more points in the grid, the better the resolution of the acoustic photograph, but the longer the processing will take.
rg = acoular.RectGrid( x_min=-0.2, x_max=0.2,
y_min=-0.2, y_max=0.2,
z=0.3, increment=0.01 )
rg.pos()
The rg object now has all information about the grid.
The actual method we will use is beamforming. Basically this works by using the combined microphones as a sound receiver with a directivity that is steered to each one of the grid points in turn. One important element in beamforming and similar methods is the steering vector implemented in Acoular in the SteeringVector class. This vector "connects" grid and microphones and takes into account the environmental conditions. These conditions are defined by the speed of sound and a possible background flow. If not set explicitely, a 'standard' environment is created in Acoular which assumes quiescent conditions (no flow) and a speed of sound of 343 m/s (air at 20°C).
st = acoular.SteeringVector( grid=rg, mics=mg )
st.env.c
Indeed the standard speed of sound is used.
Now, we define the method we want to use and set up a standard (basic) beamformer. For this need two ingredients: the CSM and the steering vector.
bb = acoular.BeamformerBase( freq_data=ps, steer=st )
Remember that Acoular uses lazy evaluation. No computation yet!
This means although we set up everything needed to perform beamforming, computation is postponed until the results are actually needed.
This will happen if we ask for the result. In this example we are interested in the sum for of all FFT frequency lines in the 3rd octave band 8000 Hz. This means we need the maps for all these frequencies and then add them together to synthetically "mimic" the result of a third octave filter for that band. The result is given as mean square sound pressure contribution at the center location of the microphone array. In Acoustics, we are usually like to have the results in the form of sound pressure levels (SPL). Acoular has a helper function L_p to compute the levels from the mean square.
pm = bb.synthetic( 8000, 3 )
Lm = acoular.L_p( pm )
The map is now stored in the array variable pm and the array Lm holds the soundpressure levels computed from this.
print(Lm.shape)
print(Lm)
The map has the same dimensions (41 x 41) as the grid. Any zero result in the map will be clipped to -350 dB level instead of -infinity.
Now lets plot the map as a color-coded image with 15 dB dynamic range between both ends of the color scale.
plt.figure()
plt.imshow( Lm.T, origin="lower", vmin=Lm.max()-15, extent=rg.extend() )
plt.colorbar();
Now we enjoy the acoustic photograph of the three sources. Although it looks a bit blurry, we can guess the location of the sources and see as well that the sources have different strength.
This post has demonstrated how to use Acoular to produce a sound map or acoustic photograph with beamforming. Here you can see how to change some parameters and here you can learn about Acoular and time domain processing.