# Sampling a 3D Sound Field

During the past six months or so I have been focusing on identifying the presence and directions of the individual sounds that make up an acoustic scene. My EigenScape recordings were recorded in fourth-order Ambisonic format, which contains more detailed spatial information than I used in my previous work.

The Ambisonic format makes it fairly simple to focus on sounds from specific directions in the scene using a process called beamforming, which is like pointing a virtual microphone. One common way of producing intensity ‘maps’ of an acoustic scene is to use beamforming to look in all directions and log the sound power in each direction to build up a picture of the scene. This is called the Steered-Response Power (SRP) method and is used in commercial ‘Acoustic Camera’ products.

Of course, this being the real world, we can’t steer the beam in all directions, just a representative set of discrete points. Doing this on a 2-dimensional plane isn’t that difficult. To sample a plane around a point, you go round in a circle, and it’s not difficult to divide a circle into equally-spaced points. In the diagram above, the sampling happens at eight points, but we can use any number of points and it will be possible to space them equally around a circle. It’s quite easy to confirm this with a simple thought experiment: Imagine a straight line of some length, and imagine dividing this line into some number of equally-spaced segments. If you then took that line and wrapped it around into a circle, those segments will always also form equally-spaced segments along the circumference of the circle. In other words, you’ll always be able to slice a pizza into equally-sized slices no matter how many guests you have. This means we can uniformly sample around a 2D plane at any arbitrary resolution. Uniform sampling around a circle.

Sometimes, sampling in 2D might be enough, but if we want a full picture of the entire acoustic scene in 3D (or indeed any 3D field viewed from a point – this is a common problem) we need to extrapolate from that circular sampling and sample round in a sphere. The question then becomes how to select points to uniformly sample across the surface of a sphere? This is a much thornier question than it initially seems.

Initially, the most intuitive way forward might be to approach this in the same way as we approached sampling the circle. We could achieve uniform circular sampling by placing points on a line and ‘wrapping’ it around, so for spherical sampling we could place points in a uniform grid on a 2-dimensional plane and then ‘wrap’ this around the sphere: Sampling a sphere using a regular grid.

Unfortunately, as you can see, when we use this regular grid and map it to the sphere, it results in the points near the poles ‘bunching up’, whereas points at the equator are more spread out. This results in us sampling the sphere in varying-resolution, with much more fine-grained detail at the poles than elsewhere. For instance, if we want a high level of detail at the equator, we end up with a ridiculously high density of points at the poles. Not only is this computationally inefficient, but can lead to the output power map being somewhat skewed towards the poles. It might appear that there’s more power in the polar regions simply as we are effectively sampling the same area over and over again, these points having equal weighting to those sampled in sparser areas (in simple models, anyway).

You’ll have seen this problem in reverse when looking at maps of the world. Land masses to the far north are stretched out and appear huge relative to those nearer the equator. Greenland, for instance, appears on a regular map to be about the same size as Africa, whereas in reality it’s only about 1/14th of the size.

So how do we fix this? Mathematically, there are actually only a very few distributions of points on a sphere that perfectly fit our criteria. These relate to the platonic solids, the regular 3D shapes made by regular 2D polygons with each vertex connecting the same number of faces. We can uniformly sample the surface of a sphere by sampling points corresponding to the vertices of these shapes! Unfortunately, only five of these shapes exist and these have relatively few vertices; the most we can get is 20, using the vertices of a dodecahedron, or equivalently the centres of the faces of an icosahedron. Uniform sampling of a sphere on the faces of an icosahedron / vertices of  a dodecahedron.

As you can see, this sampling scheme is very sparse and suitable only if we want low-resolution output. Unfortunately, it’s actually mathematically impossible to equally divide the surface of a sphere into any more points than this. So if you’re planning on serving spherical pizza to more than 20 very picky people, you’re out of luck!

So what do we do about it? We just have to settle for not-quite uniformity, and in this case good enough has to be good enough! One common approach to nearly-uniform spherical sampling is the icosphere. This is an extrapolation to a full sphere of the geodesic dome concept popularised by Buckminster Fuller as a way of constructing extremely strong dome structures for .

The method for constructing an icosphere is relatively simple. Starting with an icosahedron, we interpolate between the 12 initial points by splitting each triangular face into four smaller triangles by finding the halfway point between each pair of vertices and turning this into a new vertex. Repeating this process numerous times results in more densely-spaced points for higher-resolution sampling. Interpolating the triangular faces of an icosahedron.

Most of the code I found online to accomplish this was solid, but convoluted, involving elaborate methods for checking for pre-computed vertices (since all the triangles are adjacent, the interpolation results in lots of duplicate vertices), and requiring manual definition of which sets of vertices make up a face. I decided to write my own code within my preferred Python 3 / Numpy environment. This provided me with easier interaction with my other project code and allowed me to use some of Numpy’s features to make the code more elegant:

## Numpy Icosphere Implementation¶

In :
```# standard modules
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# distance matrix module
from scipy.spatial.distance import cdist

from utilities import normalise
```

First, I define the initial icosahedron using three intersecting orthogonal golden rectangles:

In :
```# golden ratio
R = (1 + np.sqrt(5)) / 2

vertices = np.array([[-1,R,0],
[1,R,0],
[-1,-R,0],
[1,-R,0],

[0,-1,R],
[0,1,R],
[0,-1,-R],
[0,1,-R],

[R,0,-1],
[R,0,1],
[-R,0,-1],
[-R,0,1]])

# plot
xu, yu, zu = vertices[:,0], vertices[:,1], vertices[:,2]
fig, ax = plt.subplots(1, 1, subplot_kw={
'projection':'3d', 'aspect':'equal'}, figsize=(10,10))
ax.scatter(xu, yu, zu, zorder=10, alpha=.6)
plt.show()
``` Then, I infer which sets of indices make up the faces of the icosahedron. This is not too difficult to do manually for the initial icosahedron, but the beauty of this code is that it can be used to infer faces on every level of interpolation. Scipy’s `cdist` module is used to get a matrix of distances between each point.

In :
```# find euclidian distances between all points -
# gives us a matrix of distances
euclid_dists = cdist(vertices, vertices)

# find list of adjacent vertices
sides_idx = np.where(
euclid_dists == np.min(euclid_dists[euclid_dists > 0]))

# concatenate output locations into one array
sides_idx = np.concatenate(
(sides_idx.reshape(-1,1), sides_idx.reshape(-1,1)), axis=1)

# remove duplicate sides_idx (there are many)
_, idx = np.unique(np.sort(sides_idx), axis=0, return_index=True)
sides_idx = sides_idx[idx]
```

`sides_idx` now contains a list of pairs of vertices that make the sides of the icosahedron.

In :
```sides_idx
```
Out:
```array([[ 0,  1],
[ 0,  5],
[ 0,  7],
[ 0, 10],
[ 0, 11],
[ 1,  5],
[ 1,  7],
[ 1,  8],
...]], dtype=int64)```

Then, I work out which of these sides connect to make faces:

In :
```# set up empty array
faces_idx = np.array([], dtype=int)

for i in np.unique(sides_idx[:,0]):
# extract sides_idx related to each vertex
a = sides_idx[np.where(sides_idx[:,0] == i),1]

for j in a:
for l in j:
# find 3rd adjacent vertices common to both points
b = sides_idx[np.where(sides_idx[:,0] == l), 1]
intersect = np.intersect1d(a,b).reshape(-1,1)

for m in intersect:
faces_idx = np.append(faces_idx, np.array([i,l,m]))

# output is a 1D list, so we need to reshape it
faces_idx = faces_idx.reshape(-1,3)
```

This outputs a list of the triads of points that make faces:

In :
```faces_idx
```
Out:
```array([[ 0,  1,  5],
[ 0,  1,  7],
[ 0,  5, 11],
[ 0,  7, 10],
[ 0, 10, 11],
[ 1,  5,  9],
[ 1,  7,  8],
[ 1,  8,  9],
...]], dtype=int64)```

Plugging these indices into our actual vertex co-ordinates gives us a face co-ordinates list (a three-dimensional numpy array):

In :
```# 3D matrix with xyz co-ordnates for vertices of all faces
v = vertices[faces_idx]
v
```
Out:
```array([[[-1.        ,  1.61803399,  0.        ],
[ 1.        ,  1.61803399,  0.        ],
[ 0.        ,  1.        ,  1.61803399]],

[[-1.        ,  1.61803399,  0.        ],
[ 1.        ,  1.61803399,  0.        ],
[ 0.        ,  1.        , -1.61803399]],

[[-1.        ,  1.61803399,  0.        ],
[ 0.        ,  1.        ,  1.61803399],
[-1.61803399,  0.        ,  1.        ]],

...]]])```

We can then find the midpoints of each side all in one go to generate the smaller triangles by using numpy’s `roll` function. This moves each element of an array along one space, reintroducing elements that fall off the end back at the beginning. By adding the rolled version of the `v` array to the original version and dividing by 2, we calculate the midpoints of each pair of points in one go:

In :
```midpoints = ((v + np.roll(v,1,axis=1)) / 2).reshape(-1,3)
```
In :
```midpoints
```
Out:
```array([[-0.5       ,  1.30901699,  0.80901699],
[ 0.        ,  1.61803399,  0.        ],
[ 0.5       ,  1.30901699,  0.80901699],
[-0.5       ,  1.30901699, -0.80901699],
[ 0.        ,  1.61803399,  0.        ],
[ 0.5       ,  1.30901699, -0.80901699],
[-1.30901699,  0.80901699,  0.5       ],
[-0.5       ,  1.30901699,  0.80901699],
...]])```

It is then simply a case of adding the new midpoints to the list of vertices and removing duplicates, and we have the points of our interpolated icosahedron:

In :
```# add new vertices to list
vertices = np.append(vertices, midpoints, axis=0)
# find duplicate vertices
_, idx = np.unique(vertices, axis=0, return_index=True)
# remove duplicates and re-sort vertices
vertices = vertices[np.sort(idx)]

# normalise points to surface of unit sphere
vertices = normalise(vertices, axis=1)

# plot
xu, yu, zu = vertices[:,0], vertices[:,1], vertices[:,2]
fig, ax = plt.subplots(1, 1, subplot_kw={
'projection':'3d', 'aspect':'equal'}, figsize=(10,10))
ax.scatter(xu, yu, zu, zorder=10, alpha=.6)
plt.show()
``` And there you have it. To calculate further levels of interpolations to get more vertices, you simply use the output of the first interpolation as the input in place of the icosahedron at the beginning of this process and the whole thing can be looped. It’s important not to normalise the points to the unit sphere if you intend to proceed to higher orders of interpolation as this causes the inference of faces stage (`cdist` etc.) not to work properly. Always normalise at the end of the process only.

My function that automates this whole process can be found in https://github.com/marc1701/area-beamforming/blob/SRP_dev/spherical_sampling.py. With this function you can specify your desired level of interpolation, whether or not to receive output as cartesian or spherical co-ordinates, and also whether you’d like to receive points specifying the faces of the interpolated icosahedron rather than the vertices (similarly to how a dodecahedron is the ‘dual’ of an icosahedron, with faces and vertices interchanged).

In :
```from spherical_sampling import geodesic

vertices = geodesic(3, co_ords='cart')

# plot
xu, yu, zu = vertices[:,0], vertices[:,1], vertices[:,2]
fig, ax = plt.subplots(1, 1, subplot_kw={
'projection':'3d', 'aspect':'equal'}, figsize=(10,10))
ax.scatter(xu, yu, zu, zorder=10, alpha=.6)
plt.show()``` I think this code is an elegant solution for easily calculating nearly-uniform sampling points on a sphere, and should definitely be uniform enough for most purposes. In fact, the Eigenmike itself uses a variation on this for the positioning of its microphones – the capsules are placed on the vertices of both the dodecahedron and icosahedron for 32 nearly-uniformly spaced sampling points.

## Fibonacci Sampling

I was quite happy with this code for several minutes before a colleague of mine threw a spanner in the works. He introduced me to sampling a sphere using a spiral based on the Fibonacci sequence, which is used by certain plants to evenly distribute seeds or leaves. The Fibonacci sampling method allows for the nearly-uniform distribution of any arbitrary number of points. This overcomes the main problem with the icosphere method, which is that the number of sample points, and thus the sampling resolution, is limited to the number of points that result after each interpolation of the icosahedron: 12, 42, 162, 642, 2562, etc. These values increase exponentially on each interpolation, limiting your options if you had a specific number of sampling points or resolution in mind.

In two-dimensions (as used by plants), the spiral is formulated as follows: Where λ and r represent the angle and radius in spherical co-ordinates, Φ is the golden ratio (about 1.618) and n is the index for each point up to a maximum of N points. For 500 points, this results in a pattern that looks like this: This should look pretty familiar to anyone who’s ever taken the time to look at the seeds in a sunflower. From here, it’s a pretty simple step to take this pattern into 3D on the surface of a sphere. We need to replace the radius calculation with one that will calculate the elevation of each point along a unit sphere, essentially spreading this pattern downwards over the sphere:  Sampling a sphere using the Fibonacci spiral method.

It’s a little more difficult to see the unique spiral distribution on the sphere here, but it can be seen that this is a uniform distribution relative to the regular grid (check out this cool animated version by KrazyDad for a clearer look). It can be calculated with much simpler code and massively less processor overhead than the geodesic icosphere method. I think the Fibonacci method is easily the most intuitive way of getting a simple nearly-uniform distribution of points on a sphere and should be suitable for most purposes (I’m not going to get into the nightmarish t-designs!).

For these reasons, I’ve been using the Fibonacci pattern in my work so far and it’s proven to be effective for making steered-response power maps. All of my code for spherical sampling can be found on my project github repo.