Echolocation Part 1: The Mathematics

QR Code linking to this page

Written by David Evans, first published (last updated )

Echolocation (active SONAR) is the process of creating a sound then listening for its echoes to determine the location of an object (such as how far away it is, or a more complete position through triangulation). Several species of animal (most notably bats and whales) use echolocation for navigation and hunting, and can not only determine 3D locations of objects, but also classify the object (getting information about its texture, density, and size).

Various types of specialist hardware also use echolocation; either for navigation and mapping (such as ship SONAR), or for distancing (such as security systems). However it is rarely used in general-purpose systems like laptops and phones due to hardware-related uncertainties and a lack of practical applications. In spite of this, it is still an interesting challenge to implement echolocation for general-purpose hardware.

This article is the first in a 3-part series:

Principle

At its most basic, echolocation consists of:

  1. Make a sound;
  2. Listen for the echo of the sound;
  3. Use the measured time between (1) and (2) to estimate the distance of the object.
Basic principle of echolocation.

Since the sound has to travel out from speaker to object, then back to the microphone, the distance can be calculated as approximately 12×time taken×speed of sound\frac{1}{2} \times \text{time taken} \times \text{speed of sound} (assuming speaker and microphone are relatively close together).

A sound wave travelling from a speaker, bouncing off a wall, and being reflected back to a microphone

The speed of sound in air can be approximated as 343m/s ±5%(note 1), so if our microphone picks up the sound 5ms after the speaker made it, the soundwave must have travelled about 1.72 metres in total (round-trip), meaning the object is 86cm away (±4cm)(note 2).

This equation can also tell us how granular our answer will be: if we can determine the time to within a millisecond, then we can calculate the distance with ~17.2cm granularity (plus the overall 5% margin of error from variation in the speed of sound). Most computers use an audio sample rate of either 44.1kHz or 48kHz, so if we can get to-the-sample precision in our timings (~0.02 milliseconds), we can expect a distance granularity of 3.9 or 3.4mm, plus the overall 5% margin of error(note 3).

Dedicated hardware

Hardware which performs this exact operation is readily available, such as the HC-SR04, and typically advertises accuracy of “up to 3mm” (theoretically requiring a sample rate of at least 60kHz).

HC-SR04
An HC-SR04 Ultrasonic Distance Sensor. It has a directional transmitter (speaker) on the left and a directional receiver (microphone) on the right.

This specific example sends out a series of 8 very short 40kHz beeps (which requires a digital sample rate of at least 80kHz due to folding). The beeps are inaudible to humans, since the human ear is unable to detect sounds above ~20kHz, but is audible to (for example) dogs, cats, and mice. It then listens for the first echo which matches the same pattern of 8 pulses, and reports the time elapsed.

Can we do better?

Dedicated components have the advantage of directional speakers / microphones, and very high sample rates, but they have minuscule processing power and only report a single reading: the time to first echo. Echoes contain much more information: they can reflect back from different directions, and at each boundary some of the soundwave’s energy will continue into the object rather than reflecting back. Objects will also absorb and reflect different frequencies in different proportions due to resonance. These behaviours mean we should be able to see multiple reflected echoes, and can say something about each object’s properties based on how strong the echo is at different frequencies. Doppler shift means we can even identify the relative speed of the detected objects.

Signal

The first thing to notice is that if we wish to pick up more than one echo, we’ll need to use a more complicated signal than a fixed-tone beep. The sounds are likely to be so close together that they overlap, which will be impossible to distinguish if they are all the same frequency:

Spectrogram: simulation of a 10ms 40kHz tone being reflected off surfaces 40cm and 60cm away; can you spot where the echoes begin?
sample(t)=sin⁡(2πft)wheref=40000\begin{gathered} \text{sample}(t) = \sin(2 \pi f t) \\ \text{where} \quad f = 40000 \end{gathered}

Here it is useful to look to nature. Bats use a sound known as a “chirp”. This is a sound which starts at one frequency, then sweeps to another frequency over time. As the echoes come back staggered over time, their frequency separation makes them much easier to distinguish:

Spectrogram: simulation of a 10ms 40kHz–6kHz linear chirp being reflected off surfaces 40cm and 60cm away; the echoes are easily distinguished from the immediate feedback and each other.
sample(t)=sin⁡(2π(f0+f1−f02td)t)wheref0=40000,f1=6000,d=0.01\begin{gathered} \text{sample}(t) = \sin\left(2 \pi \left(f_0 + \frac{f_1 - f_0}{2} \frac{t}{d}\right) t\right) \\ \text{where} \quad f_0 = 40000, \; f_1 = 6000, \; d=0.01 \end{gathered}

Different species of bat have different sounds; the frequency range, duration, and harmonics all vary. Some species even include a fixed tone in the middle of their call (which they can use to detect relative velocity of the target object via Doppler shift). Hunter Scott’s The DSP Behind Bat Echolocation has spectrograms of a wide selection of real bat recordings, illustrating how much variety there is.

Many bat species’ calls contain a chirp which uses a hyperbolic curve in frequency space, rather than being linear:

Spectrogram: simulation of a 10ms 40kHz–6kHz hyperbolic chirp being reflected off surfaces 40cm and 60cm away.
sample(t)=sin⁡(2πf0(log⁡(ct+1)c))wherec=1d(f0f1−1),f0=40000,f1=6000,d=0.01\begin{gathered} \text{sample}(t) = \sin\left(2 \pi f_0 \left(\frac{\log(c t + 1)}{c}\right)\right) \\ \text{where} \quad c = \frac{1}{d}\left(\frac{f_0}{f_1}-1\right), \; f_0 = 40000, \; f_1 = 6000, \; d=0.01 \end{gathered}

The advantage of such a curve is that it is not distorted by Doppler shifting, which means the echo can be detected in the same way regardless of relative velocities (very useful for a creature in flight).

It also has the nice property that higher frequencies account for proportionally less of the signal. This means that if our signal is low-pass filtered at any point (either intentionally by the speaker or microphone, or as the result of hardware limitations), only a small part of it will be lost and we should still be able to detect it reliably (though at a cost to the temporal accuracy).

We don’t have to be constrained to nature, or even to lines:

Spectrogram: simulation of a 10ms sonified image (made using the Griffin-Lim algorithm) being reflected off surfaces 40cm and 60cm away. The source image was less creepy.

So if we want to have some fun we could try locating things by bouncing smiley faces off them.

For our purposes, the properties of the hyperbolic chirp make it a good choice for getting started, even if we later revisit it to use something more fine-tuned to the situation. To avoid sharp “clicks” at the start and end of the chirp, we can use a simple linear “fade in” and “fade out”; this article uses a 1-millisecond fade in, and a 2-millisecond fade out.

Detection

Having an easily detected signal is all well and good, but how do we actually detect it?

Signal components.

One approach would be to calculate a spectrogram from the microphone (like the illustrations in the previous section) then use image analysis to locate the shape of our signal. There are well known and efficient ways to find a known shape in an image (e.g. Haar cascade classifiers), and once we find it we can easily analyse the image to see which frequencies have been absorbed (allowing us to do some rudamentary object classification). We can also subtract each signal after we find it and search again to find the next echo, even if they were overlapping. However, there is also a big disadvantage to this approach: performing frequency analysis involves a trade-off of frequency precision vs. time precision (due to the equivalent of Heisenberg’s uncertainty principle)(note 4). This significantly reduces our distance accuracy. Ideally, we would like to be able to identify echoes with sample-level precision.

The solution to this problem is simple and elegant, but needs some introduction.

Convolution

Suppose we wanted to go the other way: given a clean timeline showing the instants and intensities of the echoes, as well as our signal, we want to generate the audio which we would expect to hear. This process is called “convolution”, denoted A∗BA * B, where AA and BB are arrays of values.

To get a sense for the operation, a convolution can be implemented naïvely in Python as:

def convolve(A, B):
  result = [0 for i in range(len(A) + len(B) - 1)]
  for i in range(len(A))
    for j in range(len(B))
      result[i+j] += A[i] * B[j] # note: this * is regular multiplication!
  return result

If AA and BB have length nn, this has a computational complexity of 𝒪(n2)\text{𝒪}(n^2), which means it can slow down significantly as the length of the inputs increases. It is much more computationally efficient to use the mathematical identity:

ℱ(A∗B)≡ℱ(A)⊗ℱ(B)ℱ(A * B) \equiv ℱ(A) \otimes ℱ(B)

Where ℱℱ is the fourier transform and ⊗\otimes denotes element-wise multiplication(note 5). This identity lets us implement a convolution as:

def convolve(A, B):
  fA = fft(A)
  fB = fft(B)
  for i in range(len(fA))
    fAB[i] = fA[i] * fB[i]
  return ifft(fAB)

Where fft is some implementation of the “Fast Fourier Transform”, and ifft is the inverse; libraries with these functions exist for most programming languages.

Our AA and BB inputs must now have the same size, and depending on the Fast Fourier Transform implementation it may need to be a power-of-two, but there are techniques to avoid these limitations (discussed in part 3). For now we can pretend that it’s as flexible as the naïve approach.

The computational complexity of the Fast Fourier Transform is 𝒪(nlog⁡(n))\text{𝒪}(n \log(n)), so this is more efficient, but the part which is interesting to us is is not the speed; it is that we have turned a summation over many elements into a simple element-wise multiplication, and element-wise multiplication is reversible.

Deconvolution

Rearranging the above equation, we can see that:

ℱ(A)≡ℱ(A∗B)⊘ℱ(B)ℱ(A) \equiv ℱ(A * B) \oslash ℱ(B)

Where ⊘\oslash denotes element-wise division. This lets us “undo” the convolution which has been applied! The operation is known as “deconvolution”.

If BB is our chirp, and A∗BA * B is the audio received by the microphone, then we can use this to find the clean timeline of echoes which we are looking for (AA).

Since we need to perform this operation a lot, and the chirp (BB) does not change, we can also pre-compute Binv=1⊘ℱ(B)B_{\text{inv}} = 1 \oslash ℱ(B)} which lets us use multiplication rather than division in our loop for a small speed boost. The result of the fourier transform is actually an array of complex numbers, so calculating this inverse for each element is not quite as trivial as it could be, but it is still very easy:

1b≡b∗|b|2\frac{1}{b} \equiv \frac{b^{*}}{\left|b\right|^2}

Where:

Written as code:

def complex_reciprocal(z):
  m = 1 / (z.real ** 2 + z.imag ** 2)
  return complex(z.real * m, -z.imag * m)

(of course in real Python code you could simply write 1/z, but this is to illustrate the process using basic operations).

Dividing by 0

Division has a problem: If bb is 0, the division will explode. More generally: the accuracy will degrade as bb approaches 0. In practical terms, this means that noise at frequencies which are not used by our chirp will cause large errors in our output.

We can work around this by adjusting the equation so that it can never result in a division by zero:

1b≈b∗|b|2+c\frac{1}{b} \approx \frac{b^{*}}{\left|b\right|^2 + c}

This may look like a hack, but it actually has theoretical backing. In the optimal case, cc is a function of frequency: c(f)=1SNR(f)c(f) = \frac{1}{\text{SNR}(f)}, where SNR(f)\text{SNR}(f) is the signal-to-noise ratio calculated for a given frequency. To calculate these values we can look at the audio received by the microphone when there is no echo expected (e.g. just before we send out our signal), fourier transform it, and compare each frequency bucket with the fourier transform of our signal (both squared because we want to know the ratio of energy rather than amplitude):

SNR(f)=v×|ℱ(B)[f]|2|ℱ(inputpre-emission)[f]|2\text{SNR}(f) = v \times \frac {\left|ℱ(B)[f]\right|^2} {\left|ℱ(\text{input}_{\text{pre-emission}})[f]\right|^2}

There’s still a free constant vv (related to volume, since we don’t know how loud our signal will be when it is emitted), but this no-longer depends on the frequency so we can tune it quite easily with trial-and-error.

Combining these equations, we can simplify some terms:

1b≈b∗|b|2|b|4+(|ℱ(inputpre-emission)[f]|2v)\frac{1}{b} \approx \frac {b^{*} \left|b\right|^2} {\left|b\right|^4 + \left(\frac{\left|ℱ(\text{input}_{\text{pre-emission}})[f]\right|^2}{v}\right)}

Or if we’re feeling lazy, we can skip the SNR calculation entirely and just approximate cc as a constant: values around 1000 work quite well (this corresponds to a SNR of -30Db, reflecting the fact that many frequencies do not occur at all in our signal).

Cross-correlation

Some sources use “cross-correlation” rather than deconvolution to detect echoes. The two approaches are very closely related: the process is the same, but cross-correlation replaces 1b\frac{1}{b} with b∗b^* in the inversion above; dropping the denominator |b|2+c\left|b\right|^2 + c entirely. This means it can be turned into deconvolution by applying a collection of specifically-tuned band-passes to the input first.

In practice, using cross-correlation on its own (without these band-passes) is workable, but gives less focused results, making it more difficult to distinguish echoes from the noise and interference patterns. Deconvolution has the same performance (if we precompute BinvB_\text{inv} as noted earlier) and gives better results even when using approximated noise levels.

Simulating detection

If we can’t detect echoes in a perfect simulated environment, we don’t stand a chance in the real world, so let’s see how this deconvolution technique performs with our simulated environment from above:

Detected echoes from a simulated environment using a hyperbolic chirp with surfaces at 40cm and 60cm.

It identifies the echoes with good precision. From this data we can find the echoes programmatically with a simple heuristic:

For each sample (ordered by time):
- If the sample is larger than any others within a 0.2ms window:
  - If we do not yet have a threshold:
    - we have found the immediate feedback peak;
    - set the threshold to the current sample value * 0.01;
    - continue to the next sample.
  - Otherwise, if the sample is greater than the threshold:
    - record an echo at time = current_sample_index / sample_rate.

This finds both echoes, and positions them at 2.333ms and 3.500ms(note 6). These are very close to the true values of 2.337ms and 3.502ms (within 0.2%).

As a comparison, let’s see what the same technique finds if we use a flat tone instead of a chirp:

Detected echoes from a simulated environment using a flat tone with surfaces at 40cm and 60cm.

The echoes are still there, but it is struggling to localise them. This represents a fairly extreme example: our signal is quite long and has not finished being sent by the time these echoes return (meaning we have a lot of interference), but it demonstrates how much our choice of signal impacts the system. Let’s stick with a hyperbolic chirp for now!

Sub-sample precision

We can get even more precise distance estimates by fitting a quadratic polynomial using the samples to the left and right of each peak (using Lagrange polynomials) and differentiating it to find the maximum point of the curve:

Starting with the Lagrange polynomial definition:

ℓj(x)=∏0≤m≤km≠j(x−xmxj−xm)\text{ℓ}_j(x) = \prod_ {\begin{gathered}0 \le m \le k \\ m \neq j\end{gathered}} \left(\frac{x-x_m}{x_j-x_m}\right)

We can calculate the polynomials for a 3-point quadratic with points at xx = (-1, 0, 1):

ℓ0(x)=x−x1x0−x1x−x2x0−x2=12(x2−x)ℓ1(x)=x−x0x1−x0x−x2x1−x2=(1−x2)ℓ2(x)=x−x0x2−x0x−x1x2−x1=12(x2+x)\begin{align*} \text{ℓ}_0(x) & = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} & & = \frac{1}{2} (x^2-x) \\ \text{ℓ}_1(x) & = \frac{x-x_0}{x_1-x_0} \frac{x-x_2}{x_1-x_2} & & = (1-x^2) \\ \text{ℓ}_2(x) & = \frac{x-x_0}{x_2-x_0} \frac{x-x_1}{x_2-x_1} & & = \frac{1}{2} (x^2+x) \end{align*}

Multiplying each of these by the value at that point (v−1v_{-1}, v0v_0, and v1v_1), and collecting terms:

L(x)=(v−1+v12−v0)x2+(v1−v−12)x+v0L(x) = \left(\frac{v_{-1} + v_1}{2} - v_0\right) x^2 + \left(\frac{v_1-v_{-1}}{2}\right) x + v_0

Differentiating with respect to xx:

dL(x)dx=(v−1+v1−2v0)x+v1−v−12\frac{d L(x)}{d x} = \left(v_{-1} + v_1 - 2 v_0\right) x + \frac{v_1-v_{-1}}{2}

Equating this to 0 to find the turning point (maximum):

0=(v−1+v1−2v0)x+v1−v−12x=v1−v−12(2v0−v−1−v1)\begin{align*} 0 & = \left(v_{-1} + v_1 - 2 v_0\right) x + \frac{v_1-v_{-1}}{2} \\ x & = \frac{v_1-v_{-1}}{2(2 v_0-v_{-1}-v_1)} \end{align*}

Plugging in our measured values for vv tells us the amount we need to shift the sample by for sub-sample precision (the resulting xx will be between -0.5 and +0.5).

For this experiment, the sub-sample adjustment moves our echo readings to 2.336ms and 3.501ms(note 7); both correct to within 0.05%. In the real world we expect noise to make our accuracy much worse, but achieving this high level of accuracy in an ideal simulated environment indicates that the approach is viable.

Inferring reflection intensity

You may have noticed in the previous section that the second echo has a lower amplitude than the first, even though both represent 100% of the signal being reflected. Specifically the detected intensities are 0.028 and 0.021(note 8). This presents a problem if we want to know how much sound an object actually reflected (e.g. to check if its surface is rough or smooth), so we need to model the phenomenon and adjust for it.

As a sound wave expands, the total energy of the wave stays the same(note 9), which means the energy of the wave in any given direction decreases in proportion to the total area covered by the wavefront. For a simple scenario where the sound wave is unconstrained and can travel in all directions, this tells us that E∝d−2E \propto d^{-2}, where EE = energy in a given direction, and dd = distance travelled(note 10).

We don’t directly have the energy (EE) or the distance (dd), but for sound waves E∝A2E \propto A^2 (where AA = amplitude); and in any given environment the speed of sound in air stays relatively constant: d≈vtd \approx v t, where tt is the elapsed time, and vv is some unknown but constant velocity(note 11). Put another way: d∝td \propto t.

Combining all of these relations, we can say:

reflection intensity∝received amplitude×time since signal was emitted\text{reflection intensity} \propto \text{received amplitude} \times \text{time since signal was emitted}

The constant of proportionality does not matter to us, since the user may have an unknown gain level (volume) on their speakers and microphone anyway, so for our purposes we’ll just say the constant is 1(note 12).

Applying this per-sample scaling to our simulation (after performing the deconvolution), we get spikes which are of a consistent size:

Detected echoes from a simulated environment using a hyperbolic chirp with surfaces at 40cm and 60cm, correcting for attenuation.

The detected intensities are now 0.227 and 0.246; within 10% of each other (this remaining difference is probably due to inaccuracies in the deconvolution, and the fact that a quadratic curve is not a perfect fit for the data).

Note that for this to give accurate results, we must ensure that the measured time is “like for like”: if our deconvolution is giving us the midpoint of detected samples, then we must measure time from the midpoint of the emitted sample, not the moment we started playing it.

Latency

Outside our simulation, we have a problem: we do not know exactly when the chirp was emitted from the speakers, nor do we know exactly when the microphone picked up its echo. All we know is when we told the speakers to produce a sound, and when the microphone reported the sound to us. There is an unknown amount of latency on both sides.

There is an unknown amount of latency from processor to speaker, and from microphone back to processor

For fixed hardware, this is not a concern. We can simply measure the latencies in a known environment and hard-code them into our software. But if we want our program to work with arbitrary hardware and software stacks, then we cannot rely on knowing these latencies in advance.

The WebAudio API does include some latency reporting: if the browser supports the feature, it can tell us the predicted output and input latencies of a device. This may include software latencies (due to signal processing and queues) as well as hardware / firmware latencies (such as the time taken to send a signal to a bluetooth speaker). But these values are not guaranteed to account for everything, and they may not be available at all (especially if the user has external speakers with unknown characteristics). These are a good starting point when they are available, but they are not enough to solve the problem.

Error tolerance

So amidst all this uncertainty, how certain do we need to be?

Earlier in this article, we established that there is a fundamental error in echolocation distances of ±5% due to variations in the speed of sound. If we can’t measure the air temperature, pressure, or humidity, and we don’t have any points of comparison to calibrate against, then it doesn’t matter how good our processing is: we cannot reduce this error. This gives us something to compare our latency-induced errors with.

Also, for our purposes, we don’t actually care what the speaker or microphone latencies are individually; we only care about the sum of these latencies. This reduces our problem to 2 unknowns (total latency, and distance to object).

Total measurement error for different latency estimates.

If our latency estimate is out by 0.6ms, then we will have doubled the overall systemic error at 2 metres. At 1 metre, the threshold for doubling the error is just 0.3ms (equivalent to 14 samples at 48kHz).

We may also want to consider the impact on our reflection intensity normalisation. We can compare our example echoes at 40 and 60cm:

Reflection intensity error for objects at 40cm and 60cm with different latency estimates.

In this case, underestimating the latency is significantly better than overestimating it; an underestimation of 2ms results in an 18% error in the relative intensities, but an overestimate of 2ms introduces a 65% error.

The next article will discuss how we can use measurements from multiple components to infer latencies automatically (if we have enough components available).

Free parameters

We’re left with a few free parameters which we can experiment with:

Next steps

We have a way to get 1-dimensional data about multiple object positions by using echolocation, and we can classify these objects somewhat based on their reflectivity. Part 2 (not yet posted) takes these results a step further: inferring 2-dimensional positions of objects when multiple speakers and/or microphones are available. Part 3 (not yet posted) shows how to build these techniques into a browser-based application.

Further reading

Footnotes

  1. The speed of sound in air varies with temperature, humidity, and pressure; in typical atmospheric conditions it can be anywhere between 330m/s (cold air) and 360m/s (hot air). ↩︎

  2. There are some more uncertainties in this number which will be discussed in later sections. ↩︎

  3. Using interpolation between samples can give us a little more granularity, and an approach for this is discussed later, but it means our detection will become more susceptible to noise, and we will still not be able to discern multiple echoes which are very close together. ↩︎

  4. Mathematically, this happens because we are throwing away the phase information of the frequencies. ↩︎

  5. Conceptually, this transformation of convolution into elementwise multiplication is similar to the way that taking the logarithm turns multiplication into addition, i.e. ln⁡(A×B)=ln⁡(A)+ln⁡(B)\ln(A \times B) = \ln(A) + \ln(B). ↩︎

  6. These numbers and the graphs shown above all correspond to a sample rate of 96kHz, which is higher than most computers use, but we get the same outcome with 48kHz and a hyperbolic chirp from 20kHz to 6kHz (since we are limited to a maximum of 24kHz at this sample rate). At a sample rate of 44.1kHz the results only differ by 0.01ms. ↩︎

  7. Again we get the same values at 48kHz, and very similar values at 44.1kHz. The sample rate does affect the measurement, but only in the 3rd decimal place. We’re not seeing any significant impact from sample rate here because the simulated environment has no background noise. ↩︎

  8. To avoid aliasing artifacts, these values come from the same polynomial interpolation as the sub-sample timings. The equation is:

    intensitypeak=12(x2−x)v−1+12(x2+x)v1+(1−x2)v0=v0+(v1−v−1)28(2v0−v−1−v1)\begin{align*} \text{intensity}_\text{peak} & = \frac{1}{2}(x^2-x)v_{-1}+\frac{1}{2}(x^2+x)v_{1}+(1-x^2)v_0 \\ & = v_0 + \frac{(v_1-v_{-1})^2}{8(2 v_0-v_{-1}-v_1)} \end{align*}
    ↩︎
  9. Small amounts of energy are lost as heat, motion, and diffraction, but we can ignore those effects for this calculation. If we want to account for energy loss as heat, we can model it as Etotal∝e−αdE_\text{total} \propto e^{-\alpha d}, where α\alpha is proportional to the sound frequency squared and can be calculated using Stokes’ law. ↩︎

  10. Since we only care about proportionality rather than actual values, this relationship is valid even if there are some boundaries, such as the floor or a wall. In practice, it only starts to fail if we are in a very confined environment, such as a narrow corridor. ↩︎

  11. Probably around 343m/s although we do not care what its actual value is here. ↩︎

  12. If needed, we can infer the overall scaling factor (which contains these constants as well as the gain settings) from the signal we receive by looking at the intensity of the immediate feedback peak, as we did for detecting echoes earlier. ↩︎