mu2718 - GPS Walkthrough 1

This article is step 1 of the GPS walkthrough series. Here we give an overview over the whole post series and introduce every step required to get from a raw GPS signal to a position fix. The main results are shown while detailed analysis and explanations are following in later posts.

1 Introduction

In this post series we investigate GPS L1 C/A, the legacy GPS signal with its coarse acquisition code which is the most widely used signal for navigation and sent by all GPS satellites in operation.

We will rely in the following to an external python module developed for this series, the gps_walkthrough module. Open the collapsed code sections below to get to know how it is used an how its outputs are processed. Check the module source code if you want deeper insights into the details of the processing steps before the more detailed blog articles appear. See also the references section below for a list of sources where I have the information collected from.

To compute the following steps for yourself online, open it in Google Colab. You can even feed it your own recorded GPS radio wave to play with the satellite signal and determine your location.

Now enjoy the journey!

Code

import gps_walkthrough as gpswt  # see https://github.com/mu2718/gps-walkthrough/

# the only other packages we will rely on:
import numpy as np
import matplotlib.pyplot as plt

2 Recording Radio Wave

An active GPS antenna picks up the electromagnetic (EM) wave sent by a satellite (or ‘space vehicle’ in GPS lingo) and feeds it to an integrated signal amplifier and filter. Its output signal is directly proportional to the incoming EM wave at any time. This output is wired up with a software-defined radio (SDR) receiver. This device is able to select a frequency-band and digitize it. Using this setup, we get an unprocessed, digital representation of the incoming GPS signal.

I am using the HackRF SDR, a signal generator from Siglent delivering a high-accuracy clock to the SDR, and the u-blox ANN-MB GPS antenna.

Analyze your own GPS radio wave?

For recording the GPS radio wave, have the following equipment ready:

SDR receiver (e.g. HackRF or RTL-based) capable to receive at 1575 MHz. No precise clock reference as I use here is required.
Active GPS antenna for the L1 band with at least 20dB gain. A cheap 10$ one is sufficient (e.g. MikroTik ACGPSA).
Bias tee voltage supply if the SDR does not support powering the antenna. This is not required for HackRF with antennas accepting 3.3 Volts.

Alternatively, a synthetic recording can be generated using gps-sdr-sim.

Wait for the walkthrough step 2 to get more detailed instructions…

With the following command, the HackRF records the GPS L1 signal at 1575.42 MHz with appropriate settings for amplifiers and sampling rate of 4 MHz:

$ hackrf_transfer -r wave.dat -f 1575420000 -p 1 -a 1 -l 40 -g 40 -s 4000000

The specified file contains the recorded, raw in-phase and quadrature components (IQ) of the EM wave amplitude. We can read in this wave file and visualize a snippet of it in Figure 1.

Code

sampling_rate = 4e6  # 4 MHz sampling rate of recording

# read in the recorded radio wave file (change dtype='byte' for data from `hackrf_transfer`)
baseband = gpswt.SdrWave.from_raw_file('./wave.dat', sampling_rate=sampling_rate, 
                                       #, max_samples=100000000)  # limit samples for fast experiments
                                       dtype='float32') 
# cut the first 0.1 seconds during antenna power up
baseband = baseband.get_interval(from_time=0.1)
print(f'Wave recording read. Duration: {baseband.duration():0.2f} s')

# visualize a snippet 
t0 = 0.2   # select start time
dt = 1e-3  # time window of 1ms

signal_slice = baseband.get_interval(from_time=t0, to_time=t0+dt)
t = np.linspace(t0, t0+dt, len(signal_slice.samples))

plt.figure(figsize=(10,2))
plt.plot(t, np.real(signal_slice.samples),'-', label='In-Phase $I$')
plt.plot(t, np.imag(signal_slice.samples),'-', label='Quadrature $Q$')
plt.xlabel('Receiver Time $t$ [s]')
plt.ylabel('Amplitude [a.u.]')
plt.legend()
plt.show()

Wave recording read. Duration: 84.26 s

Figure 1: Snippet of 1 ms of the recorded radio wave baseband signal with its in-phase and quadrature amplitudes.

If you want to get more detailed instructions for recording the GPS signal, understand EM waves and SDRs, IQ components and how they can be processed using NumPy, wait for the walkthrough step 2 to get posted. In the meantime, have a look into the source code of the SdrWave class.

3 Signal Acquisition

Having the radio wave recorded, we can look for a satellite signal. Those of us having experience with pre-internet radio know: Tuning-in on a channel means essentially turning one knob to set the frequency right (and sometimes redirecting and stretching the antenna). GPS L1 C/A is not much different, but we have four knobs which have to be set, the demodulation parameters:

Frequency Shift: GPS L1 satellites send on a carrier frequency of 1575.42 MHz with very high accuracy. But due the satellites high speeds of up to 4000 m/s relative to us, the Doppler effect physically shifts the received frequency up to $\pm 5000$ Hz. Despite its relatively small magnitude, the Doppler shift has to be compensated for signal reception. Note that changes in velocity of the satellite relative to the receiver (e.g. in a fly-over) results in changes of this parameter over time and thereby need to be readjusted constantly.
An inaccurate receiver clock, i.e. the clock built into the SDR receiver, is an additional source of frequency shift: As we base all our measurement on it, a too slow clock will show up as too high frequencies in the recording and vice versa. Since I use an external high-accuracy clock, the frequency shift can be fully attributed to the Doppler effect in the following.
PRN number (C/A Code): All satellites send on the same frequency. Using a code-division multiple access (CDMA) method, where each satellite sends its own code, we can distinguish EM wave contributions from different senders. These code sequences are called C/A codes and are identified by their pseudo-random noise or PRN numbers, ranging from 1 to 37. The PRN number has therefore to be set to select a specific satellite to listen to.
Code Delay: The C/A code sequence of a satellite can only be detected if we know its time of reception to within a microsecond. The code delay is defined by the time the code is sent, as expected by the receiver based on his own inaccurate clock, and the actually observed time of reception. Note that changes in distance between the receiver and a satellite modifies this value by 1 $\mu s$ per 300 meters distance. Since satellites travel at up to 4000 m/s, we need to constantly readjust the delay to maintain reception.
Carrier Phase Angle: GPS L1 satellites send telemetry data (GPS time, orbital parameters, …). It is transmitted in the EM wave using the binary phase-shift keying (BPSK) modulation technique. In order to demodulate this information, we need to distinguish between the sine and cosine contributions of the carrier frequency, i.e. the phase angle of the received wave. Note that changes in distance to a satellite or atmospheric conditions modify the phase over time.

The GPS signal acquisition procedure tries to find the value of these parameters for optimal signal reception by going through all combinations of parameters settings and picking the optimal one.

First Signal Detection

For the very first millisecond of our radio wave recording, let’s try to find the signal sent by the satellite which emits the C/A code PRN 16 by performing the acquisition:

Code

acq = gpswt.Acquisition(prn_id=16, sampling_dt=baseband.sampling_dt)  # look for PRN 16, I know it is there :) 
acq_data = acq.search(baseband.get_interval(to_time=0.001),           # look at first 1 ms = 1 C/A code period
                      delta_freq_range=(-5000, 5000),                 # increase to 100k for HackRF w/o external clock
                      delta_freq_step=10)[-1]                         # fine-grained Doppler steps 

print(f"Signal Power:  {np.abs(acq_data['correlator'])**2 * 100:.2f} % of total")
print(f"Doppler Shift: {acq_data['delta_freq']} Hz")
print(f"Code Delay:    {acq_data['delay'] * 1e6} microseconds")
print(f"Carrier Phase: {np.angle(acq_data['correlator']):.2f} radians")

Signal Power:  1.69 % of total
Doppler Shift: -2830.0 Hz
Code Delay:    13.0 microseconds
Carrier Phase: 1.04 radians

These are the optimal demodulation parameter settings the acquisition procedure has found. Some observations we can make already:

Despite that the satellite’s signal power is below 2% of the total power that the SDR records, we can isolate it and distinguish it from noise and other satellites.
A negative Doppler shift means that the received frequency is lower than it was sent by the satellite. As we know it from everyday acoustic experiences, we can infer that this satellite is heading away from us. We will confirm this later on.

The scanning on Doppler shift and code delay performed during acquisition is shown in Figure 2. We can observe that a discrepancy of only 1 $\mu s$ of in code delay and 1000 Hz in frequency shift results in total loss of the signal. This distinct signature of the peak indicates that our demodulation works and that we have detected a signal sent by a satellite!!

Code

plt.subplot()
plt.imshow(acq._scan_data * 100, 
           extent=(-5000,5000, 0, 1000),  # scan parameter space
           cmap='OrRd', origin='lower', aspect='auto')
# center on global maximum
plt.ylim([acq_data['delay']*1e6 - 5, acq_data['delay']*1e6 + 5])
plt.xlim([acq_data['delta_freq'] - 2000, acq_data['delta_freq'] + 2000])

plt.xlabel('Doppler Frequency Shift $\Delta f$ [Hz]')
plt.ylabel('Code Delay $\\tau$ [$\\mu$s]')
plt.colorbar(label='Signal Power [%]')
plt.show()

Figure 2: Received signal power of one satellite under different demodulation parameters. The full acquisition search region is zoomed in on the peak power.

Search Satellites

We can search for all satellites in the first few milliseconds of our signal using the acquisition procedure, see Figure 3. Since we expect to see only a handful of satellites, we have to conclude that a the noise floor of around 0.5% power is present in all PRNs: We get spurious signals not originating from senders but noise on all PRNs.

Five satellites can be identified that send well above the noise floor. Since at least 4 satellites are required for a position fix, we can proceed.

Code

prn_ids = range(1, 38)  # PRN IDs to search for

signal_strength = np.zeros(len(prn_ids))
for i, prn_id in enumerate(prn_ids):
    signal_slice = baseband.get_interval(to_time=0.05)  # take first 50 ms

    acq = gpswt.Acquisition(prn_id, sampling_dt=baseband.sampling_dt)
    acq_data = acq.search(signal_slice,
                          delta_freq_step=500, 
                          delta_freq_range=(-5000,5000)) # increase to 100k for HackRF w/o external clock
    
    signal_strength[i] = np.median(np.abs(acq_data['correlator']**2))  # average power

best_prns = np.argsort(-signal_strength)[:5] + 1
print(f'Top-5 PRNs: {best_prns}')

plt.subplots(figsize=(10,2))
plt.bar(prn_ids, signal_strength * 100, tick_label=prn_ids)
plt.xlabel('Satellite [PRN #]')
plt.ylabel('Signal Power [%]')
plt.show()

Top-5 PRNs: [21 16 27  1 32]

Figure 3: Received signal power on all 37 C/A codes identified by their PRN number during first 50 ms of the recording.

If you want to know more about C/A codes, how we detect it in the recording, what the four demodulation parameters exactly mean and how we have to search for them, then wait for walkthrough step 3 to be posted… In the meantime, have a look into the source code of the acquisition and C/A code module.

4 Signal Tracking and Demodulation

In order to maintain reception over the time of the recording, we constantly need to readjust to modulation parameters for each satellite. As shown in Figure 4, the satellite’s signal is lost in a fraction of a second if not done so.

Code

# get initial doppler shift and code delay parameters using acquisition procedure
acq = gpswt.Acquisition(prn_id=16, sampling_dt=baseband.sampling_dt)
acq_data = acq.search(baseband.get_interval(to_time=0.005))  # first periods only
delta_freq_start = np.median(acq_data['delta_freq'])
delay_start = np.median(acq_data['delay'])

# tracking with disabled delay-locked loop (DLL)
tracker = gpswt.Tracking(prn_id=16, sampling_dt=baseband.sampling_dt, 
                         dll_guided=False, dll_ki=0.0,  # disable DLL               
                         delta_freq_start=delta_freq_start,
                         delay_start=delay_start)
tracking_data = tracker.track(baseband.get_interval(to_time=1.0))

amplitudes = tracking_data['correlator']
t = tracking_data['time_code_start']

fig, ax = plt.subplots(figsize=(10, 2.*1))
ax.plot(t, abs(amplitudes)**2 * 100, 'b')
ax.set_xlabel('Receiver Time $t$ [s]') 
ax.set_ylabel('Signal Power [%]')
plt.show()

Constant readjustment could be performed with the acquisition procedure shown above, by blindly searching the optimal parameters for every millisecond of the signal individually. But there are some major drawbacks using this approach.

Exploiting the fact that the parameters change continuously, we can keep them optimized using control loops: One loop locking on the frequency, one on the code delay, and one on the phase. This procedure is called tracking of the signal and has several advantages compared to the acquisition algorithm: It is much faster since it does not start search blindly from scratch, yields more accurate parameter values, and is more robust to noise.

The result of tracking is a

high-precision code delay value (for our analysis here we achieve less than 10 nanoseconds) for accurate time-of-flight estimation, and
the CDMA-demodulated signal, which allows to decode the telemetry messages sent by each satellite.

Tracking of PRN 16 is shown in Figure 5. After a short time period at the beginning where the control loops lock-in on the optimal values, the parameters develop steadily and signal power is maintained.

Code

# get initial doppler shift and code delay parameters using acquisition procedure
acq = gpswt.Acquisition(prn_id=16, sampling_dt=baseband.sampling_dt)
acq_data = acq.search(baseband.get_interval(to_time=0.005))  # first periods only
delta_freq_start = np.median(acq_data['delta_freq'])
delay_start = np.median(acq_data['delay'])

tracker = gpswt.Tracking(prn_id=16, sampling_dt=baseband.sampling_dt,
                         delta_freq_start=delta_freq_start,
                         delay_start=delay_start)
tracking_data = tracker.track(baseband.get_interval(to_time=1.0))

amplitudes = tracking_data['correlator']
t = tracking_data['time_code_start']

fig, ax = plt.subplots(4, 1, sharex=True, figsize=(10, 2.*4))
ax[-1].set_xlabel('Receiver Time $t$ [s]') 

ax[0].plot(t, abs(amplitudes)**2 * 100, 'b')
ax[0].set_ylabel('Signal Power [%]')

ax[1].set_ylabel('Code Delay $\\tau$ [$\mu$s]')
ax[1].plot(t, tracking_data['delay_used'] * 1e6, 'b', label='Used')
ax[1].plot(t, (tracking_data['delay_used'] + tracking_data['discriminator_tau']) * 1e6,
           'g', label='Measured', alpha=0.5)
ax[1].plot(t, tracking_data['delay'] * 1e6, 'r', label='Estimated')
ax[1].legend()    

ax[2].plot(t, tracking_data['delta_freq'], 'm')
ax[2].set_ylabel('Doppler Shift $\Delta f$ [Hz]')
    
ax[3].plot(t, tracking_data['phi'], 'g')
ax[3].set_ylabel('Phase $\Phi$ [rad]')
plt.show()

Over the shown time period, the code delay needs to be increased constantly. Equally as for the negative Doppler shift frequency, this is a consequence of the fact that this satellite is heading away from us and its distance is thereby increasing, i.e. its signal arrives more and more delayed over time. Over the shown time period of 1s, we observe an increase of approx. 1.5$\mu$s in delay. Multiplied by the speed of light, this corresponds to a speed of roughly 450 m/s. (Notice, that this value does not reflect the satellite’s speed, but only to its component along the line of sight from us.)

BPSK Demodulation

Each GPS L1 satellite sends navigation messages (LNAV) which contains the telemetry data. These messages are modulated onto the CDMA signal using binary phase-shift keying (BPSK). We see in Figure 6, after the lock-in phase of the tracking, that a clean digital signal can be observed. Obviously, the BPSK demodulation yielding the binary information is straightforward.

Code

fig, ax = plt.subplots(figsize=(10, 2.))
ax.plot(t, np.real(amplitudes) / np.abs(amplitudes), 'c')
ax.set_ylabel('Norm. In-Phase $I \, / \, |A| $')
ax.set_xlabel('Receiver Time $t$ [s]')

plt.show()

Figure 6: CDMA-Demodulated telemetry signal of a satellite.

A different method allowing to visualize digital modulation of a signal over longer time periods are constellation diagrams. In Figure 7, we observe that later CDMA demodulated samples (blue, green, yellow), where signal tracking locked after startup (purple), fall into two well distinguishable spots. These are the two states encoding a bit in the BPSK scheme. That means that the signal tracking and CDMA demodulation is reliable over the first 10s, and BPSK demodulation will work flawlessly.

Code

tracker = gpswt.Tracking(prn_id=16, sampling_dt=baseband.sampling_dt, 
                         delta_freq_start=acq_data[0]['delta_freq'],
                         delay_start=acq_data[0]['delay'])
tracking_data = tracker.track(baseband.get_interval(to_time=10))  # 10 s interval
amplitudes = tracking_data['correlator']

fig, ax = plt.subplots(figsize=(10,4))
ax.scatter(np.real(amplitudes), np.imag(amplitudes), 
           marker='.', alpha=0.3, c=range(len(amplitudes)))
ax.set_aspect('equal', 'box')
ax.set_xlabel('In-phase $I$')
ax.set_ylabel('Quadrature $Q$')
ax.grid(which='major')
plt.show()

Figure 7: Constellation diagram of CDMA-demodulated signal during the first 10 seconds. Color indicates reception time from start (purple) to end (yellow). A point represents an IQ sample of one pseudo-symbol or C/A code period.

Track and Demodulate All

We can now get the navigation messages being sent by performing all steps shown above for each satellite:

Acquisition, to find initial demodulation parameters.
Tracking, to follow the parameters over long recording time and get CDMA demodulation.
BPSK demodulation, to get the encoded binary stream of the navigation messages.

Code

print(f'Analyzing signal of length {baseband.duration():.1f} s:')
pseudo_symbols = {}
tracking_data = {}
for prn_id in best_prns:
    print(f'PRN {prn_id}:', end='')

    print('  Acquisition... ', end='')
    acq = gpswt.Acquisition(prn_id, sampling_dt=baseband.sampling_dt)
    acq_data = acq.search(baseband.get_interval(to_time=0.05))  # check first 50ms = 50 code periods

    delta_freq_start = np.median(acq_data['delta_freq'])
    delay_start = np.median(acq_data['delay'])

    print(f'  Tracking... ', end='')
    tracker = gpswt.Tracking(prn_id, sampling_dt=baseband.sampling_dt,
                             delta_freq_start=delta_freq_start,
                             delay_start=delay_start)
    tracking_data[prn_id] = tracker.track(baseband, progress=False)  # track full signal length

    print(f'  BPSK Demodulation... done.')
    pseudo_symbols[prn_id] = \
        np.real(tracking_data[prn_id]['correlator']) > 0  # distinguish between 0 and 1

Analyzing signal of length 84.3 s:
PRN 21:  Acquisition...   Tracking...   BPSK Demodulation... done.
PRN 16:  Acquisition...   Tracking...   BPSK Demodulation... done.
PRN 27:  Acquisition...   Tracking...   BPSK Demodulation... done.
PRN 1:  Acquisition...   Tracking...   BPSK Demodulation... done.
PRN 32:  Acquisition...   Tracking...   BPSK Demodulation... done.

This is rather lengthy calculation, taking almost 10 minutes on my laptop for the given length of the recording. Certainly, the algorithm’s are not at all optimized for speed, but are roughly 20 times faster than the acquisition procedure is.

If you want to know more about control loops, signal processing, demodulation, wait for the walkthrough step 4 to be published. In the meantime, have a look into the source code of the tracking module.

Code

# to avoid running tracking again, save the output for later experimentation...

import pickle
# if pseudo_symbols and tracking_data:  # save the data if present
#     with open('psymbols.pickle', 'wb') as f:
#         pickle.dump(pseudo_symbols, f)
#     with open('trackdata.pickle', 'wb') as f:
#         pickle.dump(tracking_data, f)

with open('psymbols.pickle', 'rb') as f:  # load it from file if available
    pseudo_symbols = pickle.load(f)
with open('trackdata.pickle', 'rb') as f:
    tracking_data = pickle.load(f)
best_prns = list(pseudo_symbols.keys())

5 Telemetry Decoding

The telemetry data sent in navigation messages by all satellites contain information about their orbital parameters as well as clock error information of the on-board atomic clocks. Both are required for calculating a satellite’s position very accurately at every moment in time.

The BPSK demodulation done above yields a binary stream. Its constituents are called pseudo-symbols. The actual information carrying bits, used for encoding the telemetry message, consist of 20 identically repeated pseudo-symbols. Messages sent by the satellite are encoded in a data structure called frame which consists of five subframes. To get the frame data, we have to perform the following steps:

Bit synchronization: Find out where in the stream of pseudo-symbols a bit starts.
Bit decoding: Decode the pseudo-symbols and generate the stream of bits.
Frame synchronization: Find the start of a frame in the stream of bits.
Frame dissection: Read the frame data structure, check its validity using checksums, and represent it in a usable form.

Let’s do that for all satellites:

Code

telemetry = {}
for prn_id in best_prns:
    print(f'PRN {prn_id}')
    
    bits, bit_start_symbol = gpswt.telemetry.bit_synchronize(pseudo_symbols[prn_id])
    bits, subframe_start_bit = gpswt.telemetry.subframe_synchronize(bits)
    telemetry_data = gpswt.telemetry.subframe_dissect(bits, subframe_start_bit, bit_start_symbol)
    
    telemetry[prn_id] = telemetry_data
    print()

PRN 21
Bits: synced, avg. pseudo-symbol errors per bit: 0.002
Subframes: synced, found 14, valid 13
Dissection: 13 subframes decoded.

PRN 16
Bits: synced, avg. pseudo-symbol errors per bit: 0.003
Subframes: synced, found 14, valid 13
Dissection: 13 subframes decoded.

PRN 27
Bits: synced, avg. pseudo-symbol errors per bit: 0.001
Subframes: synced, found 14, valid 13
Dissection: 13 subframes decoded.

PRN 1
Bits: synced, avg. pseudo-symbol errors per bit: 0.002
Subframes: synced, found 14, valid 13
Dissection: 13 subframes decoded.

PRN 32
Bits: synced, avg. pseudo-symbol errors per bit: 0.004
Subframes: synced, found 14, valid 13
Dissection: 13 subframes decoded.

We see from the output that the telemetry data of all five satellites was successfully decoded. While a few pseudo-symbols are reported to be faulty, the signal quality was good enough that almost all subframes were valid, i.e. contain no bit errors.

Let’s have a look at the telemetry data of one of the satellites. The first five received subframes read

Code

for i in range(5):  # show a full frame = 5 subframes
    print(telemetry[16][i])

{'start_symbol': 3964, 'subframe_id': 1, 'integrity': 0, 'time_of_week': 52116, 'clock': {'week_number': 213, 'sv_health': 0, 'T_GD': -1.0244548320770264e-08, 't_oc': 57600, 'a_f2': 0.0, 'a_f1': 3.183231456205249e-12, 'a_f0': -0.0005030031315982342}}
{'start_symbol': 9964, 'subframe_id': 2, 'integrity': 0, 'time_of_week': 52122, 'ephemeris': {'C_rs': 131.75, 'dn': 4.2019607428169266e-09, 'M_0': -3.0335492879613675, 'C_uc': 1.642853021621704e-06, 'e': 0.01333794859237969, 'C_us': 9.84780490398407e-06, 'sqrtA': 5153.60133934021, 't_oe': 57600}}
{'start_symbol': 15964, 'subframe_id': 3, 'integrity': 0, 'time_of_week': 52128, 'ephemeris': {'C_ic': -9.313225746154785e-09, 'Omega_0': 0.6220877159703716, 'C_is': -2.3096799850463867e-07, 'i_0': 0.965627893317696, 'C_rc': 192.03125, 'omega': 0.7866631281797375, 'Omega_dot': -7.618888785870861e-09, 'I_dot': 6.321691895270685e-11}}
{'start_symbol': 21964, 'subframe_id': 4, 'integrity': 0, 'time_of_week': 52134}
{'start_symbol': 27964, 'subframe_id': 5, 'integrity': 0, 'time_of_week': 52140}

The subframes data structure contains:

time_of_week: exact time of sending of this subframe by the satellite in seconds, as measured by the satellite’s atomic clock,
clock: satellite clock correction data, since the atomic clocks need still some correction to be accurate enough for our purposes,
ephemeris: parameters required to calculate the satellite orbit positions at any time,
meta data added by the telemetry decoder (e.g. the subframe’s start pseudo-symbol index, start_symbol, in the sequence of recorded symbols).

We can observe that a subframe is transmitted every 6 seconds (steps in time_of_week), and that one subframe is encoded using 6000 pseudo-symbols (steps in start_symbol).

If you want to know more about synchronization, decoding, dissetion and the GPS L1 navigation message data structures, wait for the walkthrough step 5 to be published. In the meantime, have a look into the source code of the telemetry module.

6 Pseudorange Measurement

In order to calculate our position, we need to perform measurements of so called observables. The one we will consider here is called pseudorange. It represents a measurable quantity which is related to the distance between the receiver and a satellite. The pseudorange to satellite i is defined by

\[ R_i = c \, (t_{i} - T_{i}) \]

with the speed of light in vacuum $c$, the receiving time $t_i$ of a signal as measured by the receiver, and the send time $T_i$ of this signal as measured by satellite i’s atomic clock in GPS time. Therefore, measuring $t_i$ and decoding $T_i$ from the telemetry allows us to calculate $R_i$.

The speed of light $c$ is identical to the speed of the radio wave propagating through space. The term $(t_i - T_i)$ is a time difference between sending and receiving the signal, which is related to the time-of-flight. The pseudorange $R_i$ is their product, speed times time difference, and therefore results in a distance related to the distance between satellite and receiver.

The reason why $R_i$ is not the actual distance from the receiver to the satellite but only related: While $T_i$ is very accurate thanks to the atomic clocks and their correction information available through telemetry, the cheap receiver clock built into our SDR, on which we base our measurements of $t_i$, is not accurate at all. Furthermore, the speed of light might not correspond to the radio wave’s effective speed of propagation since the atmosphere’s influence it, see below. These discrepancies make the pseudorange deviate from the actual, physical distance and have to be taken into account later.

Let’s have a look at pseudoranges of PRN 16 during the first few seconds:

Code

pseudo_range = gpswt.observables.pseudo_ranges(telemetry[16], tracking_data[16], 
                                               observation_interval=1.0)
for k in range(4):
    delta = pseudo_range['pseudo_range'][k] - pseudo_range['pseudo_range'][0]
    print(f"t = {pseudo_range['receive_time'][k]:.6f} s, "
          f"T = {pseudo_range['send_time_gps'][k]:.6f} s, "
          f"R = {pseudo_range['pseudo_range'][0]/1e3:.0f} km + {delta:.2f} m")

t = 0.000000 s, T = 52106.036490 s, R = -15620996756 km + 0.00 m
t = 1.000000 s, T = 52107.036488 s, R = -15620996756 km + 543.89 m
t = 2.000000 s, T = 52108.036486 s, R = -15620996756 km + 1085.70 m
t = 3.000000 s, T = 52109.036485 s, R = -15620996756 km + 1626.99 m

A comparison between the receiver time $t_i$ and the send time $T_i$ reveals that they are roughly in sync, but off by almost a day. In the last shown digits of $T_i$, we can observe that the signal received at the integer values $t_i$ are sent earlier and earlier, explained by an increase in distance. The pseudorange is negative, a non-sensical value for a distance, due to the wrong receiver time offset by almost a day. But what can be taken more seriously are the changes reported after the + sign. We can read off a receding speed of approx. 540 m/s of the satellite (again only in line of sight).

Let’s measure the pseudoranges for all satellites as we need them in our next steps.

Code

position_fix_time_interval = 0.5
print(f'Measurement Time Interval: {position_fix_time_interval} s\n')

print(f'Measuring pseudoranges for PRN: ', end='')
pseudo_ranges = {}
for prn in best_prns:
    print(f'{prn} ', end='')
    pseudo_ranges[prn] = gpswt.observables.pseudo_ranges(telemetry[prn], tracking_data[prn], 
                                                         position_fix_time_interval)

Measurement Time Interval: 0.5 s

Measuring pseudoranges for PRN: 21 16 27 1 32

If you want to know more about the measurement of the nanosecond accurate time of reception, pseudoranges, clock corrections, then wait for the walkthrough step 6 to be published. In the meantime, have a look into the source code of the observables module.

7 Satellite Orbital Position Calculation

In order to infer our position from the measured pseudoranges, we additionally need to know the position of each satellite at the time of sending a signal. As we want to know our position at meter accuracy, we better know the satellites position at the same accuracy at every instant of time.

Kepler’s equation allows to calculate orbits, i.e. position at any time, of objects (satellites) under the gravitational influence of a central mass (earth). Things would unfortunately not work out at the required accuracy, because some assumptions underlying this equation are violated:

The central mass is homogenous and spherical: As we know when climbing mountains, the earth shows some deviations to this.
No other other forces are present: The moon and sun influence the orbit with their gravitational forces, and there is a non-negligible “air” resistance due to particles in space.

These effects need to be considered and put a correction on top of Kepler’s solution. For this, all required information is sent in the ephemeris parameters in the telemetry messages.

Therefore, let’s calculate the position of the satellite sending PRN 16 at GPS time 0, at the beginning of the current week, and 12 hours later:

Code

# extract complete set of ephemeris data from subframes
ephemeris_subframes = [subframe['ephemeris'] for subframe in telemetry[16] 
                       if 'ephemeris' in subframe.keys()]
ephemeris = {}
for subframe in ephemeris_subframes:
    ephemeris.update(subframe)

t = np.array([0, 11.97*60*60])  # GPS time in seconds
pos = gpswt.ephemeris.sat_position(t, ephemeris)

for ti, posi in zip(t, pos):
    print(f't = {ti/3600:>5.2f}h, (x,y,z) = ({posi[0]/1000:>10.3f}, '
          f'{posi[1]/1000:>10.3f}, {posi[2]/1000:>10.3f}) km')

t =  0.00h, (x,y,z) = (-15948.384,   5789.441,  20242.277) km
t = 11.97h, (x,y,z) = ( 15983.496,  -5760.857,  20223.744) km

The shown $(x,y,z)$ coordinate is this satellite’s position in the GPS coordinate system WGS84. It is as earth-centered earth-fixed (ECEF) system, i.e. $(0,0,0)$ is at the earth’s center, the $z$ axis is along to rotation axis, and the whole system is rotating together with the earth, such that a fixed position on the surface does not change its coordinates with time. (Have a look at Ankur’s GPS article for nice illustrations of coordinates systems and orbits.)

We can read off the co-rotating behavior above: GPS satellite’s have a orbital period of approx. 12 hours. Therefore, we would expect to two positions reported to be identical. This is almost the case, were it not for the minus signs change in $x$ and $y$ coordinates. In 12 hours, the earth spins half a rotation around the $z$ axis, changing the coordinates of the same physical position (in a inertial frame) exactly in the observed manner, i.e. turning $(x,y,z) \rightarrow (-x,-y,z)$. In other words: From the point of view of an observer rotating with the earth, the satellite’s second position shows up in the opposite direction of the first one.

From the reported satellite position we derive:

Code

print(f'Sat. distance from earth center: {np.sqrt(np.sum(pos[0]**2))/1e3 :.3f} km')
print(f'Sat. altitude from surface:      {np.sqrt(np.sum(pos[0]**2))/1e3 - gpswt.ephemeris.EARTH_RADIUS/1e3 :.3f} km')

Sat. distance from earth center: 26412.466 km
Sat. altitude from surface:      20034.329 km

Finally, we can calculate the position of every satellite for every point in time in our wave recording.

Code

print(f'Calculating satellite positions of PRN: ', end='')
satellite_positions = {}
for prn in best_prns:
    print(f'{prn} ', end='')
    ephemeris_subframes = [subframe['ephemeris'] for subframe in telemetry[prn] 
                           if 'ephemeris' in subframe.keys()]
    ephemeris = {}
    for subframe in reversed(ephemeris_subframes):  # use newest ephemeris data
        ephemeris.update(subframe)
        
    send_time_gps = pseudo_ranges[prn]['send_time_gps']
    satellite_positions[prn] = gpswt.ephemeris.sat_position(send_time_gps, ephemeris)

Calculating satellite positions of PRN: 21 16 27 1 32

If you want to know more about solving the Kepler’s equation, its required corrections, coordinate systems and illustration of orbits, then wait for the walkthrough step 7 to be published. In the meantime, have a look into the source code of the ephemeris module.

8 Position Fixing

In this last step of our journey, we finally will get our location, a position fix. From our measured pseudoranges and the satellite’s positions, we should be able to infer it. This process is called PVT solving, since we get position (P) and velocity (V) of the receiver and an accurate time (T).

The PVT solver algorithm is sketched in Figure 8. Starting from an initial estimate of P, V and T, we need a model which generates the observables one would measure in this situation (model observables). Then a comparison of this model output with the actually measured quantities (measured observables) determines the mismatch and estimates how much the initial PVT estimate was off to get a smaller mismatch. This process is iteratively repeated until the model and measurement match good enough. We arrive at a PVT fix.

flowchart LR
  E[/PVT Estimate/] --> A[Measurement\nModel] 
  A --> D[Model\nObservables] --> B[Observables Mismatch]
  X[/Measured\nObservables/] --> B
  B --> F[PVT Correction]
  F --> E

Figure 8: PVT solving process from measured observables as input to the PVT estimate as output.

In our case, the observables are the pseudoranges, and the model contains the satellite positions. It calculates euclidean distances between them and the current position estimate. The earth center $(0,0,0)$ as initial position is good enough as a starting point and the process converges to the true position and time. Here, we neglect velocity V and assume a stationary receiver. The iterative procedure reduces to a simple least-squares optimization problem which is solved by the gradient descent method. (See Ankur’s article for illustration and derivation of the PVT solving process.)

Position Analysis

Let’s do that with our data and find the receiver’s GPS coordinates:

Code

pvt_solution = gpswt.pvt_solver.position_fix(pseudo_ranges, satellite_positions)
pvt_solution_mean = np.mean(pvt_solution, axis=0) 
print(f' Average Position: (x,y,z) = ({pvt_solution_mean[0]/1e4:.0f}X.XXX,'
      f' {pvt_solution_mean[1]/1e4:.0f}X.XXX,'
      f' {pvt_solution_mean[2]/1e4:.0f}X.XXX) km')

 Average Position: (x,y,z) = (433X.XXX, 56X.XXX, 464X.XXX) km

(Since my home is such a calm place, i censored some of the digits in order to keep it like that.) Hmm, interesting… No clue wether that makes sense, since my I am not very fluent in WGS84 coordinates. Let’s transform them to better known ellipsoid coordinates of longitude and latitude.

Code

position_ellipsoid = gpswt.coordinates.ecef_to_ellipsoid(*pvt_solution_mean[0:3])
print(f'Average Position : {position_ellipsoid[0]:0.1f}XXX°N,'
      f' {position_ellipsoid[1]:.1f}XXX°E,' 
      f' {position_ellipsoid[2]:.2f} m Altitude')

Average Position : 46.9XXX°N, 7.4XXX°E, 656.40 m Altitude

Looks not so bad, but the big question is, how good is it? Entering it in Google maps marks a spot which is only 2 meters away from my balcony where the receiver sits! The altitude matches within 10 meters with what the ‘GPSTest’ app on my smartphone shows. We have found the receiver’s location in a volume of space spanning tens of thousands of kilometers in diameter! Isn’t it cool? 😀

The above position is a long time average over the recording. In Figure 9, we can see how the position fix varies over time. Since the receiver was stationary, variations have to be attributed to errors in measurement and PVT model as discussed in a moment.

Code

delta = pvt_solution - np.mean(pvt_solution, axis=0)
delta = np.vstack([gpswt.coordinates.ecef_to_local(pt, 
                                             position_ellipsoid_mean[0], 
                                             position_ellipsoid_mean[1], 
                                            ) for pt in delta[:, 0:3]])
receiver_time = pseudo_ranges[best_prns[0]]['receive_time']

plt.figure(figsize=(15,3))
plt.plot(receiver_time, delta[:, 0], label='East')
plt.plot(receiver_time, delta[:, 1], label='North')
plt.plot(receiver_time, delta[:, 2], label='Up')
plt.ylabel('Position Fix Variation [m]')
plt.xlabel('Receiver Time [s]')
plt.ylim([-20,30])
plt.legend()
plt.show()

Figure 9: Position fix variation from mean value, in a east-north-up (ENU) coordinate system localized at the mean position.

Are you wondering at how far the satellites are away from us, too?

Code

print('Distance from receiver to satellite')
for prn in best_prns:
    delta = satellite_positions[prn][0] - pvt_solution[0][0:3]  # take first observation only
    distance = np.sqrt(sum(delta**2))
    print(f'  PRN {prn:>2}: {distance/1000:.3f} km')

Distance from receiver to satellite
  PRN 21: 21856.677 km
  PRN 16: 23475.285 km
  PRN 27: 20485.925 km
  PRN  1: 23915.182 km
  PRN 32: 24205.638 km

It’s quite mind-boggling that we received the analyzed signal over more than 20000 km distance, almost two times the earth’s diameter!

How accurately have we been able to measure these distances? The answer give the pseudorange residuals: The discrepancy between the measured and modelled pseudoranges are shown in Figure 10 for every satellite. For some of them we see less than 3 meters, for others 5 to 10 meters difference (corresponding to 10-30 nanoseconds time-of-flight error) that couldn’t be reduced by a more optimal PVT solution. These residuals cannot be explained by the model. Therefore, we can roughly assume that the above reported satellite distances are accurate to less than 10 meters, which is quite fascinating.

There is a reasons for the variation and residuals: Our model does not incorporate atmospheric effects. They influences radio propagation speed and thereby our estimation of the pseudorange. Some of them could be addressed: Ionosphere and troposphere data would be available in telemetry. (Feel free to send me a pull request to add this feature 😉) Other more localized variations of the atmosphere cannot be modelled and will remain as residuals. Taking into account a better model and other satellite system would allow a more stable and accurate position fix.

Code

pseudo_range_errors = np.vstack([gpswt.pvt_solver._position_fix_error_function(x, ti, pseudo_ranges, satellite_positions) 
                                 for ti, x in enumerate(pvt_solution)])
plt.figure(figsize=(15,3))
for prn_index, prn in enumerate(best_prns):
    plt.plot(receiver_time, pseudo_range_errors[:, prn_index], label=f'PRN {prn}')
plt.ylabel('Pseudorange Residual [m]')
plt.xlabel('Receiver Time [s]')
plt.legend()
plt.show()

Figure 10: Pseudorange residuals, i.e. not explained mismatch measured vs. model pseudorange at optimal position/time solution, for all satellites.

Receiver Clock Analysis

In addition to the position, we also get a very accurate GPS time in the PVT solving process. With this we can analyze how accurate our clock is. As mentioned in Section 2, a Siglent function generator was used as a high-precision reference clock. The spec sheet reports a precision of better than $10^{-6}$, i.e. less than a microsecond off per second. Let’s have a look wether that’s true:

Code

from scipy import stats

# bias = difference between receiver time and GPS time
clock_bias = receiver_time - pvt_solution[:, 3]
fit_bias = stats.linregress(receiver_time, clock_bias)

print(f'Average clock drift: {(fit_bias.slope) * 1e6:0.2f} microseconds/second')

Average clock drift: 0.50 microseconds/second

Indeed we find that the promise is hold, it is only slightly too fast. What about the long term stability? In Figure 11, we see slight variations of less than $\pm 0.5\mu s$ within 80 seconds recording period. The clock is rather stable.

Code

plt.figure(figsize=(15,2.5))
plt.plot(receiver_time, (clock_bias - receiver_time*fit_bias.slope - fit_bias.intercept) * 1e6)
plt.ylabel('Non-Linear Clock Bias [$\mu s$]')
plt.xlabel('Receiver Time [s]')
plt.show()

Figure 11: Receiver clock stability. The clock bias is shown with its linear component (clock drift) removed.

Position and Satellite Visualization

To conclude, we visualize our location and satellites in Figure 12. Obviously we only see satellites which are above the horizon, i.e. a direct line of sight is not blocked by the earth. Furthermore, we can observe the effect of the house wall the receiver was placed in front of: No satellites are received in the northern sky of the position fix as seen best in the side view.

And again, if you want to know more about PVT solving, then wait for the walkthrough step 8 to be published. In the meantime, have a look into the source code of the PVT solver and coordinates module.

I hope you enjoyed the journey! Feel free to contact me and please leave your comments below. If you appreciated this article, I am very happy to be invited to a coffee… 😊

Code

fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(projection='3d')
ax.set_proj_type('ortho')
    
# plot earth in spherical approximation
R = 6378137 / 1000 # [km] earth radius
u, v = np.meshgrid(np.linspace(0, 2*np.pi, 25),  # wireframe for every 15°
                   np.linspace(0, np.pi, 13)) 
x = R * np.cos(u) * np.sin(v)
y = R * np.sin(u) * np.sin(v)
z = R * np.cos(v)
ax.plot_wireframe(x, y, z, alpha=0.2, color='blue')

# add equator and null meridian
u = np.linspace(0,2*np.pi, 25)
ax.plot(R * np.cos(u), R*np.sin(u), 0*u, color='black')
ax.plot(R * np.cos(u), 0*u, R*np.sin(u), color='black')

# draw equatorial plane
x, y = np.meshgrid(np.linspace(-5*R, 5*R, 20), 
                   np.linspace(-5*R, 5*R, 20))
z = 0 * x
ax.plot_surface(x, y, z, alpha=0.2, color='gray')

# plot our position in red
ax.stem(*(pvt_solution[:, 0:3]/1000).transpose(), 
        markerfmt='or', bottom=0, orientation='z')

# plot satellites position
for prn in best_prns:    
    ax.stem(*(satellite_positions[prn][-2:-1]/1000).transpose(), 
            markerfmt='og', bottom=0, orientation='z')
    ax.text(*(satellite_positions[prn][-1]/1000 + [0,0,1000]), str(prn), None)

ax.set_box_aspect([ub - lb for lb, ub in (getattr(ax, f'get_{a}lim')() for a in 'xyz')])  # set equal aspect ratio
ax.set_xlabel(('x [km]'))
ax.set_ylabel(('y [km]'))
ax.set_zlabel(('z [km]'))
plt.show()

plt.figure(fig)
ax.view_init(90,90)
fig.set_size_inches(8,8)
plt.show()

plt.figure(fig)
ax.view_init(0,90)
fig.set_size_inches(8,8)

(a) Earth with indicated equator and prime meridian, our position (red), and satellites in WGS84 coordinate system. Stems are drawn from the equatorial plane (gray).

References

navipedia GPS L1 C/A
Wikipedia GPS L1 C/A
gnss-sdr GPS L1 C/A
gps.gov IS-GPS-200 GPS specification (PDF)
Calculating Position from Raw GPS Data, blog post by Ankur Mohan

See also more specific technical references in the source code of the gps_walkthrough python module.