## Abstract

Evil waveforms (EWFs) are anomalies in signals transmitted by a global navigation satellite system, provoked by electric malfunctions, that can significantly degrade the accuracy of the receiver’s position, velocity, and time solution. In this work, cross-correlation functions of a received signal disturbed by EWF distortion and the locally generated code signal are derived for threat models TM-A, TM-B, and TM-C, with expressions obtained for binary phase shift keying, binary offset carrier (BOC), and composite BOC pilot modulations. Closed-form expressions are offered in terms of sine integral functions for TM-A. A description is presented on how to use these results to assess the performance of EWF detectors through semi-analytic simulation, allowing the detectability and hazard regions to be determined in a computationally efficient way.

## 1 INTRODUCTION

Evil waveforms (EWFs) are anomalies in global navigation satellite system (GNSS)-transmitted signals provoked by electric malfunctions that occur in the signal generators aboard GNSS space vehicles (Phelts, 2001; Selmi et al., 2020; Thevenon et al., 2014). The first reported EWF event occurred in 1993 with Global Positioning System (GPS) pseudorandom noise 19, inducing positioning errors of several meters (Shallberg et al., 2017). These threats manifest themselves in the form of an anomalous correlation peak (Phelts & Akos, 2006). Because the peak is distorted and displaced from the expected position, the processing of EWF-distorted signals by the receiver may lead to a significant loss of accuracy in the position, velocity, and timing solution, thus preventing its use in most applications (Gómez-Casco et al., 2022). Three types of EWFs are usually considered in the literature (ICAO, 2018; Pagot, 2016; Phelts, 2001): threat models TM-A, TM-B, and TM-C, which are associated with digital, analog, and digital plus analog distortion, respectively.

Signal quality monitoring (SQM) techniques have been developed to detect distortions in GNSS signals caused by multipath or spoofing attacks (Ali et al., 2014; Irsigler, 2008). By incorporating additional early-late correlators for monitoring purpose, distortions on the autocorrelation peak can be detected, and corresponding signals can be deweighted (or excluded) (Pirsiavash et al., 2017). Similarly, SQM algorithms may also be employed for detecting the presence of different types of EWF distortion (ICAO, 2018). The performance of an EWF detector may be measured in terms of the probability of not detecting a signal anomaly for a predefined probability of false alarm.

According to Pagot (2016) and Pagot et al. (2018), the main EWF detection algorithms are based on three tests: *T*_{1} (simple ratio), *T*_{2} (difference ratio), and *T*_{3} (sum ratio), which use a bank of early and late correlators. Basically, these tests measure the distortion of the code cross-correlation shape, such as the flatness of the main peak and the correlation asymmetry, for a certain value of *C/N*_{0} and produce a decision concerning the presence/absence of an EWF. The EWF detectability in the threat space and the corresponding hazard regions can then be determined by either using experimental data or resorting to simulations (Negrinho et al., 2021).

The occurrence of EWFs is rare (Shallberg et al., 2017); thus, it is difficult to obtain experimental data that allow one to reliably evaluate SQM algorithms in the presence of EWFs and additive white Gaussian noise (AWGN). Instead, Monte Carlo (MC) simulations may be advantageously applied to assess the performance of the algorithms. MC simulations require very little mathematical analysis and can be applied to almost any communication system, although at the expense of very long simulation run times, especially when the *a priori* probability of a given effect is very small. Consequently, semi-analytic (SA) simulations, in which analytical and simulation techniques are used together, are often proposed as a rapid alternative. These techniques can be applied because the Gaussian channel noise is not affected by the nonlinearities induced by the EWF. Although SA simulations require a higher level of analysis, the payoff is a significantly reduced run time (Tranter et al., 2004).

Following our previous work (Negrinho et al., 2021; Nunes & Sousa, 2022; Sousa et al., 2022), we derive mathematical expressions for signal cross-correlation functions that can be applied to different threat models and modulations in SA simulations. As processing AWGN with the correlators is a linear operation, the result is a set of statistically dependent, zero-mean Gaussian random variables that are added to the various correlator responses to the EWF-distorted waveforms, allowing us to obtain computationally efficient SQM algorithms.

The remainder of this paper is organized as follows. Sections 2 and 3 analytically characterize the multi-correlator outputs for TM-A, TM-B, and TM-C. Section 4 discusses the metrics and tests used by the GNSS receiver for EWF detection. Section 5 presents examples of SA simulations in terms of detectability and hazard regions, whereas Section 6 presents some final remarks.

## 2 THREAT MODEL A

TM-A was originally defined for coarse acquisition (C/A) GPS signals in which the positive chips have a falling edge that leads or lags relative to the correct end-time for that chip (ICAO, 2018). The threat model is associated with a failure in the navigation data unit inside the digital partition of the satellite.

### 2.1 Threat Model TM-A1

This threat model corresponds to the digital distortion 1 described by Pagot et al. (2016): a lead/lag on every falling transition after modulation by the code signal. For this model of distortion, only the lead/lag parameter *d* is required (*d* < 0 indicates a lead and *d* > 0 denotes a lag). Typical waveforms for TM-A1 with sine-phased binary offset carrier (BOCs)(*m*, *m*) signals are presented in Figure 1, where we assume that |*d*| < *T _{c}* (chip duration). This means that the 1 and −1 state durations are not equal, resulting in a distortion of the correlation function and inducing a bias in the pseudorange measurement equal to half the difference in the durations (Raimondi et al., 2012).

Let *ρ*(*t*) ≡ *s*(*t* – *d*) – *s*(*t*). The distorted signal is , where the error signal is as follows:

1

with sign(*x*) equal to −1, 0, or +1 if *x* < 0, *x* = 0, or *x* > 0, respectively. We now define the cross-correlation operator between *ϵ*(*t*) and *s*(*t*) as follows:

2

where *T* denotes the correlation duration. Thus, *R _{ϵs}* (

*τ*) = (1/2)[

*R*(

_{s}*τ*–

*d*) –

*R*(

_{s}*τ*) + 〈|

*ρ*(

*t*)|

*s*(

*t*–

*τ*)〉], where the autocorrelation of

*s*(

*t*) is

*R*(

_{s}*τ*) = 〈

*s*(

*t*)

*s*(

*t*–

*τ*)〉 and 〈|

*ρ*(

*t*)|

*s*(

*t*–

*τ*)〉 ≈ 0 (see Figure 1), leading to

*R*(

_{ϵs}*τ*) = [

*R*(

_{s}*τ*–

*d*) –

*R*(

_{s}*τ*)]/2. The cross-correlation between and

*s*(

*t*) is . The cross-spectrum corresponding to is then as follows:

3

where is the power spectrum of *s*(*t*) and denotes the Fourier transform.

The receiver uses a bank of correlators, with the EWF-distorted signal being correlated with the set of local undistorted replicas. The processing is carried out in baseband using a low-pass filter with transfer function *H*(*f*), as shown in Figure 2. The generic correlator output is given by the following:

4

where the filter output is as follows:

5

Here, *h*(*t*) denotes the filter impulse response. The cross-correlation *R _{A}*(·) accounts for the autocorrelation function of the (baseband) GNSS signal, the distortion due to the EWF, and the filtering effect in the receiver front-end. For simplicity, we consider the receiver’s locally generated signal to be the two-level subcarrier

*g*(

*t*) =

*s*(

*t*) for binary phase shift keying (BPSK) signals and

*g*(

*t*) =

*s*(

_{BOC}*t*) for BOC and composite BOC pilot (CBOCpilot) signals. Thus, the correlation is carried out between the filtered distorted signal and the corresponding locally generated non-distorted and unfiltered version of the signal. However, in the case of transmission of the CBOCpilot signal, we perform the correlation with the two-level BOCs(1,1) carrier, with this simplified scheme resulting in a relatively small performance degradation. Thus, we have the following:

or in the frequency domain:

6

leading to:

7

There are different methods for calculating *R _{A}*(

*τ*) using Equation (7):

- Determine an analytic expression for the integral. The result will be exact, but the calculation may be complex, except for very simple transfer functions

*H*(*f*) and ideal spreading codes with an infinite period.- Calculate the integral numerically. This method is approximate, with the accuracy depending on the length of the integration step.

- Use the discrete Fourier transform (DFT) to compute from the cross-correlation samples , multiply by

*H*(*n*Δ*f*), and apply the inverse DFT (IDFT) to calculate*R*(_{A}*m*Δ*τ*). Here, Δ*τ*and Δ*f*are increments in the time and frequency domains, respectively. This method is also approximate, but its accuracy can be controlled by the size of the DFT sequence (Oppenheim & Schafer, 1975).

Consider now that *g*(*t*) = *s*(*t*) in Equation (7), which is valid for BPSK and BOCs modulations. Using Equation (3) yields the following:

8

Henceforth, for simplicity, we assume that the low-pass filter in Figure 2 is ideal, with bandwidth *B* and zero group delay, such that *H*(*f*) = 1 for |*f*|< *B* and 0 otherwise. Let *x* = *fT _{c}* (normalized frequency) and

*D*=

*d/T*(normalized EWF delay). We then obtain the following:

_{c}9

Note that *R _{A}*(–

*τ*) =

*R*(

_{A}*τ*) (even function) only if the EWF distortion is absent (

*D*= 0).

#### 2.1.1 BPSK(m) Signals

Let us define the triangular pulse Λ_{L}(*x* – *x*_{0}) as follows:

10

with *L* denoting the pulse half-duration and *x*_{0} denoting the pulse center.

If the spreading code is an ideal long spreading code, the autocorrelation function of the BPSK(*m*) signals is (Betz, 2016), where the chip duration is *T _{c}* =

*T*

_{c0}/

*m*,

*m*= 1, 2,… and

*T*

_{c0}= 10

^{−3}/1023 s is the chip duration of the GPS C/A code. The corresponding power spectrum is as follows:

11

with sinc(*x*) = sin(*πx*)/(*πx*). Applying Equation (11) in Equation (9) leads to the following:

12

The integrals in Equation (12) are of the following form:

13

yielding:

14

Appendix A explains how *I*_{s2c}(*a*, *b*; *L*) can be calculated in terms of trigonometric and sine integral functions.

It can be shown that the peak of *R*_{A,BPSK}(*τ*) occurs at *τ* = *DT _{c}*/2 and that the function is symmetric relative to this point. Consequently, the receiver’s delay lock loop (DLL) will exhibit a range error approximately equal to

*cDT*/2, where

_{c}*c*is the speed of light.

The cross-correlation *R*_{A,BPSK}(*τ*) is displayed in Figure 3 for *BT _{c}* = 1 and different values of

*D*. The selected normalized bandwidth

*BT*= 1 corresponds, for instance, to BPSK(10) with a pre-correlation filter of bandwidth 2

_{c}*B*= 20.46 MHz, which is the Galileo receiver reference bandwidth for E5a and E5b bands (European GNSS Agency, 2023).

#### 2.1.2 BOCs(m, m) Signals

The power spectrum of the BOCs(*m*, *m*) modulation is given by the following (Sousa & Nunes, 2013):

15

Substituting Equation (15) in Equation (9) leads to the following:

16

Figure 4 shows cross-correlation curves obtained from Equation (16) for different lags *D* and *BT _{c}* = 10.

#### 2.1.3 CBOCpilot Signals

Here, the distorted signal is , with the CBOC(6,1,1/11) pilot signal being defined as *s*(*t*) = [*αp*_{11}(*t*) – *βp*_{61}(*t*)]*c*(*t*) = *αs*_{11}(*t*) – *βs*_{61}(*t*), where , , *p*_{11}(*t*) and *p*_{61}(*t*) are, respectively, the BOCs(1,1) and BOCs(6,1) subcarriers, and *c*(*t*) is the code signal. Here, *s*_{11}(*t*) = *c*(*t*)*p*_{11}(*t*) and *s*_{61}(*t*) = *c*(*t*)*p*_{61}(*t*) are, respectively, the BOCs(1,1) and BOCs(6,1) signals (European GNSS Agency, 2023). The receiver correlation scheme is shown in Figure 2, with the local generator given by *g*(*t*) = *s*_{11}(*t*). As mentioned above, this correlation scheme is different from those of the BPSK and BOCs modulations, as the signal of the local oscillator is not a replica of the nominal transmitted signal, with the purpose of simplifying the receiver architecture. For this case, we have the following:

17

with (Sousa & Nunes, 2013), where Λ_{L}(|*τ*| – *τ*_{0}) = Λ_{L}(*τ* + *τ*_{0}) + Λ_{L}(*τ* – *τ*_{0}) denotes two triangular pulses of duration 2*L* centered at *τ* = ±*τ*_{0}.

The cross-correlation between signals *s*_{61}(*t*) and *s*_{11}(*t*) is given by the following (Sousa & Nunes, 2013):

18

To determine *R*_{ϵ,s11}(*τ*), we must account for Equation (1), which is equivalent to *ϵ*(*t*) = *s*(*t* – *d*) – *s*(*t*) if *s*(*t* – *d*) ≥ *s*(*t*) and *ϵ*(*t*) = 0 otherwise. Thus, we obtain the following:

19

Now, using Equation (17) yields . The cross-spectral density corresponding to is as follows:

20

The cross-correlation at the receiver is given by Equation (7):

21

with *T _{c}* = 1/1.023

*μs*or:

22

where *R*_{A,BOCs(1,1)}(*τ*) is determined by Equation (16) and:

23

Now, using Equation (18), the cross-spectrum of *s*_{11}(*t*) and *s*_{61}(*t*) is as follows:

24

Substituting Equation (24) in Equation (23) yields the following:

25

### 2.2 Threat Model TM-A2

In this threat model, the distortion in the subcarriers occurs before modulation by the code signal. This distortion corresponds to the digital distortion 2 described by Pagot et al. (2016), being applicable to signalings in which one or more subcarriers are separated from the code signal, namely, BOCs(*m*, *m*) and CBOCpilot signals.

#### 2.2.1 BOCs(m, m) Signals

Let parameter *d*_{11} express the subcarrier lead/lag. The distorted subcarrier is , and the error signal is as follows:

26

Note that for |*d*_{11}| > *T _{c}*/2, the error model in Equation (26) may not be realistic. The modulated signal is

*s*(

*t*) =

*p*

_{11}(

*t*)

*c*(

*t*), and its distorted version is , where

*c*(

*t*) is the code signal. The receiver cross-correlation with the local BOCs signal is , which can be written as follows:

27

where *R _{s}*(

*τ*) ≡ 〈

*s*(

*t*)

*s*(

*t*–

*τ*)〉. The cross-correlation

*R*

_{ϵ,11}(

*τ*) ≡ 〈

*ϵ*

_{11}(

*t*)

*c*(

*t*)

*s*(

*t*–

*τ*)〉 is depicted in Figure 5 and can be expressed by the following difference in triangles:

28

In the frequency domain, Equation (27) can be written as follows:

29

with *G _{s}*(

*f*) given by Equation (15). The cross-spectrum

*G*

_{ϵ,11}(

*f*) is determined from Equation (28), leading to the following:

30

Using Equation (29), we obtain the following:

31

with *D*_{11} ≡ *d*_{11}/*T _{c}* and

*T*= 1/1.023

_{c}*μs*. Thus, we have the following:

32

In the absence of EWF distortion, i.e., for *D*_{11} = 0, *R _{A}*(

*τ*) is an even function.

#### 2.2.2 CBOCpilot Signals

Two parameters are required to express the lead/lag of each subcarrier: *d*_{11} and *d*_{61}. These delays are assumed to be independent.

The distorted transmitted signal is , with the distorted subcarriers being for BOCs(1,1) and for BOCs(6,1). The error signals are *ϵ*_{11}(*t*), given by Equation (26), and:

33

Note that for |*d*_{61}| > *T _{c}*/12, the error model in Equation (33) may not be realistic. The receiver cross-correlation with the local BOCs signal is , which can be written as follows:

34

In the frequency domain, Equation (34) can be written as follows:

35

with *G*_{11}(*f*), *G*_{61,11}(*f*), and *G*_{ϵ,11}(*f*) given, respectively, by Equations (15), (24), and (30). Moreover, the following can be shown:

36

leading to:

37

Accounting for Equation (35) yields the following:

38

with *R*_{A,BOCs}(*τ*) and *βC*(*τ*) given, respectively, by Equations (16) and (25). This finally leads to the following result:

39

In the absence of EWF distortion, i.e., for *D*_{61} = *D*_{11} = 0, *R _{A}*(

*τ*) is an even function.

Figure 6 displays the cross-correlation *R _{A}*(

*τ*) for

*D*

_{11}= 0.15 and

*D*

_{61}= 0.05. The cross-correlation for

*D*

_{11}=

*D*

_{61}= 0 is also shown (non-distorted signal).

The curves in Figures 3, 4, and 6 show that a code discriminator with early-late spacing Δ_{EL} > 0 will exhibit synchronization errors when the EWF is present because the position of the curves is shifted, in addition to being distorted. This effect will provoke positioning errors in the receiver. Later, we will describe methods that use the curve distortion to detect the presence of EWFs.

## 3 THREAT MODELS B AND C

TM-B introduces amplitude modulation and model degradations in the analog section of a satellite. More specifically, this model consists of the output from a second-order system when the code-modulated baseband signal is the input (ICAO, 2018; Phelts, 2001). TM-B assumes that the degraded satellite subsystem can be described as a linear system dominated by a pair of complex conjugate poles. The poles are located at *σ* ± *jω _{d}*, where

*σ*is the damping factor in nepers/second and

*f*=

_{d}*ω*/(2

_{d}*π*) is the resonant frequency in Hertz. For

*t*≥ 0, the system response to the step function, with a discontinuity at

*t*= 0, is

*y*(

*t*) = 1 – exp(–

*σt*)[cos(

*ω*) + (

_{d}t*σ/ω*)sin(

_{d}*ω*] (D’Azzo et al., 2003). Let

_{d}t*ω*and define, respectively, the (undamped) natural frequency and the damping ratio (with

_{n}*ξ*< 1), such that

*σ*=

*ξω*and . The transfer function of the second-order low-pass filter describing the damped linear system is given as follows:

_{n}40

with Ω ≡ 2*πf/ω _{n}*. Next, we present expressions for the cross-correlation

*R*(

_{A}*τ*) obtained in the scenario corresponding to TM-B:

41

yielding:

42

where *ρ* ≡ 2*πB/ω _{n}* and

*μ*≡

*ω*/(2

_{n}T_{c}*π*).

### 3.1 BPSK(*m*) Signals

The power spectrum of BPSK signals is given by Equation (11). Substituting this result in Equation (42) leads to the following:

43

with:

44

In Appendix B, it is shown that the cross-correlation *R _{A}*(

*τ*) can be written as an analytic expression that involves sine and cosine integral functions. Figure 7 illustrates the correlation deformation induced by a TM-B-type distortion. The nominal signal in the figure refers to the cross-correlation in the absence of EWF distortion, which is given in Equation (14) with

*D*= 0.

In addition to deformation, the EWF also introduces a delay in the cross-correlation, which leads to range errors in the receiver’s DLL. Figure 8 exhibits the range errors in meters for different values of the damped frequency *f _{d}* =

*ω*/(2

_{d}*π*) and damping factor

*σ*. Note that the largest errors occur when

*f*and

_{d}*σ*are simultaneously small. For simplicity, we assume that the receiver’s DLL tracks the peak of

*R*(

_{A}*τ*).

### 3.2 BOCs(*m*, *m*) Signals

Using Equation (15), we can derive the following (see also the work by Sousa and Nunes (2013)):

45

Substituting this result in Equation (42) leads to the following:

46

where:

47

Figure 9 illustrates the correlation deformation induced by a TM-B-type distortion. The nominal signal in the figure refers to the cross-correlation in the absence of a TM-B-type distortion, which results from doing *D* = 0 in Equation (16).

### 3.3 CBOCpilot Signals

For a CBOCpilot signal affected by TM-B distortion, we can apply Equation (21) with *d* = 0, which yields the following:

48

We now incorporate Equations (45) and (24) in Equation (48) and utilize the following definition:

49

We use the following trigonometric equality:

50

which allows us to establish the relation between *K _{2}* and

*K*

_{0}as follows:

51

Therefore, we obtain the following:

52

The cross-correlations *R _{A}*(

*τ*) of TM-B, defined with the help of functions

*K*

_{0}(

*μ*,

*ξ*;

*x*),

*K*

_{1}(

*μ*,

*ξ*;

*x*), and

*K*

_{2}(

*μ*,

*ξ*;

*x*), can be expressed through analytic formulas using sine and cosine integral functions of complex arguments, as shown in Appendix B for BPSK signals. However, in contrast to the TM-A case, the derivation of cross-correlations for TM-B is generally complicated. An alternative straigh-forward solution resorts to the direct numerical integration of Equations (43), (46), and (52). The IDFT/DFT method, referred to in Section 2, may be used as well.

TM-C introduces both lead/lag and amplitude modulation. Specifically, this model consists of outputs from a second-order system when the modulated baseband signal suffers from lead or lag. This waveform includes a combination of TM-A and TM-B effects (ICAO, 2018). Thus, the main difference between the expressions obtained for TM-A and those to be derived for TM-C is that Equation (6) is now replaced by the following:

53

where *H*(*f*) represents the channel filter (essentially due to the receiver’s front-end filtering effect) and , given by Equation (40), is due to the degradation introduced in the analog section of the satellite.

## 4 METRICS AND TESTS

Let us consider the scheme shown in Figure 10 for the detection of EWF-distorted GNSS signals constituted by a bank of 2*N* + 1 correlators, with *N* > 1 and a correlator spacing equal to Δ. The phase *ϵ* denotes the code synchronization error. In the work by Pagot et al. (2016), the authors noted that the overall delay interval 2*N*Δ of the correlator bank need not be equal to twice the chip duration (*T _{c}*) for detection purposes. In fact,

*N*Δ ≈ 0.25

*T*is generally sufficient because (1) TM-like distortions are more visible around the prompt of the correlation function and (2) correlator outputs situated at too great a distance from the prompt are more subject to multipath.

_{c}The correlator spacing is typically lower-bounded by Δ ≈ 10 ns because the correlator outputs with very small Δ values are strongly correlated (and thus redundant), due to the effect of the front-end filter. Moreover, a time delay of 10 ns between correlators is now reachable, but lower values may be difficult to achieve by analog-to-digital converters.

In Figure 10, *r*(*t*) is the (band-pass) received signal, which includes the GNSS signal of power *P* and AWGN of power spectral density *N*_{0}/2. The corresponding carrier-to-noise ratio is (*C/N*_{0}) = *P/N*_{0}. The baseband signal is , where *h*(*t*) is the impulse response of the receiver filter and *n*(*t*) is Gaussian low-pass noise with power spectral density *G _{n}*(

*f*) =

*N*

_{0}|

*H*(

*f*)|

^{2}=

*N*

_{0}if |

*f*| ≤

*B*and

*G*(

_{n}*f*) = 0 otherwise. We assume that 2

*B*is the front-end bandwidth and that the receiver’s phase-tracking loop is closed with a sufficient rate to render the effect of the Doppler frequency negligible along the processing chain.

Because 2*BT* ≫ 1, *n*(*t*) is generally well approximated by white noise, that is, *E*{*n*(*t*)*n*(*λ*)} ≈ *N*_{0}*δ*(*t* – *λ*), with *E*{·} denoting the expectation. The correlator outputs are Gaussian random variables defined as follows:

54

where the additive noises:

55

are zero-mean, Gaussian random variables, with the following correlations:

56

For simulation purposes, the noise vector *n* = [*n*_{–N}…*n*_{0},…*n _{N}*]

^{T}may be easily generated by using the procedure presented in Appendix C.

The SQM consists of a test to determine whether the signal is affected by EWF distortion. Let us consider the hypotheses *H*_{0} (signal is not distorted) and *H*_{1} (signal is distorted). Let denote the test variable, which is normalized such that under hypothesis *H*_{0} and typically under hypothesis *H*_{1}. The result of the test is binary: if , it is assumed that the signal is not distorted; otherwise, we assume that the signal is distorted. For a given test, we define the following probabilities:

57a

57b

It is not possible to simultaneously minimize *P _{fa}* and

*P*. A workable solution to this issue is the Neyman–Pearson criterion, which consists of fixing

_{md}*P*at a preselected value and then computing

_{fa}*P*(Kay, 1998).

_{md}For a given metric *μ* (which is a function of the correlator outputs), the corresponding test variable is defined as follows:

58

where *μ*_{mea} is the measured value of the metric, which is disturbed by thermal noise and may or may not be affected by signal distortion (anomaly), *μ*_{nom} is the nominal value of the metric (without additive noise or distortion), and *λ* is a normalization parameter that depends on *P _{fa}*. First,

*λ*is adjusted to achieve a desired value of

*P*. Then,

_{fa}*P*is determined from the values of the correlator measurements. This probability depends on the amount of distortion in the autocorrelation function, the carrier-to-noise ratio, etc.

_{md}In work by Pagot et al. (2016), three tests based on ratio metrics were proposed: simple ratio, difference ratio, and sum ratio. These metrics will be analyzed next. The ratio tests (or tests based on ratio metrics) specifically attempt to detect the presence of dead zones (flat correlation peaks) and abnormally sharp or elevated correlation peaks (Phelts, 2001).

### 4.1 Simple Ratio Metric

Assuming that there are 2*N* + 1 correlators, it is possible to form 2*N* simple ratio metrics , *i* = ±1, …, ± *N*. The tests are based on random variables expressed by , where . Therefore, .

Under hypothesis *H*_{0} (no EWF), the mean and variance of can be approximately calculated using the following second-order Taylor series expansions (Papoulis & Pillai, 2002):

59

60

where and are Gaussian random variables with correlation coefficient *r* = (*E*{*XY*} – *m _{X}m_{Y}*)/(

*σ*). Considering Equations (54) and (56), we have , , , and

_{X}σ_{Y}*r*=

*R*(

_{s}*i*Δ)/

*R*(0). Thus, we obtain the following:

_{s}61

When (*C/N*_{0})*T* ≫ 1, we can simplify Equation (61), yielding , which is a convenient result, as it does not depend on the carrier-to-noise ratio. We also have the following:

62

with . Assuming that is Gaussian under hypothesis *H*_{0}, we obtain the following probability of false alarm for test :

63

where erfc(·) denotes the complementary error function:

64

For the sake of simplicity, let us consider that the 2*N* parameters are individually adjusted to achieve a common probability of false alarm (*P*_{fa}), that is, , *i* = ±1, …, ± *N*. The decision criterion consists of declaring an anomaly if at least one of the tests is positive . The corresponding probability of a false alarm is . Using the union bound approximation (Proakis & Salehi, 2008) yields the following:

65

Taking into account Equations (63) and (65), we obtain the normalization parameter for each pair of correlators:

66

### 4.2 Difference Ratio Metric

Considering 2*N* +1 correlators, i.e., *N* pairs of early-late correlators, we can form *N* difference ratio metrics , *i* = 1,…, *N*. However, . Thus, the test variables are , with and:

67

Assuming that is Guassian, the probability of a false alarm for test is given by an expression similar to Equation (63). As in the simple ratio metric, we assume that all tests are characterized by a common probability of false alarm (*P _{fa}*), that is, . For each pair of early-late correlators, we obtain a different value of the normalization parameter:

68

### 4.3 Sum Ratio Metric

Considering 2*N* + 1 correlators, we can form *N* sum ratio metrics , *i* = 1,…, *N*. The test variables are , where . In the absence of EWF distortion, the tests have zero mean and variance:

69

Assuming that is Gaussian, the probabilities of false alarm are expressed by an expression similar to Equation (63), and we obtain a different value of the normalization parameter for each pair of early-late correlators:

70

## 5 SIMULATION RESULTS

The ratio tests discussed in the previous section allow for varying performance in EWF detection, depending on the characteristics of the distorted cross-correlation *R _{A}*(

*τ*). In general, tests based on simple ratio or sum ratio metrics are used to detect deformations in the cross-correlation that result in either increased flatness or sharpness. Tests based on difference ratio metrics excel at detecting asymmetries in the cross-correlation curve but tend to perform poorly when the cross-correlation, although flattened, remains symmetric. The receiver bandwidth also plays an important role in EWF detection. Figure 11 displays the regions of detectability for TM-A1 and BPSK(1) modulation for

*P*= 10

_{FA}^{−4},

*P*≤ 10

_{MD}^{−3}, and a bank of 2

*N*+ 1 = 5 correlators with spacing Δ = 0.15

*T*. Two different bandwidths were considered:

_{c}*BT*= 2 and 5. The darkest region of each plot indicates the region of detectability. The plots show that the interval of detectability (

_{c}*D*≥

*D*

_{min}), defined in terms of the minimum value of parameter

*D*, increases with the value of (

*C/N*

_{0})

*T*in both plots. Additionally, for a given value of (

*C/N*

_{0})

*T*, the interval of detectability grows with the bandwidth.

As noted in Section 2, the presence of a TM-A1 EWF anomaly provokes a deformation and a time shift in the cross-correlation *R _{A}*(

*τ*) equal to

*DT*/2, which is responsible for the range error

_{c}*cDT*/2. If we consider that the error for a certain application must not exceed the quantity maximum-allowable error in range (MERR), EWF anomalities with

_{c}*D*values less than 2 · MERR/(

*cT*) may not be detected, but the corresponding range errors will not be hazardous (Rife & Phelts, 2006). Figure 12 displays the hazard regions (darker regions), corresponding to Figure 11, for which the range error exceeds the MERR value of 2 m and the EWF anomaly is not detected by the test.

_{c}Figure 13 presents the detectability regions (darker regions) for TM-B with BPSK(1) modulations in the threat space region of 0 < *f _{d}* < 20 MHz and 0 <

*σ*< 200 Mnepers/s, assuming two signal-to-noise ratios: (

*C/N*

_{0})

*T*= 40 dB and 48 dB. As expected, the detectability region grows with increasing signal-to-noise ratio, starting with small values of

*f*and

_{d}*σ*. To obtain the corresponding hazard regions, for a given MERR value, the plots in Figure 13 could be combined with the range error plot in Figure 8.

## 6 FINAL REMARKS

In the first part of this work, mathematical expressions were derived for signal cross-correlation functions that are applicable to different threat models and modulations. For TM-A, simple closed-form expressions were derived that can be rapidly computed via sine integral functions with real arguments. For TM-B, only results expressed in terms of integrals could be obtained with generality. When available, specific expressions for the cross-correlation allow for a significant reduction in the computational effort required by SA simulations of different EWF scenarios, with an emphasis on determining the detectability and hazard regions of the threat space. Nevertheless, even when closed-form solutions are not available, numerical integration generally constitutes a rapid computational alternative.

In the second part of this work, a bank of *N* pairs of early-late (symmetric) correlators combined with a prompt correlator was used with a test based on simple, difference, or sum ratio metrics. Alternatively, we may devise a global test in which the three metrics are used simultaneously, as proposed by Sousa et al. (2022). Although this strategy may introduce some degree of redundancy, it prevents the need to select the optimal test for each type of EWF anomaly, as the nature of the anomalies is usually not known beforehand.

This work permits us to draw the following conclusions:

The expressions obtained herein for calculating the receiver’s cross-correlation are general, providing a useful means for evaluating the performance of current and future EWF detection algorithms via SA simulations. Certain closed-form results assume ideal band-pass filtering, which is considered a good approximation to high-order real filters. The SA results presented herein provide useful insights for designing effective monitoring schemes and examining the sensitivity of receivers to specific threats and errors, without incurring a large computational effort. Note that if the correlations were calculated from signal samples (simple MC simulation), the processing burden would be higher, depending on the sample rate and number of bits per sample.

The EWF detection performance improves with the number of correlator pairs

*N*; however, for*N*> 4, we did not obtain any significant improvement. However, one must be cautious, as this condition may not hold in all scenarios. Moreover, using equal spacing between correlators may not be the optimal solution.The regions of detectability tend to increase with the growth of the coherent integration interval

*T*. In performance analyses, the product (*C/N*_{0})*T*is an important parameter to consider. Consequently, for small values of (*C/N*_{0}), large integration intervals are required. However, even for reasonably large values of*T*, anomaly detectability may not be guaranteed in significant parts of the threat space. In addition, large values of*T*imply large delays (latency) in the anomaly detection.The usefulness of these EWF detection algorithms may be jeopardized in the presence of multipath, as discussed by Nunes and Sousa (2022). Therefore, their applicability to mobile receivers in harsh environments is questionable. Note that fixed reference stations generally permit one to obtain better EWF detection performance than mobile receivers because fixed reference stations are not affected by random dynamic stress and the multipath effect can be made negligible.

Given the limitations of EWF detectors using the SQM approach, the aid of other techniques may be necessary, such as receiver autonomous monitoring techniques, which operate at the navigation solution level.

## HOW TO CITE THIS ARTICLE

Nunes, F. D., & Sousa, F. M. G. (2024). Efficient signal quality monitoring of GNSS signals disturbed by evil waveforms. *NAVIGATION, 71*(4). https://doi.org/10.33012/navi.668

## ACKNOWLEDGMENTS

This work was partially supported by FCT-Fundação para a Ciência e Tecnologia, I.P., with project reference UIDB/50008/2020 and DOI identifier https://doi.org/10.54499/UIDB/50008/2020. This research was conducted in part under the European Space Agency project MONINT (Monitoring Solution for New Integrity Applications), with GMVIS, S.A. as the prime contractor and Instituto de Telecomunicações as a subcontractor. The authors express appreciation for the collaboration of A. Negrinho, P. Boto, and P. Fernandes of GMVIS.

## A COMPUTATION OF INTEGRAL *I*_{S2C}

_{S2C}

Let us consider the following integral:

71

For BPSK(m) modulations, we have *a* = 1, and for BOCs(m,m) modulations, we have *a* = 1/2. Because sin^{2} *x* = (1 – cos(2*x*)) /2, we can write the following:

72

We can use the following relation:

73

to derive:

74

75

This result permits us to compute *I _{s2c}*(

*a*,

*b*;

*L*) in terms of trigonometric and sine integral functions, where the sine integral is defined as follows (Abramowitz & Stegun, 1970):

76

Because the arguments of the sine integrals in Equation (75) are real, these functions can be computed with very small errors using the polynomial approximations given, for instance, by Abramowitz and Stegun (1970) and Rowe et al. (2015).

## B ANALYTIC FORMULA OF *R*_{A,BPSK}(*τ*) FOR TM-B

The cross-correlation *R*_{A,BPSK}(*τ*) in Equation (43) can be expressed analytically as follows:

77

where *a* = *πμ*, *b* = 2*πμ*(*τ*/*T _{c}*),

*c*= 4

*ξ*

^{2}– 2, and:

78

79

80

The quantities *Y _{m}* are calculated in terms of sine and cosine integral functions of complex argument Si(

*z*) and Ci(

*z*), respectively (see, for instance, Abramowitz and Stegun (1970)):

81

and:

82

with:

83

We also have the following:

84

and:

85

with:

86

The complex numbers *ω _{k}*,

*k*= 1,…, 4 are the roots of

*ω*

^{4}+

*cω*

^{2}+ 1 = 0.

## C GENERATION OF ARBITRARY GAUSSIAN VECTORS

Consider the generation of a zero-mean Gaussian noise vector *U* = *GW* = [*u*_{1},…, *u _{M}*] with arbitrary covariance matrix

*C*=

*E*{

*UU*}, using the zero-mean noise vector

^{T}*W*= [

*w*

_{1},…,

*w*] with independent Gaussian components of unity variance. This problem occurs, for instance, in MC simulations. The solution consists of determining

_{M}*G*from

*C*=

*GE*{

*WW*}

^{T}*G*=

^{T}*GG*. Matrix

^{T}*C*can be written in terms of its eigenvalues and eigenvectors according to

*C*=

*V*Λ

*V*, where the columns of

^{T}*V*are orthonormal eigenvectors and Λ = diag{

*Λ*

_{1},

*λ*

_{2},…,

*λ*} is a diagonal matrix with the eigenvalues of

_{M}*C*. Finally, .

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.