Huber’s Non-Linearity for GNSS Interference Mitigation ^†

Abstract

Satellite-based navigation is prevalent in both commercial applications and critical infrastructures, providing precise position and time referencing. As a consequence, interference to such systems can have repercussions on a plethora of fields. Additionally, Privacy Preserving Devices (PPD)—jamming devices—are relatively inexpensive and easy to obtain, potentially denying the service in a wide geographical area. Current jamming mitigation technology is based on interference cancellation approaches, requiring the detection and estimation of the interference waveform. Recently, the Robust Interference Mitigation (RIM) framework was proposed, which leverages results in robust statistics by treating the jamming signal as an outlier. It has the advantage of rejecting jamming signals without detecting or estimating its waveform. In this paper, we extend the framework to situations where the jammer is sparse in some transformed domain other than the time domain. Additionally, we analyse the use of Huber’s non-linearity within RIM and derive its loss of efficiency. We compare its performance to state-of-the-art techniques and to other RIM solutions, with both synthetic and real signals, showing remarkable results.

1. Introduction

Global Navigation Satellite System (GNSS) is the technology of choice for most position- and timing-related applications when available [1]. GNSS relies on a constellation of satellites synchronously emitting known signals from which one can compute its position. GNSS encompasses the American GPS, the European Galileo, the Russian Glonass and the Chinese Beidou among others [2,3,4]. The growing dependence on GNSS within critical (and non-critical) infrastructures has posed some concerns about the potential vulnerabilities of GNSS [5]. GNSS signal reception is vulnerable to several forms of Radio Frequency (RF) interference [6] that can significantly degrade receiver performance. For this reason, several approaches have been developed to strengthen GNSS receivers against interference and jamming [5,7,8]. Jamming is a type of intentional interference where a powerful signal is broadcast in the GNSS frequency bands with the ultimate goal of preventing receiver operations. Jamming is a form of Denial of Service (DoS) [9] and can have a significant impact on critical infrastructures relying on GNSS as the location service provider [10].

In the presence of such DoS attacks, GNSS receiver performance can be significantly improved by implementing interference mitigation techniques [7]. These techniques either attempt to cancel the interfering signal or implement processing strategies that are only marginally affected by spurious signal components. Most interference mitigation techniques can be interpreted under the framework of Interference Cancellation (IC) at the receiver. At a glance, IC first detects the jamming signal and then reconstructs its waveform, which is subtracted from the received sample stream to clean the data. Popular IC methods are, for instance, pulse blanking and (adaptive) notch filtering [11,12]. This approach has the drawback of requiring interference detection and estimation, which are two potential points of failure in the chain [13].

In contrast, a novel approach was recently proposed based on the theory of robust statistics [13]. We refer to this approach as RIM, and it has the remarkable advantage of avoiding estimation of the interference waveform and its detection. RIM takes advantage of the sparse nature, in some domain, of interference, given that the desired GNSS signal is generally not sparse. Basically, received signal samples affected by interference can be treated as outliers (after a convenient transformation), and thus their contribution is rejected in a principled manner using robust statistics [14]. The concept was first implemented using a filter to process the jamming signals that were then seen as pulsed interference by the GNSS receiver [15]. The myriad non-linearity was derived starting from a Cauchy assumption on the input data corrupted by jamming. The work was extended in [16] using the complex signum non-linearity, where a Laplacian model was considered. Both works substituted the classical Gaussian assumption with a more relaxed assumption on the noise statistics that accounted for heavier tails in the noise probability density function (pdf), typically caused by large outliers in the data. In this paper, we explore the use of Huber’s non-linearity [17] as an alternative to specifying distributions, which is a popular methodology in the field of robust statistics. In practice, this choice entails a piecewise loss function that has higher outlier mitigation capabilities than standard methods based on the Gaussian assumption. Additionally, the RIM concept is extended to interference signals that are sparse in some domain other than the time one. For instance, a Continuous Wave (CW) is an outlier in the frequency domain, where only few samples are affected by the jamming signal. On the other hand, pulsed interference is sparse in the time domain. These and other situations can be formulated in a joint framework, which is one of the contributions of the present article. While this work focuses on time and frequency domain processing, other linear transforms can be used to obtain a sparse signal representation. For example, a Discrete Wavelet Transform (DWT) can be adopted and RIM can be interpreted as a generalization of wavelet excision algorithms [18,19]. The analysis of wavelet-based algorithms is outside the scope of this paper, and it is not considered here.

This work is based on the conference paper in reference [20]. With respect to reference [20], the paper provides several additional contributions. A theoretical analysis of the loss of efficiency caused by Huber’s non-linearity is provided. In particular, a theoretical formula capturing the impact of the decision threshold is derived. The proof of the properties used for the derivation of the loss of efficiency is entirely provided in Appendix A.

The performance of RIM with Huber’s non-linearity is assessed by simulations. Extensive simulations considering two types of jamming signals are provided, namely CW and swept signals. The synthetic simulations presented here are aimed at characterizing the robustness of RIM to different types of interference. This analysis is new and has not been presented before. In addition to this, new experimental results considering tests performed in a large anechoic chamber are presented. The analysis considers both acquisition and tracking stages. The findings obtained complement the results presented in reference [20] that considered narrowband GPS signals collected in a road environment and wideband Galileo E5b signals affected by jamming. Summing up, the results presented here complement both the simulations and the experimental analysis described in reference [20] and provide a complete theoretical analysis of the loss of efficiency entailed by Huber’s non-linearity.

The remainder of the paper is organized as follows. Section 2 presents the GNSS signal model considered in this work, the standard receiver processing based on the maximization of the Cross Ambiguity Function (CAF), and its robust counterpart. Then, Section 3 particularizes the robust CAF concept to the use of Huber’s non-linearity, while Section 4 provides an efficiency analysis of such approach. The proposed methodology is validated with both synthetic and real-world experiments in Section 5 and Section 7, respectively. Section 6 details the experimental configuration for the real data gathering. Finally, the paper concludes with final remarks in Section 8.

2. Signal and System Model

The signal collected by the antenna of a GNSS receiver can be modelled as [2]:

$$y(t)=\sqrt{2C}d\left(t-{\tau }_0\right)c\left(t-{\tau }_0\right)\cos \left(2\pi (f_{RF}+f_0)t+{\phi }_0\right)+\eta (t)+i(t),$$

(1)

which is the sum of three terms—the useful signal component, a noise and an interference term. In (1), C is the useful signal power, $d(·)$ is a navigation message and $c(·)$ is a pseudo-random code from a family of quasi-orthogonal sequences. The useful signal is delayed by the communication channel that also introduces a Doppler shift, $f_0$, with respect to the nominal signal RF, $f_{RF}$. ${\tau }_0$ and ${\phi }_0$ denote the delay and phase shift introduced by the communication channel. It is noted that a GNSS receiver collects several signals from several satellites. Due to the quasi-orthogonality of the codes, $c(·)$, used to broadcast the different signals, the receiver processes them in an independent way. For this reason, a single useful component is considered in (1).

The noise term, $\eta (t)$, is a zero-mean Additive White Gaussian Noise (AWGN). The term, $i(t)$, can assume several forms [21,22] depending on the type of interference considered. The receiver amplifies, filters, down-converts and samples $y(t)$, producing the baseband complex sequence

$$y[n]=\sqrt{C}\tilde{d}\left(n{T}_s-{\tau }_0\right)\tilde{c}\left(n{T}_s-{\tau }_0\right){\mathrm{e}}^{j2\pi f_{0}n{T}_s+j{\phi }_0}+{\eta }_{BB}[n]+i_{BB}[n],$$

(2)

where n is the time index, and square brackets are used to denote discrete time sequences sampled at a frequency of $f_s=\frac{1}{T_s}$. The subscript “BB” denotes filtered signals down-converted to the baseband. The symbol $\widetilde{·}$ indicates the impact of front-end filtering on the useful signal component.

Noise, ${\eta }_{BB}[n]$, is modelled as an AWGN with independent and identically distributed (i.i.d.) real and imaginary parts, each with variance, ${\sigma }^2$. A model commonly adopted for ${\sigma }^2$ is

$${\sigma }^2=N_0{B}_{Rx},$$

(3)

where $B_{Rx}$ is the front-end, one-sided bandwidth, and $N_0$ is the Power Spectral Density (PSD) of the input noise, $\eta (t)$. The ratio between the carrier power, C, and the noise PSD, $N_0$, defines the Carrier-to-Noise power spectral density ratio ($C/N_0$).

2.1. Standard Processing

The goal of the signal processing stage of a GNSS is to estimate the signal parameters, ${\tau }_0$, $f_0$ and ${\phi }_0$. In standard processing, signal parameters are determined by first computing the CAF, defined as [2]

$$C(\tau ,f_d)=\sum _{n=0}^{N-1}y[n]c\left(n{T}_s-\tau \right){\mathrm{e}}^{-j2\pi f_{d}n{T}_s},$$

(4)

where N is the number of samples used for the computation of the CAF that is integrated over $T_c=N{T}_s$, the coherent integration time. The CAF is a bi-dimensional function that depends on $\tau$ and $f_d$, the locally tested code delay and Doppler frequency. The signal parameters are then computed as

$$\begin{array}{cc}\hfill \left\{\widehat{\tau },{\widehat{f}}_d\right\}& =\arg \underset{\tau ,f_d}{\max }|C(\tau ,f_d)|\hfill \\ \hfill \widehat{\phi }& =\angle C(\widehat{\tau },{\widehat{f}}_d)\hfill \end{array}.$$

(5)

The maximization process described by (5) is implemented by acquisition and tracking algorithms that employ a bank of correlators to compute and maximize the CAF. The correlators are also used to estimate the $C/N_0$ that is thus determined post-correlation. The interference term, $i_{BB}[n]$, can significantly distort the CAF, preventing proper receiver operations. In this way, the estimated $C/N_0$ is strongly influenced by the interference term, $i_{BB}[n]$, and by the interference mitigation strategies implemented by the receiver. In this respect, the estimated $C/N_0$ can be used to quantify the effectiveness of interference mitigation approaches. Pre-correlation mitigation techniques, such as those based on Huber’s non-linearity, are applied before the computation of (4) on the digital samples, $y[n]$, which are, at first, pre-processed to reduce the impact of $i_{BB}[n]$.

2.2. Robust Interference Mitigation

RIM is a pre-correlation technique, i.e., applied directly to the samples, $y[n]$, before correlation, which can be implemented according to the three steps depicted in Figure 1 [20,23]. A linear tranform, ${\mathbf{T}}_1$, is, at first, used to project the jamming component into a domain where it admits a sparse representation, i.e., it affects a limited number of samples. For example, a low-pass filter can be used to make the jamming component appear as a sequence of time pulses [15,24]. This approach is effective for instantaneously narrowband jamming signals that sweep a large frequency band. When filtered, these signals result in a sequence of pulses. Transform Domain (TD) processing was suggested by [23], where the DFT was used to project the samples, $y[n]$, into the frequency domain where $i_{BB}[n]$ admitted a sparse representation. Transform ${\mathbf{T}}_1$ produces the TD samples

$$Y[k]={\mathbf{T}}_1(y[n]).$$

(6)

The change of index, from n to k, is a notational convention adopted to indicate that the input samples, $y[n]$, have been brought to a different representation domain. When ${\mathbf{T}}_1$ is the identity operator, the change of index used for the enumeration of the TD samples is no longer required.

Following ${\mathbf{T}}_1$, a Zero Memory Non-Linearities (ZMNL) is used to reduce the impact of outliers in the TD. The complex signum and myriad non-linearities were considered in reference [15,24]. In this paper, Huber’s non-linearity [17,25] is analysed. While this ZMNL is widely used in robust statistics, its application to complex samples for interference mitigation is new, and it is one of the main contributions of this paper. A generic non-linearity is denoted here as $\psi (·)$ and produces the samples

$$Y_{\psi }[k]=\psi (Y[k]).$$

(7)

Finally, a second linear transform, ${\mathbf{T}}_2$, is applied to the samples, $Y_{\psi }[k]$, to obtain new filtered time domain samples. ${\mathbf{T}}_2$ inverts the effects of ${\mathbf{T}}_1$ and brings back the samples from the TD to the time domain. In time domain processing, ${\mathbf{T}}_2$ is simply the identity. The output of ${\mathbf{T}}_2$ is denoted here as

$$\tilde{y}[n]={\mathbf{T}}_2(Y_{\psi }[k]).$$

(8)

As a consequence, two special cases of RIM can be identified:

• Time domain processing: when ${\mathbf{T}}_1$ is a low pass-filter and ${\mathbf{T}}_2$ is the identity operator. In narrowband receivers, low-pass filtering is not needed. and ${\mathbf{T}}_1$ can be replaced by the identity operator,
• TD processing: ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$ are inverse operators and

  $${\mathbf{T}}_1\circ {\mathbf{T}}_2=\mathbf{I},$$

  (9)

where $\mathbf{I}$ is the identity operator. If ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$ can be represented as invertible square matrices, they are orthogonal projection operators that change the signal representation basis. In this paper, we consider only orthonormal transformations that preserve power relationships.

RIM aims to reduce the impact of $i_{BB}[n]$ on the final cleaned samples, $\tilde{y}[n]$, which are used for the computation of the robust CAF [15]:

$$C_{\psi }(\tau ,f_d)=\sum _{n=0}^{N-1}\tilde{y}[n]c\left(n{T}_s-\tau \right){\mathrm{e}}^{-j2\pi f_{d}n{T}_s}.$$

(10)

In this respect, after RIM, a standard receiver architecture can be adopted where $\tilde{y}[n]$ are used instead of the original samples, $y[n]$.

3. Huber’s Non-Linearity

As mentioned in the previous section, the focus of this paper is Huber’s non-linearity, defined as [26]

$${\psi }_H(Y[k])=\left\{\begin{array}{cc}Y[k]\hfill & \mathrm{for} |Y[k]|\le T_h\hfill \\ T_h\mathrm{csign}(Y[k])\hfill & \mathrm{for} |Y[k]|>T_h\hfill \end{array}.\right.$$

(11)

where $T_h$ is a decision threshold. The ZMNL implements outlier detection on the TD samples, $Y[k]$; if their amplitude is lower than $T_h$ they are considered inliers and they are left unchanged. Otherwise, their amplitude is clipped to $T_h$. The phase of the samples, $Y[k]$, is preserved through the complex signum operator, defined as [27]

$$\mathrm{csign}(Y[k])=\left\{\begin{array}{cc}\frac{Y[k]}{|Y[k]|}\hfill & \mathrm{for} Y[k]\ne 0\hfill \\ 0\hfill & \mathrm{for} Y[k]=0\hfill \end{array}\right..$$

(12)

ZMNL (11) is the derivative of Huber’s penalty function that was originally designed to introduce robustness to the standard Least Squares (LS) estimation [17,25]. The original derivations from Huber considered real input samples [17,25]. We extended the definition of Huber’s ZMNL to complex samples using the phase decoupling approach [20,28].

The Complex signum non-linearity arises when modelling the joint pdf of $i_{BB}[n]$ and ${\eta }_{BB}[n]$ as Laplace [24]. It can be seen that when $|Y[k]|>T_h$, the applied non-linearity is equivalent to the one resulting from the Laplace assumption [16]. The processing implemented by Huber’s non-linearity can thus can be interpreted as a switch between two regimes: the Gaussian and the Laplace ones. This fact is highlighted in Figure 2 which shows the two branches of Huber’s non-linearity.

The amplitude of the input samples, $Y[k]$, determines the switch between Gaussian and Laplace regimes. The decision threshold can be set according to different criteria. A possibility is to consider the passage between the two regimes as a test on the amplitude of the input samples, $Y[k]$. In this respect, (11) can be considered to be a form of interference detection: if the input samples have a magnitude larger than $T_h$, then interference is detected and the complex signum non-linearity is applied. The decision threshold, $T_h$, can be set in order to guarantee a constant false alarm rate [29]:

$$P_{fa}(T_h)=P(|Y[k]|>T_h|{\mathcal{H}}_0)$$

(13)

where ${\mathcal{H}}_0$ denotes the null hypothesis consisting of the absence of interference. Under ${\mathcal{H}}_0$, $Y[k]$ follows a complex circularly symmetric Gaussian distribution, and its amplitude follows a Rayleigh distribution.

In Section 5, we set $T_h$ according to the empirical value determined by Huber [30] for real random variables, whereas we performed a parametric analysis in Section 7 where different threshold settings are considered. In all cases, the threshold, $T_h$, should be obtained by scaling the noise standard deviation, $\sigma$. This fact is better discussed in the next section where the loss of efficiency is analysed. This loss only depends on the normalized threshold, $\frac{T_h}{\sigma }$.

4. Efficiency Analysis

In this section, the loss of efficiency caused by Huber’s non-linearity is theoretically evaluated. The loss of efficiency is the performance degradation caused by the non-linearity in the absence of interference, and it is measured with respect to the optimal Gaussian estimator, i.e., with respect to standard GNSS processing [17]. The loss of efficiency is given by

$$L_0(T_h)=\frac{{\mathrm{SNR}}_{out}^{\psi }}{{\mathrm{SNR}}_{out}}$$

(14)

and is the ratio between the post-correlation Signal-to-Noise Ratio (SNRs) obtained in the absence of interference, i.e., for $i_{BB}[n]=0$. The post-correlation SNR measures the quality of the CAF, and it is defined as [31,32]

$${\mathrm{SNR}}_{out}=\underset{\tau ,f_d}{\max }\frac{{\left|\mathbb{E}\left\{C(\tau ,f_d)\right\}\right|}^2}{\frac{1}{2}\mathrm{Var}\left\{C(\tau ,f_d)\right\}}.$$

(15)

When Huber’s non-linearity is applied, the standard CAF, $C(\tau ,f_d)$, is replaced by its robust counterpart, defined in Section 2.2. Thus, ${\mathrm{SNR}}_{out}^{\psi }$ is obtained according to definition (15), where $C_{\psi }(\tau ,f_d)$ replaces the standard CAF, $C(\tau ,f_d)$.

In order to evaluate ${\mathrm{SNR}}_{out}^{\psi }$, it is necessary to determine the first two moments of the robust CAF, $C_{\psi }(\tau ,f_d)$. In particular,

$$\mathbb{E}\left\{C_{\psi }(\tau ,f_d)\right\}=\sum _{n=0}^{N-1}\mathbb{E}\left\{\tilde{y}[n]\right\}c\left(n{T}_s-\tau \right){\mathrm{e}}^{-j2\pi f_{d}n{T}_s}.$$

(16)

The first moment of the robust CAF is proportional to $\mathbb{E}\left\{\tilde{y}[n]\right\}$, which depends on the transformations, ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$. It is possible to show [23] that when ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$ are orthonormal transformations, the statistical properties of the input samples, $y[n]$, are preserved. Orthonormal transformations preserve the Gaussian nature of the noise affecting the input samples and preserve their variance. Recall that, for the loss of efficiency analysis, we are interested in a nominal, no-jamming situation where the Gaussian assumption is plausible. In addition to this, no correlation is introduced by orthonormal transformations among the input samples. For this reason, the analysis provided in the following is conducted on the input samples, $y[n]$, and ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$ are set equal to the identity. The results obtained are, however, general and apply to orthonormal transformations, such as the DFT and Inverse DFT (IDFT). The case of non-orthonormal transformations is outside the scope of this paper.

When ${\mathbf{T}}_1={\mathbf{T}}_2=\mathbf{I}$, $\tilde{y}[n]={\psi }_H(y[n])$ and

$$\mathbb{E}\left\{{\psi }_H(y[n])\right\}=\mathbb{E}\left\{y[n]\right\}\left[1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}+\frac{\sqrt{\pi }}{2}\frac{T_h}{\sqrt{2}\sigma }\mathrm{erfc}\left(\frac{T_h}{\sqrt{2}\sigma }\right)\right]=\mathbb{E}\left\{y[n]\right\}A\left(\frac{T_h}{\sqrt{2}\sigma }\right)$$

(17)

with

$$A\left(\frac{T_h}{\sqrt{2}\sigma }\right)=1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}+\frac{\sqrt{\pi }}{2}\frac{T_h}{\sqrt{2}\sigma }\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{T_h}{\sqrt{2}\sigma }\right).$$

(18)

Function $\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}(·)$ in (17) and (18) is the complementary error function [33]. The proof of (17) is provided in Appendix A and shows that, on average, the samples processed with Huber’s non-linearity are proportional to the input samples, $y[n]$. On average, the non-linearity causes a scaling that is proportional to (18). The scaling only depends on the normalized threshold, $\frac{T_h}{\sqrt{2}\sigma }$. Combining (16) with (17), it also follows that the robust CAF is, on average, proportional to the standard CAF:

$$\mathbb{E}\left\{C_{\psi }(\tau ,f_d)\right\}=A\left(\frac{T_h}{\sqrt{2}\sigma }\right)\mathbb{E}\left\{C(\tau ,f_d)\right\}.$$

(19)

This result implies that, on average, the CAF is only scaled by Huber’s non-linearity. For this reason, no biases are introduced on the signal parameters that are estimated as the values corresponding to the maximum value of the CAF.

The variance of the robust CAF can be computed using a similar approach and exploiting the independence of the noise samples. In more detail,

$$\begin{array}{cc}\hfill \mathrm{V}\mathrm{a}\mathrm{r}\left\{C_{\psi }(\tau ,f_d)\right\}& =\mathrm{V}\mathrm{a}\mathrm{r}\left\{\sum _{n=0}^{N-1}\tilde{y}[n]c\left(n{T}_s-\tau \right){\mathrm{e}}^{-j2\pi f_{d}n{T}_s}\right\}\hfill \\ & =\sum _{n=0}^{N-1}\mathrm{V}\mathrm{a}\mathrm{r}\left\{\tilde{y}[n]\right\}=\sum _{n=0}^{N-1}\mathrm{V}\mathrm{a}\mathrm{r}\left\{{\psi }_H(y[n])\right\}.\hfill \end{array}$$

(20)

The absence of cross-correlation terms in (20) is due to the independence of the input noise samples and to the fact that ${\psi }_H(·)$ does not introduce memory. The variance of ${\psi }_H(y[n])$ is computed in Appendix A for weak signal conditions and is given by

$$\mathrm{Var}\left\{{\psi }_H(y[n])\right\}=2{\sigma }^2\left[1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}\right]=\mathrm{Var}\left\{y[n]\right\}\left[1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}\right].$$

(21)

Note that the weak signal condition is very plausible in the context of GNSS signals, which are known to be received with very low power levels compared to the noise floor [2]. By combining (20) and (21), it follows that

$$\mathrm{Var}\left\{C_{\psi }(\tau ,f_d)\right\}=\mathrm{Var}\left\{C(\tau ,f_d)\right\}\left[1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}\right]$$

(22)

and

$${\mathrm{SNR}}_{out}^{\psi }=\frac{A^2\left(\frac{T_h}{\sqrt{2}\sigma }\right)}{1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}}\underset{\tau ,f_d}{\max }\frac{{\left|\mathbb{E}\left\{C(\tau ,f_d)\right\}\right|}^2}{\frac{1}{2}\mathrm{Var}\left\{C(\tau ,f_d)\right\}}=L_0\left(T_h\right){\mathrm{SNR}}_{out}.$$

(23)

From the above results, it follows that the loss of efficiency is given by

$$L_0(T_h)=\frac{{\left[1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}+\frac{T_h}{\sqrt{2}\sigma }\frac{\sqrt{\pi }}{2}\mathrm{erfc}\left(\frac{T_h}{\sqrt{2}\sigma }\right)\right]}^2}{1-{\mathrm{e}}^{-\frac{T_h^2}{2{\sigma }^2}}}.$$

(24)

This result was obtained under weak signal conditions, and it is valid for orthonormal transformations. In [23], the validity of (24) was evaluated through simulations considering time domain processing, i.e., with ${\mathbf{T}}_1={\mathbf{T}}_2=\mathbf{I}$. In the following section, the validity of (24) is verified through simulations for frequency domain processing. ${\mathbf{T}}_1$ implements a DFT of size N, while ${\mathbf{T}}_2$ implements the inverse transform, the IDFT.

Validation through Simulations

In order to support the validity of (24) for orthonormal transformations, Monte Carlo simulations were used to compute the post-correlation SNRs and estimate the loss of efficiency. Frequency domain processing was implemented where the DFT and its inverse were used as linear transformations, ${\mathbf{T}}_1$ and ${\mathbf{T}}_2$. The two transforms were normalized in order to preserve norms and signal energy. Notice that in reference [20], the special case of ${\mathbf{T}}_1={\mathbf{T}}_2=\mathbf{I}$ was considered.

The parameters of the signal used for the purpose of showing the loss of efficiency can be consulted in

Table 1. Simulation parameters used for the evaluation of the loss, $L_0(T_h)$, for the case of frequency domain processing.

Parameter	Value
Sampling Frequency	$f_s=4$ MHz
Signal	GPS L1 C/A
Integration Time	1 ms
$C/N_0$	40 dB-Hz
DFT size	4000 samples
No. of Simulation Runs	$4\times {10}^5$

. Particularly, GPS L1 Coarse Acquisition (C/A) signals were generated according to model (2) with $i_{BB}[n]=0$. Simulated signals were used to compute several realizations of the CAF which, in turn, were used to estimate the post-correlation SNR, defined in (15). Simulated signals were processed and integrated using a local code and carrier with the delay and Doppler frequency matching the parameters of the simulated components. When the parameters of the local signal matched those of the incoming samples, the CAF is maximized. The CAF was computed with and without RIM, and both ${\mathrm{SNR}}_{out}$ and ${\mathrm{SNR}}_{out}^{\psi }$ were evaluated. Finally, loss (24) was determined as a function of $T_h$.

The asymptotic analysis of (24) showed that, for large values of $T_h$,

$$L_0(T_h)\to 1,\mathrm{for}{T}_h\to \infty .$$

(25)

This result is consistent with the fact that for large values of $T_h$, ${\psi }_H(·)$ becomes the identity and no processing is applied to the TD samples. If the threshold is too large, the presence of interference is never detected and no processing is applied. In this way, no loss is introduced. For small values of $T_h$, ${\psi }_H(·)$ degenerates to the complex signum non-linearity [24] and

$$L_0(T_h)\to \frac{\pi }{4}=0.7854(-1.049 \mathrm{d}\mathrm{B}),\mathrm{for} T_h\to 0.$$

(26)

A good agreement between simulations and theoretical results is observed in Figure 3. The loss, which is provided in dB scale, monotonically increases as a function of the normalized threshold, $T_h/\sigma$, from a minimum of $-1.049$ dB to a maximum of 0 dB. When the loss is above $3\sigma$, the effect of Huber’s non-linearity is negligible, and the loss approaches unity. The analysis confirms the validity of the loss derived in (24), where the loss of efficiency under nominal conditions is relatively small.

The loss of efficiency in the absence of jamming is expected to be the same for both time and frequency domain processing. This fact is confirmed by the curve obtained by simulations, and is shown in Figure 3. The loss of efficiency determined using DFT/IDFT transformations is the same as that shown in [20] where time domain processing was considered.

5. Simulation Analysis

RIM was validated using a synthetically controlled interference which allowed for adjustment of its parameters. The interference was injected on a real capture of a GPS L1 C/A signal. The signal was captured with a Universal Software Radio Peripheral (USRP) front-end, with a sampling frequency of $f_s=4$ MHz and a front-end bandwidth of 2 MHz. The length of the capture was about 10 s. In the absence of interference and using a standard receiver approach that does not implement RIM, the total number of available satellites was 5, with an average $C/N_0$ of 46 dB-Hz. The Phase Lock Loop (PLL) and Delay Lock Loop (DLL) bandwidths were 10 Hz and 2 Hz, respectively.

Two kinds of jamming signals were injected into the signal capture. On one hand, a CW interference with a fixed frequency at 1 kHz with respect to the central frequency (e.g., the baseband in these experiments) was used, and this is mathematically defined as

$$i_{\mathrm{CW}}[n]=A_J{e}^{j2\pi f_{J}n{T}_s},$$

(27)

where $A_J$ is the amplitude of the jamming signal, and $f_J=1$ kHz is the jamming frequency. On the other hand, the second interference generated was a Saw-Tooth (ST) signal with frequency sweeping from 0 to 4 MHz, repeated periodically every 10 μs. The ST is defined as

$$i_{\mathrm{ST}}[n]=A_J{e}^{j2\pi T_s{\sum }_{j=0}^n{f}_J[j]},$$

(28)

where $f_J[n]$ is the instantaneous frequency.

The performance of the different robust methods was compared with adaptive notch filtering and a standard receiver implementing (4) with no anti-jamming. The assessment of RIM methods included the one proposed in this article, based on Huber’s non-linearity, but also those based on the Laplace and Cauchy assumptions. In the figures, we denote the latter non-linearities by complex signum [16] and myriad [15], respectively. We use the PLL discriminator variance for satellite 11 in the dataset versus Jamming to Noise power ratio ($J/N$) as a metric of tracking quality. Huber’s threshold was set to $T_h=1.345\sigma$, which is known to be optimal for the real-valued signal case [30]. The myriad non-linearity used in this paper was studied in [15] and takes the form

$${\psi }_M(Y[k])=\frac{K}{K+{|Y[k]|}^2}$$

(29)

where $K=6{\sigma }^2$ was selected for the processing performed. This value was selected according to the findings detailed in [15]. Note that all of the robust methods considered here operate in the frequency domain. The $J/N$ was determined by the total noise variance, $2{\sigma }^2$, and the amplitude of the jamming signal:

$$\frac{J}{N}=\frac{A_J^2}{2{\sigma }^2}.$$

(30)

In our approach, we first estimated the noise variance using the clean samples collected in the absence of interference. The $J/N$ was then fixed, and $A_J$ was determined by inverting (30).

Given the superior performance of DLL under high noise conditions, we focused the study on evaluating the impact on the PLL. The results are provided in Figure 4 and Figure 5 for the two jamming signals being tested. The adaptive notch filter used in this paper was the Infinite Impulse Response (IIR) one-pole filter analysed in [34,35]. The transfer function of the filter was as follows:

$$H(Z)=\frac{1-z_0[n]z^{-1}}{1-k_{\alpha }{z}_0[n]z^{-1}},$$

(31)

where $z_0[n]$ is the filter zero, and $k_{\alpha }$ is the pole contraction factor that controls the width of the notch. Note that $z_0$ determines the current notch centre frequency, $\mathsf{\Phi }(n{T}_s)=\frac{f_s}{2\pi }\angle z_0[n]$. Ideally, this frequency should correspond to the interference centre frequency. The adaptive notch filter is indeed designed to adjust $z_0[n]$ in order to track the jamming frequency [34]. In this paper, $k_{\alpha }$ was chosen as $0.7$ for the ST jamming signal and $0.9$ for the CW jamming signal.

Figure 4 shows the results obtained when considering the case of CW interference. The black line represents the threshold for the PLL loss of lock, which occurs when the standard deviation of the discriminator is larger than ${15}^{\circ }$ (corresponding to a variance of $0.068$ rad²) [2]. For low $J/N$ conditions, Huber’s non-linearity and myriad ZMNL performed similarly to standard processing, which is optimal in the absence of interference. As the $J/N$ increases, the variance of the discriminator output, obtained in the absence of mitigation, progressively increased, leading to a loss of lock for a $J/N$ of about 15 dB. At low $J/N$, the adaptive notch filter introduced a small variance degradation. This is due to the fact that the jamming signal was too weak, and the notch filter was unable to estimate its frequency. Note that the adaptive notch filter removesd a portion of the signal spectrum. In this respect, a part of the useful signal component was also removed. When the CW interference was characterized by a $J/N$ greater than 5 dB, the adaptive notch filter was able to track it and effectively remove it. For $J/N$ greater than 20 dB, the notch filter effectively estimated the interference frequency, and the variance of the discriminator output remained almost constant. The notch filter allowed the receiver to operate for larger $J/N$ values than standard processing. In this respect, the receiver did not lose lock under the tested conditions. Robust approaches operate differently from the notch filter and, in particular, they do not need to estimate the frequency of the CW interference. In the frequency domain, the CW has a sparse support, and only a very limited number of frequency samples is affected. The support of the CW signal in the frequency domain does not depend on the $J/N$, and for this reason, the variance of the discriminator output is constant. This effect is clearly visible in Figure 4.

Figure 5 shows the tracking results obtained in the case of a ST jamming signal. For low $J/N$ conditions, the adaptive notch filter performs slightly better with respect to the CW jamming signal case. As the $J/N$ increased, the variance of the discriminator output obtained in the absence of mitigation progressively increased, leading to loss of lock for $J/N$ of about 8 dB. A similar trend was observed for the adaptive notch filter, which was not able to track the variation in the jamming signal. The residual interference components not removed by the notch filter caused a progressive increase in the discriminator variance.

From Figure 5, it is possible to observe that RIM methods provided almost invariant performance in the $J/N$ range analysed. Besides, the discriminator variance was always kept below the loss of lock threshold. RIM allowed receiver operations for all of the conditions tested

6. Experimental Setup

The effectiveness of Huber’s non-linearity for RIM was experimentally evaluated using the data collected during the anechoic chamber experiment described in reference [22,36]. These data had not been used before for the analysis of Huber’s non-linearity and thus they provide results complementary to those obtained in reference [20], which considered different experiments and signals.

As described in [22], the experiments were performed inside a large anechoic chamber available at the Joint Research Centre (JRC) in Ispra, Italy. The tests involved an in-car jammer and a Spirent GSS8800 GNSS hardware simulator. The jammer broadcast a ST signal able to span a frequency range of about 12 MHz in about 9 μs. Its power was varied through a programmable attenuator placed between the jammer antenna connector and the transmitting antenna. The attenuator was used to create a variable jamming power profile—at the beginning of the experiment, the attenuator was set to its maximum value with an attenuation equal to 81 dB. The attenuation was progressively reduced in 2 dB decrements to a minimum value of 45 dB. Finally, the attenuation was again increased in 2 dB increments to its maximum value. A view of the anechoic chamber where the experiment was conducted is provided in Figure 6a, whereas Figure 6b provides a schematic representation of the experimental setup.

The data collected in the anechoic chamber were used to test different interference mitigation techniques [22,36] that could be directly compared with Huber’s non-linearity analysed here. The attenuation profile, the corresponding $J/N$ and additional information on the setup can be found in [22,36], and thus, they are not repeated here.

The Spirent GSS8800 simulator was used to generate GPS L1 C/A and Galileo E1 signals that were recorded along with the jamming signal. A high-fidelity National Instruments (NI) RF signal analyser (NI PXI-5663) was used for the data collection. The parameters adopted for data recording are provided in

Table 2. Parameters of the NI PXI-5663 front-end adopted for the data collection.

Parameter	Value
Centre Frequency	1575.42 MHz
Sampling Frequency	10 MHz
Sampling Type	Complex I&Q
Number of bits	16

.

7. Experimental Results

In this section, the experimental results obtained using the data collected according to the setup described in Section 6 are analysed. Both GPS L1 C/A and Galileo E1c signals were considered. More specifically, the raw In-phase Quadrature (I/Q) data collected with the NI PXI-5663 front-end were processed with a custom Matlab software receiver, implementing RIM and using Huber’s non-linearity. The receiver adopted a standard architecture based on standard acquisition and tracking [2]. The parameters adopted for the processing of the GPS and Galileo signals are provided in

Table 3. Summary of the parameters adopted for the processing of the GPS L1 C/A and Galileo E1c signals.

Parameter	GPS L1 C/A	Galileo E1c
Coherent integration time	1 ms	4 ms
Coherent integration time after bit/secondary code synchronization	20 ms	20 ms
DLL order	2	2
DLL bandwidth	5 Hz	5 Hz
PLL order	3	3
PLL bandwidth	10 Hz	10 Hz

.

RIM was implemented as a pre-correlation technique on the input samples before further processing. Acquisition and tracking were analysed separately. Acquisition performance was assessed at first by comparing the CAFs obtained in the presence of jamming with and without RIM. Tracking was analysed by considering the effective $C/N_0$ [31,32] estimated using the correlator outputs obtained using standard DLLs and PLLs. During tracking, the receiver continuously estimates the signal $C/N_0$, which is determined post-correlation, and it is thus proportional to the post-correlation SNR introduced in Section 4. In particular, the $C/N_0$ is obtained by first estimating (15) and scaling it with respect to the coherent integration time, $T_c$:

$$\frac{C}{N_0}=\frac{1}{2T_c}{\mathrm{SNR}}_{out}.$$

(32)

For this reason, the effective $C/N_0$ reflects the impact of interference, which is perceived for wideband jamming signals as a noise increase, and interference mitigation. In this respect, many research papers [36,37,38,39] express the performance improvements obtained using interference mitigation in terms of $C/N_0$ gain. The $C/N_0$ estimated by the receiver quantifies the ability of a specific technique to mitigate the impact of jamming, and it is used in the following text for performance evaluation.

The $C/N_0$ values obtained considering different mitigation techniques are compared in the following in order to evaluate the performance of Huber’s non-linearity. Both time domain and frequency domain processing are considered for both GPS and Galileo signals. The Galileo E1c signal is a pilot channel. For this reason, a pure PLL with a four-quadrant arctangent is used after achieving secondary code synchronization.

7.1. Time Domain Processing

Time domain processing was implemented using a low-pass filter to make the jamming signal appear as a sequence of pulses [36]. The low-pass filter materialized the transformation, ${\mathbf{T}}_1$, while ${\mathbf{T}}_2$ was set equal to the identity. The characteristics of the low-pass filter were selected according to the properties of the signal to be processed. In this respect, the spectrum of the GPS L1 C/A modulation was narrower than that of the Galileo E1c signal. For this reason, it was possible to use a narrower filter for the processing of the GPS L1 C/A signal.

The CAFs obtained using time domain RIM for the GPS L1 C/A signal are shown in Figure 7, which considers the standard case without mitigation, the usage of the complex signum non-linearity and two configurations adopting Huber’s ZMNL with $T_h=0.5\sigma$ and $T_h=\sigma$, respectively. The CAFs were computed using the dataset collected during the anechoic chamber experiment described in Section 6. In particular, the data used for the CAF computation were extracted from the dataset after 1000 s from the beginning of the experiment. At this point, a $J/N=12$ dB was experienced. A coherent integration time, $T_c=1$ ms, and 10 non-coherent integrations were used for the CAF computation. As shown in Section 5, standard processing and notch filtering were not able to maintain the signal lock for this $J/N$ level. Similarly, for this $J/N$ level, standard processing was unable to reveal the signal presence. This fact clearly emerges from the CAF shown in Figure 7a—no clear signal peak emerges from the CAF noise floor when interference mitigation is not implemented. When RIM is implemented, in this case a low-pass filter with a 1 MHz cut-off frequency is used, a clear signal peak can be identified in the CAFs. The CAFs shown in Figure 7b–d have similar shapes and clearly reveal the presence of the useful signal. Similar results were obtained from the three non-linearities considered in Figure 7. While the CAF analysis provided in Figure 7 is qualitative in nature, it clearly shows the benefits of RIM at the acquisition stage. In the following section, we will focus on tracking performance, evaluated using the effective $C/N_0$ that was estimated post-correlation. This metric allows a more qualitative analysis of RIM and its ability to remove interference terms when computing the correlator outputs. This allows a better comparison between different non-linearities.

The results obtained for the time domain processing of the GPS L1 C/A signal are provided in Figure 8, which analyses the $C/N_0$ values obtained by considering different parameters for Huber’s non-linearity, compared with standard processing and with the complex signum non-linearity. In this case, a low-pass filter with a 1 MHz cut-off frequency was used to protect the main lobe of the GPS L1 C/A modulation. When standard processing was adopted, i.e., without interference mitigation, the receiver was unable to maintain the signal lock for high jamming power levels. In Figure 8, standard tracking lost lock after about 960 s when the $J/N$ was about 12 dB, and front-end saturation started to occur [22]. Figure 8 also provides the $J/N$ estimated from the collected samples. In particular, the noise power, N, was, at first, estimated using the samples collected at the beginning of the experiment when the jamming power, J, could be neglected. The noise power was considered to be constant throughout the whole duration of the experiment, and the remaining samples were then used to estimate the total received power given by the sum of J and N. From the total received power and N, the $J/N$ shown in Figure 8 was finally estimated. After about 900 s from the start of the test, a significant reduction in the steps describing the $J/N$ behaviour was observed. This phenomenon is due to front-end saturation. Indeed the jamming signal was so powerful that it was clipped by the Analog-to-Digital Converter (ADC) quantization function. This led to an apparent reduction in the received jamming power. However, front-end saturation generates harmonics that spread the jamming power on different frequencies reducing the effectiveness of RIM. This phenomenon is equivalent to increasing the number of outliers that affect the received samples. The estimated $J/N$ is reported also in the following figures to provide a better evaluation of the improvements led by RIM. Despite the impact of saturation, RIM provided significant protection against jamming, and the receiver was able to maintain signal lock for all the configurations tested except for the case with $T_h=2\sigma$.

When the decision threshold was too high, as for $T_h=2\sigma$, significant levels of interference passed through undetected and thus were unmitigated. In these cases, the processed samples, $\tilde{y}[n]$, were still affected by interference that compromised the receiver performance. This fact clearly emerges in Figure 8 for the $T_h=2\sigma$ case—while Huber’s non-linearity was able to improve receiver performance for medium-to-high levels of interference, the large detection threshold, $T_h=2\sigma$, did not allow the receiver to maintain signal lock for the whole duration of the test. The receiver performance further increased by reducing $T_h$. This fact clearly emerges from Figure 8—for $T_h=0.5\sigma$ and $T_h=\sigma$, the receiver was able to maintain lock for the whole duration of the test. Small values of $T_h$ increased the robustness to jamming. On the other hand, the efficiency of RIM, i.e., the receiver performance in the absence of interference, decreased with $T_h$. This fact was discussed from a theoretical point of view in Section 4, and it is further analysed in the box in the bottom left part of Figure 8. This box shows the $C/N_0$ values obtained at the beginning of the experiment when jamming levels can be neglected. As expected, standard processing provided the best performance as no loss of efficiency occurred. The $C/N_0$ loss increased with the reduction of $T_h$. For $T_h=2\sigma$, $C/N_0$ values close to those obtained for standard processing case were observed. For small threshold values, $C/N_0$ losses close to 1 dB were observed. These results are in agreement with the theoretical results discussed in Section 4. Figure 8 also provides the $C/N_0$ obtained considering the complex signum non-linearity. As discussed in Section 4, this ZMNL can be considered to be the limit case of Huber’s non-linearity for a vanishing small threshold, i.e., when all the input samples are considered to be affected by interference. In this case, only the phase information is preserved. As expected, the complex signum ZMNL provided the best performance in terms of interference mitigation, but, at the same time, it had the highest loss of efficiency. The threshold of Huber’s non-linearity allows one to select the best compromise between loss of efficiency and robustness to interference. These results are in agreement with the findings presented in reference [20].

The $C/N_0$ values obtained for the processing of the Galileo E1c signal are provided in Figure 9. In this case, the low-pass filter used to isolate the main frequency components of the useful GNSS signal had a cut-off frequency equal to 2 MHz. This design choice preserved the main lobes of the Binary Offset Carrier (BOC) modulation adopted for the Galileo E1c signal.

The results shown in Figure 9 are similar to those discussed above for the GPS L1 C/A modulation. The pure PLL adopted for the processing of the Galileo E1c signal provides improved performance with respect to the Costas loop used for the processing of the GPS L1 C/A modulation. In this respect, the receiver was able to maintain lock for weaker signal conditions. Also, in this case, time domain RIM provided significant performance improvements with respect to standard processing. The selection of the decision threshold, $T_h$, allows one to obtain different compromises between loss of efficiency and robustness to jamming. As for the previous case, resilience to jamming improved for large values of $T_h$, while the efficiency in the absence of interference progressively decreased. The loss of efficiency was, however, limited and in the order of 1 dB.

7.2. Frequency Domain Processing

Frequency domain processing was obtained using the DFT and its fast implementation, the Fast Fourier Transform (FFT), as the first linear transformation, ${\mathbf{T}}_1$, and the IDFT for the materialization of ${\mathbf{T}}_2$. The size of the two transformations was set here to be equal to 1 ms (${10}^4$ samples) for the GPS L1 C/A case and to 4 ms ($4\times {10}^4$) for the processing of the Galileo signals. These values correspond to the length of the code periods of the two GNSS signals.

The results obtained using the frequency domain RIM are provided in Figure 10 for the GPS case. Similar to the time domain case, RIM provided significant improvements with respect to standard processing. In this case, all frequency domain mitigation techniques allowed the receiver to maintain signal lock for the whole duration of the experiment.

The case with $T_h=2\sigma$ had the best efficiency at the expense of a reduced robustness to jamming. On the other side, the complex signum non-linearity provided the highest robustness and the lowest efficiency.

In this case, frequency domain processing was more effective than its time domain counterpart. This fact can be clearly seen in Figure 11 that compares $C/N_0$ time series obtained using time domain and frequency domain processing for $T_h=0.5\sigma$.

The $C/N_0$ curve obtained using frequency domain processing is practically always above the curve obtained using time domain processing. The benefits of frequency domain processing are more evident for medium levels of jamming, before front-end saturation. In such cases, a gain of about 2 dB can be observed. The improved performance of frequency domain processing is due to the fact that the jamming signal affects less samples when projected into the frequency domain than when filtered in the time domain. Increased performance is however paid by the implementation cost of the FFT operations.

The results obtained processing the Galileo E1c signal are provided in Figure 12. The results are consistent with the previous findings where different compromises between robustness and efficiency are obtained by varying the decision threshold, $T_h$.

8. Conclusions

Interference mitigation is crucial to protect GNSS from both intentional or unintentional jamming signals. Due to the widespread use of GNSS technology for personal use, but also in managing critical infrastructures, there is a need to develop techniques that make GNSS resilient to those relatively easy to occur, threats. This paper presented RIM, an interference mitigation framework that exploits the sparsity of jamming signals in a convenient domain (e.g., time or frequency) before correlation with the local codes. The methodology leverages results in robust statistics and expands it when appropriate, for instance to derive the Huber’s non-linearity for complex-valued signals. This paper provided analytical expressions for the loss of efficiency, that is, the degradation of performance caused by the proposed method under nominal conditions where the jamming signal is not present, showing negligible losses. The RIM approach has the desirable feature of avoiding the estimation of the jamming signal waveform, which drastically reduces its usability. Performance was analysed using both synthetic data and an experimental setup, validating the remarkable mitigation results with respect to state-of-the-art techniques.