Abstract
This work explores the use of a low-cost global navigation satellite system (GNSS) antenna array including front-ends and a global navigation satellite system (GNSS) software receiver to receive signals of opportunity (SoO) whose pseudorandom noise (PRN) code is unknown. The front-ends are only loosely synchronized in time and frequency via hardware elements, and precise synchronization or calibration is achieved by using open service global navigation satellite system (GNSS) signals. After calibration, the raw received signals from all antenna elements are added coherently, which allows the pseudorandom noise (PRN) codes of the unknown signals of opportunity (SoO) to be estimated. The pseudorandom noise (PRN) sequences are then fed into a test receiver with a single antenna element that uses the sequences to acquire and track the signals of opportunity (SoO) in a conventional way. The process of chip estimation combined with the use of these sequences in a test receiver is called blind processing. The paper discusses the used algorithms, limitations, the expected performance in the chip error rate (CER), and effective loss of signal power when tracking the signals of opportunity (SoO) in a test receiver. An experimental setup with an array of 40 antenna elements is described, and results from simulated data and from one real global positioning system (GPS) M-code signal used as the signals of opportunity (SoO) show the feasibility of this concept. Among the types of global navigation satellite system (GNSS) signals of opportunity (SoO), the GPS M-code is more difficult to estimate than its Galileo or BeiDou counterparts due to its high chipping rate. A chip error rate (CER) of 15.1 % is achieved for the M-code signal. Applications of blind processing include receiver prototyping, signal quality monitoring of the signals of opportunity (SoO), and server-side processing for the purpose of signal authentication.
1 INTRODUCTION
In blind signal processing (BSP) problems, we do not have access (or have only partial access) to the structure of the transmitted signal (Cichocki & Amari, 2002). In satellite navigation, the received signal power levels are typically well below the noise floor; thus, a receiver must either know the structure of the transmitted signal or apply BSP. The most popular use of BSP is to track encrypted signals by exploiting similarities in the navigation signal structures in order to estimate signal parameters that are relevant for positioning, such as the delay, Doppler, or carrier phase (Betz & Cerruti, 2020; O’Driscoll & Curran, 2016; Woo, 1999). Other applications include signal analysis by timely deconvolution of spreading code sequences to raise the signal above the noise floor (York et al., 2014). A very direct method of blind processing is the use of a high-gain dish antenna, which directly raises the signal above the noise floor and thus reveals the signal structure. Setups with dish antennas are often used for signal quality monitoring (SQM), such as that reported by Thoelert et al. (2018 ,2019). Antenna arrays represent a potential alternative to dish antennas, potentially providing a high gain with the advantage that signals from different directions can be processed simultaneously. Antenna arrays are normally used as a countermeasure to jamming or spoofing (Ward, 1994), but antenna arrays have also been used for SQM, as described by Gunawardena et al. (2019). Additionally, dish antennas can obviously be used to retrieve unknown pseudorandom noise (PRN) code sequences, which has been discussed in detail by van der Merwe, Bartl, et al. (2020). Another use of BSP is to distinguish genuine signals from spoofing signals, possibly making use of artificial intelligence methods (van der Merwe, Nikolikj, et al., 2020).
If an antenna array is fully calibrated, i.e., if all group delay and phase delay values are known with high precision, it is generally straightforward to coherently superimpose the signals from the different elements while compensating for geometric delays corresponding to the direction of arrival (DOA) of the signal under consideration. Under ideal conditions, such as additive white Gaussian noise (AWGN) and negligible coupling between the antenna elements, the gain resulting from this kind of geometric (or nonadaptive) array processing is 10 log M, where M is the number of the antenna elements. Unfortunately, the calibration of such an array may require significant hardware efforts, and typically, M ≪ 10 in satellite navigation.
An alternative is to calibrate the antenna array on the fly, which can also be seen as a type of BSP. On-the-fly calibration of an antenna array can be achieved if a signal transmitted at the same (or similar) frequency with a known signal structure is received together with the signals of opportunity (SoO). The channel parameters (i.e., group and phase delay) for the individual antenna plus front-end paths are then retrieved from the known signal. After the parameters have been estimated, they can be applied to the signals of opportunity (SoO), thereby calibrating the array. This method is inherently robust and imposes almost no synchronization requirements on the used hardware, which can be consequently built with standard components. This design idea has been brought to hardware and software by the authors and forms the basis for this publication. A block diagram of the prototype system is shown in Figure 1.
The outline of this paper is as follows. Section 2 describes the method of calibration and the retrieval of pseudorandom noise (PRN) sequences and discusses the achievable chip error rate (CER) as well as important limitations when applying this approach to real-world signals. Furthermore, a link budget analysis is conducted to gain insight into the number of antenna elements required to be able to retrieve global positioning system (GPS) M-code chips. Section 3 covers the antenna array, front-ends, and software implementation on a many-core central processing unit (CPU) for chip retrieval, plus a description of how the sequences are applied within a test receiver. In Section 4, the software is first verified with simulated signals, and then, the entire demonstrator is applied to a real-world GPS C/A + M-code signal. Finally, the concluding section summarizes the lessons learnt and indicates future research directions as well as possible applications of this technique.
2 THEORY
2.1 Blind Chip Sequence Estimation
We consider a generic signal model for two code division multiple access signals arriving at one receiving antenna. One signal shall be called the “open” signal, and the other signal is denoted the “blind signal.” For each component, we assume a standard global navigation satellite system (GNSS) signal model (cf., Equation (8.19) of the well-known global navigation satellite system (GNSS) textbook of Kaplan & Hegarty (2017)) and superimpose the signals onto each other. Both signals are received together with AWGN background. Interference and contributions from other signals are neglected at this point. We assume that both signals have the same nominal center frequency. The possible presence of a data message is also ignored, as it is irrelevant for the following discussion. After a down-conversion to the baseband, the received signal model reads as follows:
1
where:
·OS, ·B are open and blind labels
s(t) is the received signal model expressed in the receiver time scale t
C is the signal power
fD is the Doppler frequency
ϕ is the carrier phase
k is the chip index
cOS,k, CB, k are PRN code chips (open, blind)
m is the modulation waveform
Tc is the chip duration
τ is the propagation delay
n(t) is the noise
We assume that the parameters mOS, Tc,OS, and cOS,k are known for the open signal. For the blind signal, we assume that only the modulation waveform mB and chipping rate Tc,B are known. These parameters can be retrieved either from interface documents or by measuring the power spectral density (PSD) of the blind signal with a high-gain antenna. We assume that mB and Tc,B are fixed parameters and do not vary in time. In contrast, the blind PRN sequence cB,k is unknown and can have an infinite duration.
For our proposed processing strategy, it is important that a relationship between the blind parameters fD;B,ϕB, and τB and the open parameters fD;OS, ϕOS, and τOS can be established. We distinguish two cases:
Identical transmitters: fD;B = fD;OS; ϕB = ϕOS + ϕoff; τB = τos
Two different transmitters: fD;B = fD; B(fD;OS); ϕB = ϕB(ϕOS); τB = τB (τOS)
The first case, with identical transmitters, is easily explained. Both received signals will exhibit the same Doppler shift and propagation delay. We only allow for a slowly varying phase difference ϕoff that may occur because of slightly different delays in the transmitting or receiving hardware chain or because of different antenna phase centers for the open and blind signal. Furthermore, we assume that both modulation waveforms are aligned at the transmitter side, meaning that, e.g., the open and blind chip edges coincide.
The second case is more generic and requires that fD;B, ϕB, and τB can be calculated from fD;OS, ϕOS, and τOS. This case arises if the positions of both transmitters and the receiver antenna are known. Then, one can infer the receiver clock offset and receiver clock drift from the open signal parameters fD;OS, ϕOS, and τOS. Once the receiver clock offset and receiver clock drift are determined, the blind parameters fD;B, ϕB, and τB can be calculated. This is similar to backward and forward modeling of code pseudorange, Doppler, and carrier phase observation equations.
Notably, the assumption of identical center frequencies is relevant for our discussion, as only in this case can hardware delays be assumed as identical for both signals. For the case of different center frequencies, further functional dependencies need to be established. It should be emphasized that different carrier frequencies might imply different local oscillators in the front-end, which would introduce further unknown phase offsets. Thus, hereafter, we will only consider the case of identical transmitters, with the blind and open signals having the same center frequency.
Considering that the open signal is tracked, the parameters fD;OS, ϕOS, and τOS can be known with high accuracy. Knowing τOS allows us to construct the satellite transmit time scale tSV:
2
Multiplying Equation (1) with the complex conjugate of the open signal phase gives the following:
3
where ñ represents the noise after multiplication and is the contribution from the open signal. Both terms will be discussed later. Equation (3) realizes a phase-aligned and transmit-time-aligned representation of the blind signal. Thus, the open signal allows one to construct, e.g., a “side channel attack” to the blind signal, bringing it to a form that allows an easy retrieval of the PRN code sequences.
With M antennas/receivers tracking the same satellite, their samples are coherently summed in order to obtain a good chip estimate for cB,k. We denote the signal model for the m-th antenna/receiver chain as and define the sum signal as follows:
4
It should be noted that Equation (4) does not contain any array factors, as the open signal phase of Equation (3) implicitly assumes the role of the array factor, i.e., the open signal phase factor depends on the relative position of the antenna elements.
If mB(t) and Tc,B are known (which we assume here), an optimal matched filter can be applied to demodulate the code values cB,k from . Assuming that the noise in ñ is AWGN and that mB and mOS are orthogonal, i.e.:
5
and that the chip values assume only the values +1 and −1, then the soft-decision estimated chip ĉs,B,k value is given by the following:
6
The hard-decision estimate ĉB,k is obtained by taking the sign of the real part, i.e.:
7
We interpret the hard-decision estimates as the conventional binary chip values represented by the sign, and the soft-decision estimates correspond to the real part of the complex-valued signal amplitude, thus indicating the reliability of those estimates. The blind sequences ĉB,k can then be stored in a file together with the respective time stamps for further processing of the blind signals within the test receiver.
When combining signals from M antenna elements, one must ensure that the polarity of the blind code chips is identical for all elements. This requirement is clearly fulfilled if the pilot component of an open signal is used to calibrate the array. If a data-only signal, such as the GPS L1 C/A code, is used to calibrate the array, the preamble must be decoded in order to resolve the polarity of the open signal, and after preamble decoding, the open signal phase-locked loop (PLL) needs to remain in phase lock. In this context, the occurrence of half-cycle slips is important, and our processing scheme continuously checks the polarity of the preamble and reports a half-cycle slip if a polarity change is detected.
2.2 Chip Error Rate and Correlation Loss
Under the assumptions of the last section, it is straightforward to calculate the likelihood that the estimated blind chip value is correct. The inverse of this likelihood is the chip error rate (CER). Assuming uncorrelated AWGN with an identical PSD N0 for all antenna/receiver paths, the noise power in equals MN0. The effective total blind signal power in is obtained by adding the individual blind signal amplitudes, i.e.:
8
The CER is given according to Equation (8.108) of (Kaplan & Hegarty, 2017) by the following:
9
with:
10
Under the assumption of equal signal power in each antenna/receiver path (CB;m = CB;1), we have the following:
11
and
12
When the blind sequences are applied within the test receiver for acquisition and tracking, the number of incorrectly estimated chips will reduce the correlation function of the received signal with the blind replica signal. This degradation is proportional to the CER. The effective power loss L is calculated as follows:
13
The interested reader is directed to the paper by Betz & Cerruti (2020), in which Section 3.2 “Model for dual-channel codeless processing” shows many similarities to the presented approach in terms of calculating L.
We conclude that the use of open signals to obtain synchronization and coherency avoids the use of a complex, cost-intensive coherent phased array system and allows the use of unsynchronized front-ends. Furthermore, there is no need for system calibration. The open signals perform an online calibration and therefore ensure an optimal performance. With the help of the open signal tracking, a complete phase-coherent synchronization of the antenna array is achieved. This synchronization applies to any point at which the open signal can be received; thus, the antennas can be located anywhere within the line-of-sight to the same transmitter. This includes both local antenna arrays and more spatially distributed arrays. The method can be applied in real time and during postprocessing.
Following this approach, the antenna+front-end chains have no stringent requirements regarding the sample stream synchronization. Synchronization errors that occur for our approach arise from code tracking noise of the open signal. This noise is typically very small compared with the chip length, even when the open signal, which is used for tracking, has a lower bandwidth and code rate than the blind signal.
For a proper constructive accumulation of signal streams in Equation (8), the carrier of the open signal must be wiped off for all signal streams. Having a sufficiently high carrier-to-noise ratio, which gives only a few degrees of phase error, and a low PLL bandwidth will allow one to constructively add the single signal streams into a combined stream.
2.3 Exemplary CER Assessment for GPS M-code
The feasibility of the chosen approach for estimating unknown PRN code chips is illustrated with an exemplary calculation for typical low-cost antenna/front-end parameters and the GPS M-code signal transmitted from GPS Block IIR-M satellites. Due to the high chipping rate, the GPS M-code signal represents a benchmark for this method; other encrypted GNSS signals or SoOs typically have lower chipping rates and are thus easier to estimate. Furthermore, it should be emphasized that the presented method actually estimates the BOC(10,5) component of the GPS signal, which is assumed to coincide with the M-code. The extent to which this assumption is true is unknown to the authors.
The following assumptions were made for the receiver antennas and front-ends:
Receive antenna gain pattern from approximately −5dBic at horizon to +5dBic at zenith
Noise figure of 1 dB for the antenna low noise amplifier (LNA)
Loss of 1 dB before the antenna LNA
Gain of 28 dB for the antenna LNA
Noise figure of 2.5 dB for the front-end
Antenna temperature of 140 K
Receiver noise temperature of 170 K, depending on quantization; the noise PSD was finally set between −203.0 and −203.7 dBW/Hz (2- and 4-bit quantization)
The losses included the following:
Atmospheric losses between 0.18 and 1.65 dB
Implementation losses of 0.2 dB (4-bit) and 0.69 dB (2-bit)
We assume the following at the transmission side:
Effective transmission power, including an antenna gain between 28.6 dBW and 30.9 dBW total equivalent isotropic radiated power within an off-boresight angle of 14°
Satellite antenna de-pointing loss of 0.25 dB
Further losses due to imperfect open signal codes and phase tracking are not considered. Phase errors would lead to amplitude losses because samples are not perfectly rotated in-phase in Equation (4). The errors grow toward lower elevations with lower carrier-to-noise ratios. Code tracking errors would lead to incorrect mapping of samples to neighbor chips, which, in the case of a different chip sign or binary offset carrier (BOC) modulation, leads to losses in both neighboring chips/BOC subchips, i.e., to errors in the superposition of the signals in Equation (4) or in the integration boundaries of Equation (7). This loss is assumed to be very small, because the code synchronization error of the open signals will be a small fraction of the chip/subchip size. To further reduce the synchronization error, carrier-aided code tracking can be applied. A further unknown loss could be possible coupling effects of the antenna elements, including correlated noise of the antennas, which is ignored here.
Figure 2 shows the required array gain in decibels (= 10 log M) to achieve a CER of 16 %, i.e., L = −3.35 dB. It can be seen that with 40 antenna elements (gain of 16 dB), it is theoretically possible to exceed the CER of 16 % down to an elevation angle of approximately 33°−40°. At the horizon, the received power decreases due to the satellite and receiver antenna patterns (resulting in an increasing required number of antenna elements). Furthermore, mutual shadowing of the antenna elements will further reduce the received power, unless the antennas are placed at a sufficient distance.
2.4 Limitations of the Approach
This section discusses a few constraints of the methodology, going beyond the previously discussed hardware limitations that will eventually increase the CER. The following points are discussed:
Estimation of the phase offset ϕoff
Interference from another blind signal source of the same type
Non-orthogonal modulation waveforms
The open signals of satellites/transmitters of interest can be acquired and tracked by a central processor for all antenna/front-end channels of an antenna array platform. In this way, a synchronization between the signal streams of all antenna channels is established. The code numerically controlled oscillator (NCO)s deliver the necessary timing information to synchronize the streams regarding the chip edges. The carrier NCOs individually synchronize their carrier to the common open satellite signal. In this way, all dynamics of the carrier can be wiped off from the signal. The final step is to account for the phase difference ϕoff between the open signal and the blind signal in Equation (6). This step can be performed by applying a fixed known phase offset or by estimating the phase offset. One efficient option for estimating ϕoff is to apply the noncoherent maximum likelihood (ML) phase estimator derived in (Borio & Lachapelle,2009) to a batch of the complex chip amplitudes (Equation (6)) as follows:
14
For the software receiver and the results within this work, the following algorithm was used:
15
Because the dynamics of the phase offset are very slow, the number of chips used for calibration can be very large, and the choice of the used algorithm does not influence the phase estimation or chip estimation. The length of the averaging window in Equation (15) and the rate of estimation depend on the technical circumstances resulting in the variable phase offset. If, for example, the variability is caused by two spatially separated transmission antennas (one for the open and one for the blind signal), then the parameters can be calculated by geometrical considerations. Estimation errors in will increase the CER and depend on the total received power CB and the length of the averaging window in Equation (15). If the window can be made sufficiently long, the errors can be neglected. According to Thoelert et al. (2019), the specific reason for the presence of a time-varying phase offset for GPS III satellites is the use of a separate amplifier and antenna chain. The necessity of offset estimation is evident, as will be shown later in Figure 13.
Thus far, only signals from a single transmitter have been considered. Of course, this method can also be applied to scenarios in which multiple transmitters broadcast open and blind signals of the same type. We assume that the open signals can be well-distinguished from each other, as in GNSS. In this case, the methodology provides two levels for focusing on the blind transmitter corresponding to the considered open signal. First, the integration boundaries in Equation (6) are defined by the time scale of the open signal under consideration, which are assumed to be aligned with the corresponding blind signal. The chip boundaries will generally be different for the other blind signal sources, and thus, the contributions of the other blind signals have a tendency to average out in Equation (6). However, in the case of two blind signals with the same Doppler shift, situations may occur in which the blind chip of the other signal falls for a longer time span into the time slot of the considered blind signal. In that case, the estimation process will be distorted. The second level for mitigating the influence of other blind signals on the chip estimate is the DOA. Equation (3) and Equation (4) represent a beam-forming maximized array gain toward the DOA of the open signal (and thus toward the DOA of the considered blind signal). Other blind signals with a similar DOA will negatively impact the CER; however, if the DOA difference is large, the impact can be neglected. This is very similar to the case of a dish antenna observing signals from two close satellites; in this case, even a dish antenna cannot separate the signals from each other. Both levels of focus (integration boundaries and DOA) are, in principle, independent of each other and should ensure that blind signals from different sources can be well discerned from each other. However, it should be mentioned it will be difficult to discriminate two blind sources moving along nearby trajectories.
Considering that there are at least two signals present, one open and one blind, there are different possibilities regarding the complexity for retrieving the blind chips without interference from the open component or any further components. The more simple cases occur when the blind signal is spectrally separated or waveform-orthogonal to the open signal, as indicated by Equation (5). A more difficult case occurs if the signals overlap in the spectrum or are non-orthogonal. One example for this is the BeiDou B1 signal, which has an overlap of the legacy B1I (as an open signal) with the lower sideband of B1A BOC(14,2) (taken as the blind signal). To avoid a spillover from the open signal to the blind PRN sequence estimate, an orthogonalization procedure can be applied to obtain a modified blind modulation waveform :
16
It should be noted that this equation accounts for the possibility that the modulation waveforms are complex-valued, which is the case if the waveform also represents a frequency offset. The orthogonalization procedure causes a reduced signal power for the blind chip estimation process. For the case of B1I and B1A, the orthogonalization process is very similar to single-sideband filtering, because the signals are very similar in the lower sideband. It should be mentioned that further signal components present in the navigation signal may require orthogonalization of the blind waveform. Finally, components sharing the same modulation, e.g., the BOC(1,1) component of Galileo E1B and E1C, cannot be separated at all.
3 DEMONSTRATOR SYSTEM
To demonstrate the feasibility of this approach, a prototype system making extensive use of GNSS software radio technology was set up, as described in this section. The system uses low-cost GNSS antennas and front-ends, collects raw intermediate frequency (IF) samples on a server, and eventually runs a multi-threaded GNSS software-defined radio (SDR) implementing the estimation algorithms. The estimated blind chips are then fed into a commercially tailored GNSS SDR.
3.1 Hardware Setup for Blind Processing
The antenna array consists of 40 Tallysman VSE6137 antennas with VeroStar Technology by Tallysman Wireless, Inc., as shown in Figure 3. This type of antenna allows the reception of GNSS signals with a homogeneous gain from horizon elevation up to zenith. The antennas have inbuilt LNAs with a gain of 37 dB. By omitting housings and placing the antennas on a common ground plate covered by a single radome, the antennas are aligned in a dense arrangement at a separation distance of only 106 mm. The whole array can be tilted up to 45° to allow the reception of satellites near the horizon without shading from neighboring antennas.
Data acquisition is performed by 20 low-cost customer off-the-shelf (COTS) SDRs of the type bladeRF 2.0 micro xA4 by NUAND, LLC, allowing a sampling rate of 61.44 MS/s in-phase quadrature (IQ) data for a single channel and up to 35 MS/s and a resolution up to 12 bits for 2×2 multiple-input/multiple-output streaming over a USB3.0 connection, each collecting data from two antennas (cf., Figure 4 left). The core is the latest-generation Cyclone V field programmable gate array (FPGA) by Intel. All SDRs are synchronized by a GPS disciplined oscillator (GPSDO) in a daisy chain.
The data stream from each SDR consists of intermittent IQ data of two antennas with a resolution of 12 bits. These data are reduced in real time to the selected resolution of 8, 4, or 2 bits before being stored on an Intel next unit of computing (NUC) PC with a Core i3 processor, equipped with a solid state drive (SSD) with 1-Tb capacity, thus allowing 2 hours of recording (cf., Figure 4 right). More details of the system can be found in (Lesjak et al., 2022).
The platform of the test receiver for blind signals is equipped with similar hardware for data acquisition, but different COTS antennas. Two Intel NUC PCs record the data from two bladeRF SDRs connected to four Tallysman TW7972 triple-band GNSS antennas. For comparison purposes, a low-cost mass-market COTS receiver (μBlox M8T) connected to a Tallysman TW2710 antenna is attached to the common carrier plate as well, which by itself can be mounted on a tripod or on a car’s roof rack. Both SDRs and recording NUC PCs are housed in a compact, lightweight, easy-to-transport 19″ flight case.
After acquisition, the data are copied via 10-Gb data links to a blind GNSS processing module (BGPM) server PC (Figure 5 right) equipped with a 32-Tb redundant array of inexpensive disks consisting of SSDs for blind processing. A scheme of the overall setup is displayed in Figure 5 left. Because each of the data streams contains intermittent data from two separate antennas, the data must be split to separate files before processing. It is not possible to split and write out the streams to separate files during acquisition, because both serial advanced technology attachment (SATA), and USB interfaces share the same bus and would not allow the necessary bandwidth.
The data acquisition planning (selection of data rate, recording resolution, and processing profiles for chip estimation and use within the test receiver, downloading of assisting data, and timed start of acquisition) is done within a single graphical user interface. The user can select manual step-by-step processing or fully automatic processing in batch mode, starting at any point in the sequence.
3.2 Software Considerations for Blind Processing
Processing raw IF signals from 40 antennas is a computationally demanding task for a GNSS SDR. Figure 6 presents an illustration of the chip estimation algorithm and data flow. According to their synchronization points (e.g., NCO-based code phase), the streams are added to one final stream accounting for the carrier phase, as indicated in Equation (4). In the case of chip estimation, the stream is further compressed into one complex-valued amplitude per chip through correlation along each chip with the known BOC or any other modulation, cf., Equation (6). Finally the binary chip estimates (Equation (7)) are stored in a file together with the transmit time epochs. The synchronization and chip estimation are illustrated in Figure 6 and work as follows:
All antenna streams are tracked by delay-locked loop (DLL) and PLL individually for each satellite.
During the replica generation for the open signal, replicas for the blind signal are also generated. Based on the open signal code NCO, this gives a chip index (per signal sample) of the blind signal, as well as a replica (per sample) of the assumed blind modulation waveform mB. The carrier replica of the open signal tracking can be reused.
In the correlation process for open signal tracking, the carrier replica is applied, including a constant phase shift in order to wipe off the carrier and rotate to the blind signal phase.
Additionally, the assumed modulation waveform of the blind signal is multiplied.
For each sample, the resulting chip estimate of the sample is added to the chip stream of the specific antenna at the chip index calculated in step 2.
After all receiver channels have indicated a correlation dump of the open signal tracking, all antenna chip streams are summed per sample and written to the output file.
All IF streams are processed in parallel. Because of inevitable differences in the start times of the streams, the algorithm foresees a certain ring buffer that needs to account for the maximum start time differences (up to 1 s in the worst case). Related to this comparable large time difference, we manually assessed the time differences between the network-time-protocol-synced NUCs and observed a maximum difference between individual NUCs of less than 10 ms. Thus, we conclude that the time difference was caused by variable software latency to launch the sample recording. It would be possible to share a trigger by using a special FPGA image and chaining multiple BladeRFs via a mini expansion header (cf., https://github.com/Nuand/bladeRF/issues/221), but this has not been implemented in our system.
The estimation algorithm was implemented in the multi-sensor navigation analysis tool (MuSNAT) described by Pany et al. (2019). The development of the core of MuSNAT goes back to times of single-core CPUs without real parallelism for processing different signals of different sources and frequencies and their respective components. Nowadays, multi-core and many-core CPUs are fundamental for energy efficiency and for addressing the increased processing demands of applications. To exploit this feature for multi-antenna processing, the MuSNAT was profiled and optimized for the Threadripper architecture, featuring up to 64 cores (Advanced Micro Devices (AMD), 2022).
Figure 7 shows the timing diagram of the MuSNAT running on an AMD Ryzen Threadripper 3970X with 32 cores, tracking 40 C/A code signals as open signals with an 80-MHz sampling rate from 40 different signal streams. Within the upper part of this figure, the most time-consuming functions are listed, which are “generateCodeDP” for baseband replica signal generation and “generateCarrier” for carrier replica generation. The Open Multi-Processing (OMP) library was used to parallelize the processing of different PRNs. The main part of the figure indicates the well-distributed work load of the different OMP threads as continuous brown lines. The integrate and dump epochs for the GPS C/A code signals are marked as color ticks within the thread diagrams. The time axes shows nearly 300 ms of computational time. In contrast to the behavior for only a few streams, in additional to the primary work of replica generation and correlation, the work of stream reading from the file and processing become more dominant and are parallelized.
MuSNAT processes the IF streams on a packet-per-packet basis. The packet length ranges from a few milliseconds to a few hundreds of milliseconds and has a substantial impact on the processing performance. Figure 8 compares the effective CPU utilization for two different stream packet sizes. The lower histogram shows that for longer packets, the CPU moves into a state in which all 40 channels are working in parallel, and sometimes for file reading and sample bit format conversion with additional threads, more than 60 logical cores are used (more than 38 on average). For a smaller packet (upper plot), only 26 cores are used on average. Tracking several satellites per antenna stream or frequency band will decrease the thread waiting times because of the additional work that must be performed per packet and will also allow efficient use of the full CPU capabilities.
Overall, the algorithm runs at a comparable speed on the Threadripper CPU. Sole tracking of two GPS C/A code signals for each of the 40 antenna elements under the settings described in Section 4.2 runs within approximately 93 % of real time. If the blind chips are also estimated, the processing time is approximately 480 % of real time.
3.3 Considerations in the Test Receiver
An IGASPIN SX3 software receiver was used as the test GNSS receiver to acquire and track the retrieved blind sequence from the BSP. Here, three different approaches were implemented and tested for acquiring the blind sequence. In the first approach, the blind sequence is acquired by using a fast Fourier transform (FFT) with the help of a blind assistant data file containing information about the traveling time (or transmission and reception time) of the related signal. In the second approach, the blind signals are acquired by utilizing the results from the navigation processor. In this case, no FFT acquisition is required, and the true traveling time of the signals is directly calculated. Thus, the software receiver must use a hot start. In the third approach, the open signal of the same satellite must be acquired and tracked by the receiver beforehand, and then, the required information can be obtained.
As soon as a blind signal is acquired, the acquisition information is passed forward to the tracking unit of the receiver. In this part, a new implemented correlator called a “blind correlator” correlates the infinite retrieved blind sequence with the input signal. The raw measurements are then fed to the navigation processor for positioning, if the number of measurements is sufficient.
4 RESULTS
Simulated GPS C/A and M-code-like signals, i.e., a BOC(10,5) signal and real-world GPS signals, were used to test the demonstrator and to prove the feasibility of this approach.
4.1 Chip Sequence Estimation for Simulated Data
To verify the software described in Section 3.2 for estimating the blind PRN code sequence, simulated GPS C/A+BOC(10,5) signals were generated and processed. As the values for BOC(10,5) PRN code sequences are known, the CER can be calculated directly. A MATLAB-based tool was used to generate 40 IF sample streams, each corresponding to one antenna element. Only one satellite was simulated. The signal structure contains the C/A code component, including the navigation message plus a BOC(10,5) pilot signal representing the M-code. The spreading code for the BOC(10,5) signal was chosen to be periodic, with a period of 1 ms. The signal power for each component, C/A and BOC(10,5), was defined to be 50 dBHz, as measured after filtering and before analog/digital conversion (ADC). The phase offset ϕoff was 90°. The sample rate was 80 MHz, and IF 20-MHz and 2-bit sampling was simulated, corresponding to a quantization loss of 0.6 dB. It should be noted that each of the IF sample streams has a slightly different starting time (within ± 100 ms) to test the synchronization capability of the software.
In the first run, the phase offset ϕoff was assumed to be known. The estimated chips were compared against the true values, and the CER was compared with the theoretical value. In the first evaluation, the blind chips were estimated from only one of the 40 simulated IF sample streams. The CER is computed by comparing batches of 10230 chips to the true values and is shown in Figure 9 for a single antenna/receiver path. If Equation (12) is evaluated for this case, we need to account for Tc;B = 0.196 μ s, CB;1/N0 = 49.4 dBHz, and M = 1. This results in a CER of 42.7 %, in good correspondence with Figure 9. The I/Q histogram of the soft-decision chip estimates shown in Figure 10 exhibits a nearly symmetric structure corresponding to the high CER and shows a pattern related to the 2-bit quantization of the IF sample stream. An I/Q histogram shows the likelihood of occurrence of the complex-valued soft-decision chip estimates (Equation (6)) as a color figure (blue … value never occurs, red … value frequently occurs; I … real part, Q … imaginary part).
If all 40 streams are processed in a combined manner, M = 40 and the CER decreases to 12.2 %. The CER obtained from processing the simulated signals is shown in Figure 11, and the I/Q histogram of Figure 12 clearly shows an elongated structure. The obtained CER is slightly higher (i.e., worse) (≈4 %) than the ideal. We treat this as a blind chip estimation implementation loss of 1.3 dB, which might be explained by DLL jitter that causes minor misalignment between the streams or other quantization effects.
In the second run, the ability to estimate the phase offset ϕoff was tested. The offset was estimated each full second by evaluating the sum of Equation (15) over an interval of 100 ms. The exactly same CER was finally obtained as listed in Table 1. The exact coincidence is attributed to the fact that the search in Equation (15) over all phase offsets is quantized to steps of 1°. Obviously, the BOC(10,5) signal was sufficiently strong for the applied estimation interval, with the correct offset of 90° being identified.
4.2 Chip Sequence Estimation for Real Data
To further verify the method, numerous tests with live GNSS signals were conducted. The antenna array in Section 3.1 was used to record signals, and then, PRN sequences of the blind signal were extracted using the software described in Section 3.2. To have a ground truth, a dish antenna with a diameter of 2.4 m was pointed to the satellite of interest, and the corresponding IF signal was recorded in parallel. The blind processing described in Section 3.2 was then applied to the dish antenna stream. Because of the high gain of ≈ 30 dB of the dish antenna, the resulting CER can be considered to be negligible. Finally, the CER is obtained by comparing the dish antenna PRN sequences to the sequences from the antenna array.
In the following, one exemplary result is presented for the GPS satellite PRN18. This Block III satellite is one of the few that seem to broadcast the M-code over Europe (Steigenberger, 2020; Steigenberger et al., 2018). Tests with other GPS satellites and other systems have been performed. The results presented in the following are representative of the other cases.
The antenna array measurement was conducted in Graz, Austria on April 7, 2022. Each stream was stored in an 8-bit I/Q format with a sample rate of 35 MHz. The streams are at baseband, i.e., IF = 0 MHz. The recording lasted approximately 30 min, but only 20 s of the recording is presented here as a representative result. Standard GPS C/A code receiver processing was applied to each stream individually to assess the received GPS C/A code power. Summing all power estimates according to Equation (8) gives 61.68 dBHz for GPS PRN18. The lowest GPS C/A code power for a single antenna was 38.27 dBHz, whereas the highest power was 49.49 dBHz. The high variability of signal power for the different antenna elements requires further investigation and is beyond the scope of this paper. The sampling start was synchronized via software, and a maximum start time difference of 0.8 s was observed between the streams. The platform is, however, able to provide millisecond-accurate synchronization. The dish antenna was located in Munich, Germany.
The soft-decision I/Q plot for the dish antenna shown in Figure 13 demonstrates that the CER is near zero. The two binary chip values of the BOC(10,5) component broadcast by this satellite can be easily discriminated. Phase offset estimation via Equation (15) was not employed, as the phase misalignment was compatible with the good separation of both chip values. Thus, the data from the dish antenna can serve as a reference.
An I/Q plot for the BOC(10,5) soft-decision chip estimates obtained for PRN18 via the antenna array is shown in Figure 14. During processing, the phase offset between GPS C/A and the BOC(10,5) component is estimated each second by evaluating the sum of Equation (15) over an interval of 100 ms. The plot demonstrates that phase offset estimation obviously aligned the elongated side parallel to the I-axis.
The CER is shown in Figure 15. The mean value is 15.1 %. As the M-code transmit power is unknown for this satellite, one cannot compare the CER against any known true value. A rough estimate can be obtained by assuming that the combined received power of the BOC(10,5) component equals the C/A code power of 61.68 dBHz, which, according to Equation (12), would result in a CER of 22.4 %. An additional 2.6 dB would be required to explain the observed lower CER. In parallel, the GPS PRN23 (another Block III) satellite was processed with the antenna array, but not with the dish antenna. For PRN23, we noted a less-elongated structure, indicating a lower transmit power for that satellite. Other M-code-transmitting satellites were not available during the measurement interval.
Finally, the PRN sequence obtained from the antenna array was used to generate a BOC(10,5) replica signal, and this replica signal was correlated against a received GNSS signal from an antenna of the test receiver. The replica signal was adjusted in carrier phase, frequency, and Doppler phase to match the corresponding parameters of a tracked PRN18 GPS C/A code signal within the IF sample stream of the test receiver antenna. The correlation time was 100 ms, and the resulting correlation function is shown in Figure 16. One can clearly observe the BOC(10,5) auto-correlation function, indicating that the estimated blind PRN sequence can be found within another independently captured GNSS signal. The observed correlation amplitude corresponds to a C/N0 value of approximately 48 dBHz. We also observe that the peak occurs at an offset of 400 ns compared with the GPS C/A code, but this is explained by an inconsistent chip indexing scheme within the MATLAB script used to compute the correlation function. In contrast, the test receiver described in Section 4.3 does not show a 400 ns offset between pseudoranges of the C/A-code and the BOC(10,5) signal (see Figure 20).
4.3 Tracking Blind Signals Using an Estimated Chip Sequence
Estimated blind sequences of the BOC(10,5) component for the GPS PRN18 satellite were used to acquire and track the respective component within the test receiver, as described in Section 3.3. It is emphasized that the single antenna of the test receiver is not used for sequence estimation. The presented data set was collected at Graz, Austria on April 25, 2022 (different date than in the previous subsection). The signal from the one GNSS antenna was sampled at 35 MS/s (I/Q) and was converted within the SX3 to a real-valued IF sample stream with 70 MS/s. To track the BOC(10,5) component, a standard early-late tracking scheme with a bump-jumper was configured in the test receiver using an early-late correlator spacing of 8 IF samples, corresponding to 114.29 ns or 0.58 chips (assuming a chipping rate of 5,115,000 chips/s), a second-order DLL with a bandwidth of 1 Hz, and a second-order PLL with a bandwidth of 9 Hz. The test receiver was static in open-sky conditions. A coherent integration time of 4 ms was used.
Figure 17 shows the prompt correlator values as a function of time. The overall behavior shows clearly stable tracking. The blind signal is effectively a pilot signal, as the possible presence of any navigation message is absorbed into the chip estimates. One can also see a variation in signal power, which could stem from either variations in the received power, variations in the CER, or slowly varying local multipath effects. A first investigation to separate those effects did not show any clear result, e.g., the GPS C/A code C/N0 values in Figure 21 exhibit a different pattern. Further analysis, e.g., by investigating the soft-decision estimates (Equation (6)), is needed. Nevertheless, the discriminator values shown in Figure 18 further confirm stable tracking. The variations in the discriminator noise seem to correspond to the blind signal power variations. Figure 19 shows the difference between the resulting code pseudorange and carrier phase expressed in meters. This difference is a measure of the code multipath effects and the code thermal noise. The root mean square value for t = 300 – 320s is 25 cm, and during this time, the C/N0 is estimated to be approximately 45 dBHz by the SX3. It is noted that this value is higher than the theoretical expectation for a BOC(10,5) signal with AWGN Betz (2000), and the specific reasons for this deviation might stem from the high correlator spacing, an incorrect discriminator gain, or numerical approximations within the test receiver. These considerations are, however, beyond the scope of this paper, whose main focus is the chip estimation process. For the sake of comparison, Figure 21 and Figure 22 show related plots for the GPS C/A code signal received simultaneously from the same satellite. The same tracking settings were used, but the correlator spacing was reduced to 57.14 ns or 0.058 GPS C/A code chips. The variability is much higher in the BOC(10,5) code noise in Figure 19 than in Figure 22, indicating that different loop bandwidth values are used. Figure 20 shows the difference between the pseudoranges of the C/A and blind codes, which confirms that there is no bias between these two pseudoranges. Overall, the BOC(10,5) tracking results match our expectations for this early development stage of BOC(10,5) usage and are inline with the GPS C/A code data.
5 CONCLUSIONS
In this work, the ability to estimate PRN sequences by a self-calibrating antenna array has been theoretically described and verified with simulated and real-world data. Limitations of the approach have been pointed out. The main purpose of the demonstrator system is to support the development of GNSS receiver algorithms for higher-order BOC signals and to test these algorithms with real-world signals in an unclassified domain. The algorithms include acquisition, tracking, and, if a low CER is achievable, even SQM. The SQM use case would benefit from further research on indications regarding the extent to which the obtained PRN sequences are free from the limitations discussed in Section 2.4. Currently, degradation in the estimation process can only be observed on an empirical basis and is far from being repeatable. Another application of the array could be found in using the sequences for server-side signal authentication (Rügamer et al., 2016). By using the array, the sequences can be retrieved from a signal in space at a trusted, spoofing-free location, thereby avoiding the need for a classified server room.
AUTHOR CONTRIBUTIONS
Dominik Dötterböck implemented and optimized the blind sequence estimation and performed various theoretical analyses and simulations. Thomas Pany provided the system concept, including the mathematical description. Thomas Prechtl was responsible for the system architecture, system implementation, user software implementation, system testing, and data evaluation. Roman Lesjak coordinated the research project and contributed to the system design, system implementation, and testing. Amir Tabatabaei was responsible for the test GNSS software receiver section to acquire and track the blind sequence retrieved by BSP.
HOW TO CITE THIS ARTICLE
Dötterböck D., Pany T., Lesjak R., Prechtl R., Tabatabaei A. (2023). PRN sequence estimation with a self-calibrating 40-element antenna array. NAVIGATION, 70(4). https://doi.org/10.33012/navi.600
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.