GNSS tri-band, quad-constellation sky recording

I’d like to share my first tri-band GNSS sky recording along with scripts for acquisition and tracking.

The recording is 50 seconds long with three channels of simultaneous 2-bit complex samples at 70 Msa/s (approx. 50 MHz usable bandwidth per channel) centered on L1, L2, and E5 (~1191 MHz). It’s therefore guaranteed to contain a complete data frame for each satellite—the length of an ephemeris frame is 30 seconds for all the GNSS systems, at least for the legacy data formats (though Galileo cuts this in half when using dual frequency via interleaving). An autonomous solution can be obtained for GPS, GLONASS, and Galileo (all of which have at least 4 satellites present), and nearly for BeiDou (only 3 satellites present).

Here is the recording in tcpdump / libpcap format. The utility packet2wav_3ch extracts one of the three RF channels from this packet stream and converts to a raw stream of signed 8-bit values for a downstream software receiver or other tool.

gnss-20170427-L1L2L5.pcap (5.4 GB)

By my count there are 34 satellites in this recording that are emitting 64 distinct GNSS signals (more if you include P/Y and M and CS/PRS codes and count I and Q separately). The satellites with active signals are as follows:

GNSS system PRNs or FDMA channels
GPS 2 4 5 12 13 15 18 20 21 25 26 29
GLONASS -7 -5 -3 -2 1 2 3 5 6
Galileo 7 12 14 19 20 24 26
BeiDou 7 14 34
SBAS 133 135 138

and the signals themselves:

GNSS system Band PRNs or FDMA channels
GPS L1 2 4 5 12 13 15 18 20 21 25 26 29
GPS L2C 4 5 12 25 26 29
GPS L5 25
GLONASS L1 -7 -5 -3 -2 2 5 6
GLONASS L2 -7 -5 -3 -2 1 2 3 5 6
GLONASS L3OC 9 26
Galileo E1 7 12 14 19 20 24 26
Galileo E5b 7 12 14 19 24 26
Galileo E5a 7 12 14 19 24 26
BeiDou B1I 7 14 34
BeiDou B2I 14
SBAS L1 135
SBAS L5 133 135 138

Below are plots of some of the acquisition metrics (but I’ve not yet automated any sort of detection threshold). Click for larger image:

acquisition-composite

The L5/E5 channel should support tracking of the full Galileo E5 signal (E5ab). I haven’t implemented this tracking mode yet, but it involves using a combination of a pure PLL for the pilot signals and a Costas loop for the data signals, optimally weighted and combined. (These would be the “C8X” and “L8X” observables in the taxonomy of RINEX 3.0.)

A makefile is provided in the GNSS-DSP-tools repository that will download the recording, run all the acquisition routines, and track the signals. Acquisition of GPS L1 is quick, but some other signal types are slower because they search all PRNs exhaustively and with a fairly high sensitivity.

Some miscellaneous notes on the signals:

  • I’m not sure why SBAS 133 and 138 are not showing up on L1. They have strong signals on L5, and other recordings in the past have shown them on both L1 and L5 with no trouble. Perhaps there’s some flaw in the acquisition related to the higher data rate on SBAS (although I do only 1 ms integrations, so it should not matter).
  • Same for GLONASS channels 1 and 3 on L1 vs. L2.
  • BeiDou PRNs 7 and 34 are apparently not transmitting on B2I. (PRN 34 is probably transmitting the BeiDou phase 3 L1 MBOC signal, incidentally.)
  • The GLONASS L3OC PRNs nominally run from 1 to 24, corresponding to the slot number, at least according to this GPS World article. Not sure why there’s a PRN 26 being transmitted by one of the GLONASS-K satellites. As I recall, there were other PRNs outside of the range 1–24 also used by GLONASS-K; perhaps the code assignments have changed over time.

The antenna is a homemade 3-turn helical with highpass filter and LNA, both from Mini-Circuits. It’s pretty much the simplest thing that could possibly work while still having acceptable bandwidth (1150–1610 MHz) and acceptable spatial pattern and circular polarization. Some slight interference from LTE is visible in the waterfall plots (using baudline), but it’s weak enough that the GNSS signals are not affected. Probably better preselection filtering is needed, or a stronger LNA, or both.

GNSS Firehose update

Some updates on the GNSS Firehose system as it approaches general usability:

  • now supports three 50 MHz RF channels (nominally covering L1/L2/E5), 840 Mbit/s total payload
  • supports command and status over Ethernet (in addition to UART)
  • firmware now in C, running on a RISC-V soft CPU (picorv32 from Clifford Wolf)
  • change over to KiCad for schematic and PCB layout

The enclosure is the same as before, a Hammond 1455L1201, with fancier artwork this time:

f-perspective

f-front

f-back

A run of 20 boards was assembled, with a very pleasant yield of 100% (thanks, Sparqtron!). Here’s a shot of some of them. They have yet to receive the small 3D-printed plastic cap for the TCXO. That turns out to be quite important to reduce air currents near the device—it dramatically improves the performance of the local clock, so that, for a reasonably strong signal, there’s only the occasional excursion of 1 Hz or so over timescales of a few seconds. (I should make an Allan-variance plot from one of the PLL tracks.)

f-board-6up

GNSS Firehose status, example L1/L2 sky recording

Just a quick update on the GNSS Firehose digitizer project. I’ve decided to get a few systems professionally assembled; they will be similar to this prototype unit:

gnss-firehose-proto-2

Here is some sample data, a sky recording taken on May 6 at around 13:38 UTC:

gnss-3.dat (743 MByte)

This represents 10.2 seconds of data simultaneously sampled on bands centered at ~1584.8 MHz and ~1227.7 MHz; each channel has a useful bandwidth of about 50 MHz. The format is raw Ethernet packets as written by tcpdump. The packet2wav utility unpacks the various sample streams from these packets and checks timestamps.

To use this file with my software receiver, here’s a command to acquire all the GPS L1 C/A signals in view:

$ <gnss-3.dat packet2wav | ./acquire-gps-l1.py /dev/stdin 69984000 -9334875
prn   1 doppler  3200.0 metric  3.45 code_offset  863.2
prn   2 doppler  1800.0 metric  1.45 code_offset  603.4
prn   3 doppler  4200.0 metric  2.33 code_offset  456.6
prn   4 doppler  1000.0 metric 10.23 code_offset  134.4
prn   5 doppler -3600.0 metric  1.54 code_offset   57.2
prn   6 doppler -4800.0 metric  1.46 code_offset  433.3
prn   7 doppler  -800.0 metric  1.45 code_offset  222.5
prn   8 doppler  3000.0 metric  1.46 code_offset  488.3
prn   9 doppler   400.0 metric  1.57 code_offset  362.4
prn  10 doppler  3400.0 metric  1.54 code_offset  506.3
prn  11 doppler  1000.0 metric  7.62 code_offset  728.8
prn  12 doppler  -400.0 metric  1.46 code_offset   93.4
prn  13 doppler -2600.0 metric  1.44 code_offset  595.4
prn  14 doppler  -800.0 metric  7.23 code_offset  558.5
prn  15 doppler  -600.0 metric  1.51 code_offset   90.7
prn  16 doppler     0.0 metric  1.43 code_offset  772.2
prn  17 doppler  3200.0 metric  1.44 code_offset  301.7
prn  18 doppler -1200.0 metric  2.20 code_offset  760.5
prn  19 doppler -1600.0 metric  1.59 code_offset  657.4
prn  20 doppler -4400.0 metric  1.50 code_offset  910.4
prn  21 doppler  2200.0 metric  1.42 code_offset  674.6
prn  22 doppler  -800.0 metric  6.54 code_offset  923.6
prn  23 doppler  4200.0 metric  2.08 code_offset   94.4
prn  24 doppler  2200.0 metric  1.55 code_offset  406.4
prn  25 doppler  3000.0 metric  4.29 code_offset  513.2
prn  26 doppler  1600.0 metric  1.47 code_offset  628.4
prn  27 doppler   600.0 metric  1.51 code_offset  135.4
prn  28 doppler   200.0 metric  1.45 code_offset  631.4
prn  29 doppler  1600.0 metric  1.47 code_offset  378.4
prn  30 doppler -2800.0 metric  1.52 code_offset  737.5
prn  31 doppler  3000.0 metric  6.61 code_offset  367.1
prn  32 doppler  2600.0 metric  5.98 code_offset  176.8
prn 133 doppler  1600.0 metric  3.20 code_offset  481.3
prn 135 doppler  1400.0 metric  3.82 code_offset   86.7
prn 138 doppler  1400.0 metric  2.67 code_offset  824.9
$ 

A plot of these acquisition metrics across the various PRNs:

acquisition

Here’s a summary of the visible signals on L1 after acquisition is done across all the GNSS services:

GPS L1 C/A:      1 3 4 11 14 18 22 23 25 31 32 133 135 138
GLONASS L1 C/A:  (none found)
Galileo E1b:     14
Galileo E1c:     14
BeiDou B1I:      11 14

and on L2:

GPS L2CM:        1 3 25 31
GLONASS L2 C/A:  -2 -1 3 5 6
GLONASS L3I:     (none found)
GLONASS L3Q:     (none found)
Galileo E5bI:    14
Galileo E5bQ:    14
BeiDou B2I:      11 14

My L1/L2 antenna, an AeroAntenna AT2775-42, doesn’t quite cover GLONASS on L1. Occasionally a signal on channels -7 or -6 is strong enough to show up, but no luck for this capture. GLONASS L2 is fine though, as are the newer GLONASS L3 CDMA signals (not present in this capture unfortunately). Of course other signals such as L2CL and P will be acquirable and trackable as well, but they are more difficult to acquire blindly.

Just for fun, here are the raw samples from the L1 channel. Clearly not much can be seen from this, but it can be useful as a quick check. Samples are 2 bits (represented on the output stream as two’s-complement 8-bit for processing convenience), alternating I and Q, at 69.984 Msa/s.

$ <gnss-3.dat packet2wav | od -Ad -tx1 | head
0000000 ff 03 03 03 01 ff fd 01 ff 01 03 01 ff 01 03 fd
0000016 01 fd ff fd ff fd ff ff ff fd ff ff ff 03 ff 01
0000032 01 03 01 03 ff 03 03 03 01 ff ff ff ff 03 01 01
0000048 03 01 01 fd fd 01 01 03 01 01 01 01 03 01 01 03
0000064 fd fd fd ff fd 03 fd ff ff 03 01 ff ff fd fd 01
0000080 ff 01 03 ff ff ff 03 03 03 01 ff ff ff ff 01 ff
0000096 fd 01 ff ff 01 01 ff ff ff fd 03 fd ff fd ff fd
0000112 01 01 fd 01 ff fd ff 01 fd ff ff ff 01 01 ff ff
0000128 01 01 01 01 ff 01 03 01 03 ff 01 ff 01 01 01 01
0000144 01 fd fd 01 ff 01 01 ff 01 03 ff 01 01 03 01 01

And the raw samples from channel 2, the downconverter channel looking at L2:

$ <gnss-3.dat packet2wav 2 | od -Ad -tx1 | head
0000000 ff 01 fd 01 01 01 01 ff ff ff fd 01 ff 03 01 03
0000016 03 01 fd 01 fd ff ff 01 ff 03 fd 03 fd ff ff 01
0000032 03 ff 01 fd ff fd 01 01 ff 01 fd 01 01 fd ff ff
0000048 ff 01 03 ff 01 fd 01 ff 01 ff 03 ff ff ff ff ff
0000064 ff ff 01 fd 03 01 03 ff 01 fd ff 01 ff ff fd 01
0000080 ff 01 01 ff 03 01 01 03 03 01 01 ff ff ff 03 ff
0000096 ff ff 01 ff 01 ff 01 01 03 01 01 fd ff fd ff ff
0000112 fd 03 ff 01 ff ff 03 ff 01 ff fd 01 01 fd 01 fd
0000128 01 fd ff ff ff fd fd fd 01 01 01 03 ff 03 01 03
0000144 ff 01 ff 01 01 01 01 ff 03 01 01 ff ff ff ff fd

I plan to add the FPGA and software support for the third RF channel (set to cover L5 by default) over the next few weeks while the boards are being manufactured and assembled. I now have a simple helical antenna suitable for L5, so I should be able to do tri-band experiments.

Pointers to the GitHub repositories containing hardware design and software receiver:

https://github.com/pmonta/GNSS_Firehose (contains packet2wav.c)
https://github.com/pmonta/GNSS-DSP-tools (contains acquire-gps-l1.py etc.)

Squaring receiver for Galileo blind search

I’ve been following the drama of the two recently-launched Galileo FOC satellites. To my knowledge, nothing is yet known about the status of the transmitters, nor have the PRNs been revealed, nor does an acquisition search over all 50 codes in the Galileo ICD yield any signals when the satellites are in view (though PRNs 11, 12, and 19 show up as expected).

So if the satellites are transmitting at all, perhaps it’s with a nonstandard code. One way to check is to take a trip back in time and revisit one of the very oldest GNSS receiver techniques, that of squaring a BPSK signal to recover a tone, then driving a PLL with this tone. (The Macrometer V-1000 receiver used this scheme.)

Now Galileo E1 when viewed in a 4 MHz bandwidth is not BPSK—it is, I believe, a three-level signal, {-1,0,1}, being the sum of equal-power E1-B and E1-C with negligible contribution in quadraphase from PRS. But squaring still works on a signal of this type, though not quite as well as with BPSK.

Here is a waterfall plot with several satellites visible, and indeed Galileo PRN 11 is one of them. Its identity is confirmed by acquiring and tracking in the conventional way, then verifying that the Doppler is identical (modulo the factor of 2 from squaring).

squaring

As of September 20, though, there is still nothing from the new satellites. A predicted Doppler curve can be derived from orbital elements, but there is nothing in the spectrogram at the anticipated Doppler. So we wait.

Code for this “receiver” (all of a dozen lines) is at the GitHub repository (squaring.py). Its output can be piped to the indispensable baudline spectrum analysis tool for interactive viewing.

GNSS stamp collecting

Here’s a fun montage of the various GNSS signals-in-space. I think this accounts for all extant open-access signals that are likely to remain so (ruling out Galileo E6 for example) and that are intended for systems having global coverage (so no QZSS or IRNSS). A few comments on the signals:

  • Each waveform is 400 ms in length and is sampled at 1 ms intervals after code and carrier wipeoff
  • Secondary codes have not been removed
  • Carrier-to-noise ratios vary. The weakest signals are Galileo E5a-I and E5a-Q; by eye there is not much to see, but the acquisition metric is clearly well above threshold, and the histograms are clearly bimodal. My L1/L2 antenna receives some L5 but with about 15 dB of attenuation (!). I really need an antenna that covers down to L5. The GPS L5 signals must have been blisteringly strong to come through as well as they did.
  • The pilot signals with no secondary codes, GPS L2CL and GLONASS L2 P, are shown correctly offset from the (undrawn) centerline. I’m not sure much is known about the data modulation on GLONASS P. Apparently L1 has some data at 50 Hz (no secondary/meander code?) and L2 is unmodulated.
  • So far there are two GLONASS satellites transmitting L3OC CDMA signals. The signals shown are from PRN 30, but PRN 33 is also active and of comparable quality.

The tools for acquiring and tracking these signals are in my GitHub repository:

https://github.com/pmonta/GNSS-DSP-tools

They are command-line-ish and not very polished yet.

track-collection

Height-based multipath mitigation

Here’s a crazy idea which I might as well put on the blog.

Multipath is an important error source for GNSS reference stations. Monuments for antennas are nearly always placed close to the Earth’s surface, so the ground will act as a reflector with a grazing geometry that generates short-delay multipath. Usually other objects contribute as well (nearby buildings or fences for example).

Many solutions exist for multipath mitigation, both at the antenna and correlator levels. Another possible system technique, though, would seem to be to move the antenna upwards, far enough away from local objects that any reflected signals have delays larger than the support of the autocorrelation function for any signal of interest.

Conceptually, an antenna could simply be placed at the top of a tall tower a few hundred meters in height. The tower would ideally be transparent to RF (perhaps of lightweight dielectric construction). Of course there are many practical problems with this, but the environment around the antenna would be nearly ideal.

Another possibility is to place the GNSS antenna on a UAV, which would keep station above the reference monument and several auxiliary sensors (whose location and stability are not critical). The UAV would simultaneously maintain links with visible GNSS satellites (aided by an on-board inertial system) and with sensors on the ground using any of a variety of accurate (~0.1 mm say) short-range ranging techniques. In this way the pristine airborne GNSS signal environment is transferred to the reference monument despite the relative movement.

The UAV could execute certain maneuvers to continuously calibrate its antenna, similar to robot absolute antenna calibrations on the ground. The craft could spin slowly around the vertical axis, or tilt slightly, or both. The attitude from the inertial system would become part of the observation stream to close the calibration loop. By contrast, there is great reluctance to move ground-based GNSS reference antennas to carry out any sort of ongoing calibration program on them. With flyers, continuous monitoring comes for free.

Small rotations and tilts on an airborne platform are impossible to completely avoid, and high-wind situations may force some loss of observing time. But for most of the time, the environment should permit an accurate tie from UAV track to ground network. Depending on the choice of flying craft, several may be needed, spelling each other for charging or refueling. Automatic fleet management will probably have many good solutions over the next few years.

It’s hard to know whether there’s a benefit without more detailed study, but the prevailing trends in GNSS system accuracy seem to be increasing the relative importance of multipath. If we assume progressively better satellite orbit and clock estimates and ionosphere and troposphere sensing, then multipath may well loom as the last remaining large, difficult, uncertain systematic error. (Perhaps a UAV could help with estimating the wet troposphere delay as part of normal operation, to the extent that measurements on the very bottom segment of the troposphere are predictive of the full path.)

Finally, I wonder whether UAVs could help with urban canyons or tree-canopy issues. A surveyor might deal with an awkward situation by tossing a UAV into the air and replacing a bad-GNSS-signal problem with a perhaps easier-to-solve UAV-to-ground-sensor problem using vastly stronger optical or RF links.

Semicodeless P(Y)-code processing using high-rate aiding

Present schemes for semicodeless P(Y) processing assume an autonomous receiver that estimates W-code bits directly from the received signal. Given the limited C/N0 available, naturally the signals are noisy, resulting in squaring loss.

One simple way around this is to use more reliable estimates of the W bits acquired with a medium-gain antenna. These estimates can be published continuously in near-real-time by a suitable Internet-based service, then used by any client receiver, which could be a conventional real-time receiver or a recorded-waveform software receiver.

This need not be all that expensive. All that’s required is approximately one medium-gain antenna per satellite per hemisphere. A one-meter dish of 20 dB gain would reduce squaring loss to practically zero. Of course a reasonable Internet uplink is needed of 480 kbit/second for each satellite. Storage costs could be controlled by retaining the W bits for a limited time (a few weeks or months).

Another possibility, in a world where every receiver is uploading full RF waveforms to a central service, is to sum together the signals from many receivers before estimating the W bits. (If the summing is done prior to detection, this is effectively a phased-array antenna.) Quite a few signals would be needed, though, if each receiver has the usual omni antenna. Better to rely on dedicated medium-gain antennas.

The W-code bits would not be of any use to spoofers since they would be tens or hundreds of milliseconds out of date.

The same trick can be played with other unknown codes, such as the GPS M code or the PRS codes on other GNSS services. The bit rates of the published reference waveforms would be much higher than those of W-code, but perhaps the effort would be repaid by observables that could be obtained in no other way, helping with multipath and ambiguity resolution. In the limit of a totally unknown waveform, this is just VLBI.

Parallel processing of recorded GNSS signals

Most GNSS receivers process signals serially. This is natural for tracking loops based on PLLs and DLLs, as they have a feedback structure. If signals are recorded and stored, however, another viewpoint might be more flexible.

Let’s regard the recorded waveform as a series of chunks of length, say, 5 minutes. All these chunks can be processed in parallel, though at the cost of ambiguities in whole cycles of carrier phase for each chunk. (Let’s assume that acquisition or aiding has already allowed each chunk processor to start with good estimates of code phase and doppler, and that suitable guard intervals allow the tracking loops to converge somewhat in advance of the start of each chunk, so that effectively the chunks overlap a little.) Once all chunks are processed, whole cycles of carrier phase are simply cumulatively summed. This reduces the ambiguity set to the normal case of just a single ambiguity for the whole interval of the satellite pass (assuming no cycle slips or loss of lock).

So an attractive GNSS processing scenario might be:

  • deposit all waveforms in a central place, such as one of the cloud computation environments like Amazon’s S3 and EC2
  • do all processing of interest in parallel, by allocating as many processors as needed; place intermediate results as annotations on a common scoreboard
  • coalesce the results, obtain observables, and post-process

By having the entire waveform accessible at once to a common pool of processors, a kind of annotation-based processing becomes possible. First, acquisition might be performed at fixed intervals, possibly aided by a location estimate and orbit estimates from IGS. Once the file has been annotated with acquisition results, each chunk can be tracked as outlined above. Vector tracking, differencing at the correlator level, quality monitoring, etc. can all be included as additional workflow options.

GPS P code exploration

As a first step towards obtaining GPS P-code observables, it seems prudent to verify that the P code is detectable in a test recording with high C/N0.

Here’s the result with an L1 recording containing a strong signal from PRN 30 (about 50 dBHz). The peak is smeared over an extent somewhat larger than two chips, the result of some residual code doppler. Also the peak is offset by about 1 chip, again the result of residual error from the C/A acquisition (about 0.1 C/A chip). It was great to see it at all though. Certainly it proves out the P code generator.

p_code

The recording is about 1.4 seconds long and contains these data bits:

000100010110111101111011110010001011000011000001100010101110010100100100101

It turns out I was a bit lucky and got a complete TOW word at the end:

0001000101101111011110111100
10001011                        # Preamble
00001100000110                  # TLM Message
0                               # Integrity Status Flag
0                               # Reserved
101011                          # Parity
10010100100100101               # TOW Count (inverted: 01101011011011010)

This TOW count corresponds to Wednesday 19:40:12 in the GPS week, so the first bit of the preamble is at 19:40:06. Adding up the C/A code phase and the previous bits gives the time of the first RF sample of the file as 19:40:05.450822469 or 3375955761914 P code chips. That was the offset given to the script that produced the above plot.

I’ve put the code generators and some preliminary acquisition and tracking software in a new GitHub repository:

http://github.com/pmonta/GNSS-DSP-tools/

I’m sure it’s not as efficient as low-level code in C, but it’s nice to have concise scripts for prototyping, with everything in the Matlab-like python / numpy environment.

New “GNSS Firehose” board

I’ve finally gotten around to updating the GNSS front-end digitizer. Along with a new Ethernet PHY chip (the old one from Vitesse seems to be no longer available), there is an external clock option, an expanded auxiliary header, and a number of small improvements in signal integrity. The external-clock header can accept an external OCXO or rubidium signal, for example; and multiple boards can be driven with a common clock.

new-board

Here’s a spectrum at L1. Despite the poor antenna placement (almost surrounded by tall trees), the GPS C/A signal shows up quite well as a broad peak of 2 MHz bandwidth. There is substantial ripple in the antenna’s ~35 MHz passband, and unfortunately the antenna filtering cuts off around 1595 MHz, so GLONASS signals are suppressed. The signals near 1557 MHz are probably satellite downlinks, and the peak near 1584 MHz is the receiver’s DC spur.

spectrum

The usable alias-free bandwidth of the system is about 50 MHz per channel. At L1 this is enough to cover all the services, from BeiDou B1 starting at 1559 MHz to GLONASS extending to 1610 MHz.

Here’s a C/A correlation peak from this recording (PRN 13). The nice sharp corners are a result of using all of the C/A bandwidth:

peak

Next steps are to clean up the software and HDL and to test the other two channels. See previous blog posts for a pointer to the GitHub repository containing the newly-updated design files.