Autonomous phase-sensitive radio-echo sounder theory

This page describes some of the theory behind the Autonomous phase-sensitive radio-echo (ApRES) sounder, including

• a description of the linear chirps the system emits,

• how individual and multiple reflectors are represented in the returned signal,

• how the range to these reflectors is encoded in the frequency content of the returned signal, and

• how to extract the range to reflectors using a fourier transform.

Overview

Many of the material below can be found in the ApRES manual, [Brennan et al., 2014], and [Nicholls et al., 2015].

ApRES emits ‘chirps’, which consist of continuous radio waves lasting one second. During each chirp the frequency of the emitted radio wave increases linearly with time from \(f1\) to \(f2\). The bandwidth of the system is defined as \(B = f2-f1\). The signal is transmitted downwards into the ice sheet through one antenna (referred to by ‘Tx’ for transmit) and is partly reflected back to the radar’s receiving (or ‘Rx’) antenna, where it is compared to the transmitted signal to determine the range to sub-surface reflectors.

To be more specific, throughout each chirp, the radar continually compares the received signal to the transmitted signal. Because the transmitted signal is increasing in frequency throughout the chirp and because the received signal was transmitted a short time earlier than it is received, the received signal is always lower frequency than the transmitted signal. In fact, the difference in the frequencies of the two signals is linearly proportional to the time it took the signal to travel to the reflector and back to the radar. Finally, if we know the speed of radio waves in the ice, we can use this travel time to compute the range to the reflector.

We will demonstrate this for the simple case of one reflector, before explaining how it applies when there are multiple reflectors.

import numpy as np
import matplotlib.pyplot as plt
import scipy
from numpy.random import default_rng
import numpy as np
from scipy.signal import argrelextrema

1. Define a reflector

Suppose that there is just one reflector beneath you when you deploy ApRES. In reality there are hundreds, but we will assume there is just one.

Let’s define the range to the reflectors as

R = 110  # units [m]

2. Compute the two-way travel time to the reflector

If we assume uniform dielectric properties in the ice, we can consider how long it would take for a transmitted signal to travel to a reflector which has a range \(R\), and back to the receiving antenna. This is called the ‘two-way travel time’ to a reflector and it is given by $\( \tau = 2R\frac{\sqrt\epsilon}{c}, \)$

where \(c\) is the speed of light in a vacuum and \(\epsilon\) is the relative permittivity of ice.

We’ll use constants from [Brennan et al., 2014], to compute \(\tau\) for our reflector:

c = 299792458
ep = 3.1
def tau(R):
    return 2*R*np.sqrt(ep)/c
print(f'two-way travel time, τ = {tau(R):.3e} s.')

two-way travel time, τ = 1.292e-06 s.

3. Plot the frequency of the transmitted signal

ApRES transmits one-second-long radio-wave chirps with a frequency that increases linearly with time from \(f_1 = 200\) MHz to \(f_2 = 300\) MHz:

T=1                 # chirp duration
f_1 = 200e6         # starting frequency
f_2 = 400e6         # ending frequency
f_c = (f_1+f_2)/2   # center frequency
B = f_2 - f_1       # bandwidth

The bandwidth \(B\) and the center frequency \(f_c\) are set by the design of the radar.

Let’s define a simple function that returns the frequency of the signal as a (linear) function of time. So that later we can demonstrate what it looks like when the returned signal is delayed by a given time, we add to the function an optional delay in seconds. The sample frequency is set at 40 kHz by default.

sampling_frequency = 40000    # [Hz]
t = np.linspace(0, T, sampling_frequency*T)   # time vector
def chirp(delay = 0):
    return f_1 + (t-delay)/T*(f_2-f_1)

Next we use this function to create and plot the frequency of the transmitted signal:

tx = chirp()

fig, ax = plt.subplots()
ax.plot(t,tx/1e6,'-', label ='transmitted signal')
ax.set_xlabel('time, $t$ [s]')
ax.set_ylabel('frequency, $f$ [MHz]')
ax.set_title('Frequency of transmitted signal during one chirp')
ax.set_xlim(0, T);
plt.show()

4. Plot the frequency of the received signal

The signal received by ApRES is the same as the transmitted signal except that it has travelled to the reflector and back. This means that the received signal is delayed by the two-way travel time \(\tau\) to the reflector.

We can use chirp, defined above, to generate a plot representing the frequency of the received signal and plot it with the frequency of the transmitted signal.

rx = chirp(tau(R))

fig, ax = plt.subplots()
ax.plot(t,tx/1e6, '-C0', label ='transmitted signal')
ax.plot(t,rx/1e6, '-C1', label ='received signal')
ax.set_xlabel('time, $t$ [s]')
ax.set_ylabel('frequency, $f$ [MHz]')
ax.set_title('Frequency of transmitted and received signals during one chirp')
ax.legend()
ax.set_xlim(0, 5*tau(R));
ax.set_ylim(f_1/1e6, f_1/1e6+5*tau(R)*B/T/1e6);

Note that we have had to zoom into a the first few microseconds of the chirp to be able to see the delay in the received signal because the two-way travel time, \(\tau\), is so small compared to the total duration of the chirp.

We can see by eye that the horizontal distance between the two signals is the two-way travel time, \(\tau\), as required:

tau(R)

1.292060425871348e-06

5. Compute the frequency difference

From the plot above you can see that the delay in the received signal causes a difference in the frequency of the transmitted and the received signals at each moment throughout a chirp. This frequency difference is equal to the vertical distance between the two curves.

Given the simple geometry of the plot, it is straightforward to relate this frequency difference, \(f_d\), to the two-way-travel time, \(\tau\):

\[ f_d = K\tau, \]

where \(K\) is slope of the curves in the plot above, i.e., the rate of frequency increase in Hz/s, given by \(K=B/T\), where \(T\) is the chirp duration (1 s) (Equation 1 from Brennan et al.).

This can be understood by looking at the plot above and recognizing that \(K\) is the slope of the curves, the vertical distance between them is the frequency difference, \(f_d\), and the horizontal distance between them is the time delay, \(\tau\).

Let’s compute the frequency difference between the transmitted signal and the signal received from our reflector at a range of 110 m:

K = B/T             # [Hz/s]
f_d = K*tau(R)         # [Hz]

print(f'The frequency difference for a reflector {R} m deep in ice is {f_d:.1f} Hz.')

The frequency difference for a reflector 110 m deep in ice is 258.4 Hz.

Note that this frequency is much lower than the transmitted and received signals. It is in the range that is audible to humans, rather than the MHz range of the transmitted and received signals. In fact, the frequency is very close to middle C.

6. How does ApRES measure this frequency difference?

Above we defined the range to a reflector, \(R\), and computed the difference, \(f_d\), between the frequencies of the transmitted and received signals that this reflector would yield. In reality, we don’t know \(R\) apriori, we want to use ApRES to estimate \(f_d\) so we can compute \(R\). To estimate \(f_d\), ApRES mixes the received and transmitted signals. The two signals interfere with each other to generate a new signal consisting of the original two signals plus a component with a frequency equal to the difference frequency \(f_d\). The original two signals (and a higher frequency component that also arises from the mixing) are filtered out, leaving what is called the ‘deramped’ signal. This deramped signal is what is saved to disk by ApRES. We refer to these second-long deramped signals as ‘chirps’.

To compute the range to our single reflector we simply need to compute the frequency of the deramped signal (next section) and plug this into an expression that relates \(f_d\) to \(R\). We can derive such an expression by combining the two equations above, \(f_d = K\tau\) and \(\tau = 2R\frac{\sqrt\epsilon}{c}\), to give

\[ R = \frac{c}{2\sqrt\epsilon} \frac{f_d}{K}. \]

7. Compute the frequency of this deramped signal with a fourier transform

In our simple case, the deramped signal looks like a sine wave:

\[ s = \sin{2\pi f_d t} \]

s = np.sin(2*np.pi*f_d*t) 

Let’s plot this deramped signal:

f, ax = plt.subplots(figsize=(18,5))
ax.plot(t, s)
ax.set_title(f'deramped signal (what would be recorded by ApRES if there was one reflector at {R} m)')
ax.set_xlabel('t [s]')
ax.set_xlim(0, 0.05)

(0.0, 0.05)

As mentioned above, we need to compute the frequency of this signal. We do this using a fourier transform. This video is an incredibly clear explanation of how fourier transforms achieve this.

The cell below computes the discrete fourier transform of s and the frequency bins. By default np.fft.fft produces as many frequency bins as there are time-domain samples (=len(s)), evenly distributed between 0 Hz and the sampling frequency. Therefore the values of the frequencies are computed by multiplying the indexes by sampling_frequency/no_of_samples. Note also that the frequency spectrum is normalized by the number of samples.

def fft(s):
    no_of_samples = len(s)
    S = np.fft.fft(s)/no_of_samples         
    indexes      = np.arange(no_of_samples) 
    frequencies  = indexes * sampling_frequency/no_of_samples
    return S, frequencies
S, frequencies = fft(s)

Next we compute the absolute value of the complex numbers in S and plot them against the frequencies yielding the usual depiction of frequency domain.

fig,ax = plt.subplots(figsize=(18,3))
ax.set_title('frequency domain')
ax.plot(frequencies, np.abs(S))
ax.set_xlabel('frequency [Hz]')
ax.set_ylabel('amplitude')
ax.set_xlim(0, 2*f_d);

Notice that there is a peak at the frequency of the deramped signal, \(f_d = 258.4\) Hz

8. Compute the range to the reflector

We can now use the equation we derived above (\(R = \frac{c}{2\sqrt\epsilon} \frac{f_d}{K})\) to convert these frequencies to range,

def range(frequencies):
    return frequencies * c/(2*K*np.sqrt(ep))
r = range(frequencies)

and plot the result

fig,ax = plt.subplots(figsize=(18,3))
ax.set_title(f'amplitude-range plot showing our reflector at {R} m')
ax.plot(r, np.abs(S))
ax.set_xlabel('range [m]')
ax.set_ylabel('amplitude')
ax.set_xlim(0, 2*R);

This range-amplitude plot is the ApRES equivalent of the usual time domain plot you would get from an impulse radar system: i.e. one which sends out a single pulse of radio-wave energy and records the echo from the reflector(s) beneath.

As expected, in the plot above we see a peak at 110 m. We can detect the position of the peak using argrelextrema, a function from the package scipy.

peaks = argrelextrema(np.abs(S), np.greater)      # this function finds local maxima  
print(f'reflector detected at {r[peaks[0][0]]:.3f} m.')

reflector detected at 109.825 m.

This is very close to the prescribed range of 110 m. However, because frequencies is quantized and (in our case) increases in 1 Hz increments, the retrieved peaks are restricted to be an integer number of hertz. This is why Brennan et al. call this a coarse range measurment.

9. The effect of bandwidth \(B\) and chirp duration \(T\)

The peak in the plot above has a non-zero width that is related to length of the deramped signal. For example, increasing the bandwidth, \(B\), would decrease the width of the peak.

The function below collects together the important bits of the code above and generates a range-amplitude plot, and allows us to change key parameters like the bandwidth and chirp duration to explore their impact.

For example, using a much smaller bandwidth of only 30 MHz results in a much wider peak in the range-amplitude plot.

T = 1  # [s]
f_1 = 200e6  # [Hz]
f_2 = 230e6  # [Hz]
sampling_frequency = 40000   # [Hz]
R = 80   # [m]
t = np.linspace(0,T,int(sampling_frequency*T))   # time vector
B = f_2 - f_1       # bandwidth
K = B/T             # [Hz/s]

def amp_range_plot():

    tx = chirp()
    f_d = K*tau(R) 
    s = np.sin(2*np.pi*f_d*t) 
    S, frequencies = fft(s)
    r = range(frequencies)
    
    f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18,5))
    ax1.plot(t, tx/1e6,'-', label ='transmitted signal')
    ax1.set_xlabel('time, $t$ [s]')
    ax1.set_ylabel('frequency, $f$ [MHz]')
    ax1.set_title('Frequency of transmitted signal during one chirp')
    ax1.set_xlim(0, T)

    ax2.plot(t, s)
    ax2.set_title('deramped signal')
    ax2.set_xlabel('t [s]')
    ax2.set_xlim(0, 0.05)

    ax3.set_title('amplitude-range plot showing our two reflectors')
    ax3.plot(r, np.abs(S))
    ax3.set_xlabel('range [m]')
    ax3.set_ylabel('amplitude')
    ax3.set_xlim(0, 2*R)
    peaks = argrelextrema(np.abs(S), np.greater)     

    print(f"starting frequency: {f_1} Hz")
    print(f"ending frequency: {f_2} Hz")
    print(f"B: {B} Hz")
    print(f"K: {K} Hz/s")
    print(f"T: {T} s")
    print(f"R: {R} m")
    print(f"sampling_frequency: {sampling_frequency} Hz")
    print(f"two-way travel time: {tau(R):.3e} s")
    print(f"frequency of deramped signal: {f_d:.1f} Hz")
    print(f"number of cycles in deramped signal: {f_d*T:.1f}")
    print(f'reflector detected at {r[peaks[0][0]]:.3f} m.') 
    
    plt.show()
amp_range_plot()    

starting frequency: 200000000.0 Hz
ending frequency: 230000000.0 Hz
B: 30000000.0 Hz
K: 30000000.0 Hz/s
T: 1 s
R: 80 m
sampling_frequency: 40000 Hz
two-way travel time: 9.397e-07 s
frequency of deramped signal: 28.2 Hz
number of cycles in deramped signal: 28.2
reflector detected at 79.460 m.

The lower resolution we see when \(B\) is reduced is a result of their being fewer cycles in the deramped signal than in cases when \(B\) is larger. The number of cycles in each chirp is given by $\( N = f_d T. \)$

The more cycles in a signal, the more precisely the fourier transform can determine its frequency. This is a fundamental property of fourier transforms and its crops up in many applications, including in quantum mechanics as the Heisenberg uncertainty principle. This is a great video on the subject from the same source as the video linked above about fourier transforms: https://www.youtube.com/watch?v=MBnnXbOM5S4. This video graphically demonstrates the idea that it is difficult for a fourier transform to determine the frequency of a short signal containing few cycles.

Interestingly, you will notice that if you play around with the numbers in the cell above (all else being equal) the width of the peak does not depend on the chirp duration, \(T\). At first sight this is surprising because you would expect that a longer duration chirp would contain more cycles and would therefore result in a narrower peak. However, notice that as you vary \(T\), the number of cycles, which is printed out beneath the cell, does not change. This is because increasing \(T\) while keeping \(B\) constant decreases the rate of change of the frequency, \(K\). This means that by the time the signal arrives back at the radar, the transmitted signal has not increased as much as it would have done if \(K\) were larger. Therefore the frequency of the deramped signal, \(f_d\), is lower than it would have been.

A longer chirp duration tends to increase the number of cycles per chirp, but the decrease in frequency of the deramped signal \(f_d\) counteracts this. In fact, the two effects balance exactly, resulting in the number of cycles per chirp being independent of \(T\).

Mathematically, we can show this by substituting \(f_d = K\tau\) into the equation above,

\[ N = K\tau T. \]

then substituting in \(K=B/T\), which describes how \(K\) decreases if \(T\) increases (as described in words above) gives

\[ N = B\tau. \]

From this we can see how the bandwidth \(B\) effects the peak width, while \(T\) does not; \(T\) cancelling out in the last step above corresponds to the two effects descrbed above balancing each other exactly.

10. Add more reflectors

So far we have considered just one reflector. In general we simultaneously receive signals from reflectors at a whole range of depths. To give a feel for what this looks like, we will plot the frequency of ten signals with their delays selected randomly, to represent signals from ten reflectors at different depths.

First let’s redefine a few functions to make sure to avoid any issues caused by varying parameters in the previous section.

c = 299792458
ep = 3.1
def tau(R):
    return 2*R*np.sqrt(ep)/c
T=1                 # chirp duration
f_1 = 200e6         # starting frequency
f_2 = 400e6         # ending frequency
f_c = (f_1+f_2)/2   # center frequency
B = f_2 - f_1       # bandwidth
K = B/T             # [Hz/s]

sampling_frequency = 40000    # [Hz]
t = np.linspace(0,T,sampling_frequency*T)   # time vector
def chirp(delay = 0):
    return f_1 + (t-delay)/T*(f_2-f_1)
def fft(s):
    no_of_samples = len(s)
    S = np.fft.fft(s)/no_of_samples         
    indexes      = np.arange(no_of_samples) 
    frequencies  = indexes * sampling_frequency/no_of_samples
    return S, frequencies
def range(frequencies):
    return frequencies * c/(2*K*np.sqrt(ep))

Next, let’s define our ten reflectors at random depths and plot the frequency of the received signals from each reflector.

R = 110
rx = chirp(tau(R))

fig, ax = plt.subplots()
ax.plot(t, tx/1e6, '-C0', label ='transmitted signal')
ax.plot(t, rx/1e6, '-C1', label ='signal received from 110 m', linewidth = 4)
# generate some randomly delayed signals
rng = default_rng(seed = 4321)
ranges = np.absolute(rng.uniform(low = 0, high = R*10, size = 10)) # np.array([110, 3, 30]) #

print(f"Plotting the frequencies of signals reflected from the following depths:")
for i, r in enumerate(ranges):
    print(f"{r:.3}")
    rx_rnd = chirp(tau(r))
    ax.plot(t,rx_rnd/1e6,'-', label =f'signal received from {r:.1f} m', linewidth = 1)

ax.set_xlabel('time, $t$ [s]')
ax.set_ylabel('frequency, $f$ [MHz]')
ax.set_title('Frequency of transmitted and received signal during one chirp')
ax.legend()
ax.set_xlim(0, 5*tau(R));
ax.set_ylim(f_1/1e6, f_1/1e6+5*tau(R)*B/T/1e6);

Plotting the frequencies of signals reflected from the following depths:
3.38
8.89e+02
5.86e+02
1.36e+02
1.01e+02
9.62e+02
58.4
9.64e+02
5.56e+02
2.89e+02

In practice, signals from all the reflectors are combined together and arrive back at the radar simultaneously, thoughout the chirp. This means that the signal received by the radar is the sum of signals reflected by many reflectors - hundreds or thousands in reality, compared to the ten we plotted above. Moreover, the signals all have slightly different frequencies differences, \(f_d\), depending on the range to each reflector. Let’s generate the deramped signal that these ten reflectors would yield and plot it.

We first compute the frequency differences for each reflector:

freq_differences = K*tau(ranges)         # [Hz] 

For each frequency difference we then generate a signal signal with that frequency:

chirp_list = [np.sin(2*np.pi*freq_difference*t)  for freq_difference in freq_differences]

Finally, we sum them together to create one signal.

chirp_array = np.stack(chirp_list, axis=1)
s = chirp_array.sum(axis=1)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18,5))

ax1.plot(t, s, label ='chirp from multiple reflectors')
ax1.set_xlabel('time, $t$ [s]')
ax1.set_ylabel('chirp amplitude [V]')
ax1.set_title('the full chirp')
ax1.set_xlim(0, 1)
ax1.legend()

ax2.plot(t, s, label ='chirp from multiple reflectors')
ax2.set_xlabel('time, $t$ [s]')
ax2.set_ylabel('chirp amplitude')
ax2.set_title('the first 1/10 of a second')
ax2.set_xlim(0, 0.01)
ax2.legend()

<matplotlib.legend.Legend at 0x139453210>

The plot above shows what ApRES would record to disk in our simple 10-reflector case. The signal has ten frequency components, each correspondong to one reflector. The challenge now is to estimate these frequency componets to determine the range to all the reflectors. As above, we use a fourier transform.

Using the same approach as we did in the single-reflector case, we apply a fourier transform to the deramped signal to extract the frequency components, then convert the frequencies to range:

S, frequencies = fft(s)
r = range(frequencies)

fig,ax = plt.subplots(figsize=(18,3))
ax.set_title('amplitude-range plot showing our ten reflectors')
ax.plot(r, np.abs(S))
ax.set_xlabel('range [m]')
ax.set_ylabel('amplitude')
#ax.set_xlim(0, np.max(ranges)*1.1);

Text(0, 0.5, 'amplitude')

FInally, we can extract the peaks as we did before.

len(S)

40000

S_first_half = S[:int(len(S)//2)]
peaks = argrelextrema(np.abs(S_first_half), np.greater)     
print(f'reflectors detected at {r[peaks[0]]} m.')    

reflectors detected at [  3.40541349  58.31770596 101.31105123 135.79086279 289.03446969
 555.93375172 586.15679641 889.23859674 962.02931002 963.73201676] m.

These are all close to the ranges we prescribed:

np.sort(ranges)

array([  3.37514327,  58.38920621, 101.12516631, 135.63533655,
       288.82637573, 556.00188035, 586.31305979, 889.26122168,
       961.89128611, 963.510313  ])

np.sort(ranges)-r[peaks[0]]

array([-0.03027021,  0.07150025, -0.18588493, -0.15552623, -0.20809396,
        0.06812864,  0.15626338,  0.02262495, -0.13802391, -0.22170375])

11. Summary

• ApRES emits one-second-long chirps.

• ApRES detects the range to sub-surface reflectors by combined the frequency of a transmitted signal, which continually increases throughout each chirp, to the frequency of the signal received from the reflectors.

• The transmitted and received signal are in the MHz range.

• This combination results in a ‘deramped’ signal which is in the audio frequency range.

• The deramped signal is saved to disk by ApRES.

• The frequency components of the deramped signal are extracted using a fourier transform to determine the range to sub-surface reflectors.

On the next page we apply this theory to real ApRES data collected in Antarctica.