## Abstract

Low-cost global navigation satellite system and global positioning system(GPS) receivers require reliable baseband processing to guarantee accurate positioning. However, classic baseband performance is limited in challenging cases due to the characteristics of traditional loop filters. Accordingly, a snapshot baseband maximum likelihood estimator (MLE) using super-long coherent integration (S-LCI) in a fractional Fourier domain (FrFD) is proposed to upgrade the traditional frequency/phase/delay lock loop tracking algorithms. First, applying the S-LCI correlation in an FrFD increases the accuracy of a weak and dynamic signal estimation. Tolerance of the initial guess error in the snapshot baseband processing is then relaxed by the MLE. Finally, a gradient descent algorithm accelerates the convergence of signal estimation. Moreover, we derive the Cramer-Rao lower bound for the proposed MLE. Both numerical simulations and real-world experiments based on this GPS receiver prototype verify the effectiveness of its high-accuracy estimations of weak signals, strong tolerance for large initial guess errors, and prompt responses to converging.

- baseband processing
- Cramer-Rao lower bound
- fractional Fourier transform
- maximum likelihood estimator
- snapshot GNSS receiver
- super-long coherent integration

## 1 INTRODUCTION

Accurate positioning with low-cost global navigation satellite system (GNSS)/global positioning system (GPS) devices currently plays an important role in our daily lives. Examples include navigation for pedestrians via smartphones (Lachapelle & Gratton, 2019; Luo et al., 2019c; Ng et al., 2021) and wearable sensors (Faragher et al., 2018, 2019) and autonomous ground vehicles (Humphreys et al., 2020). These applications depend on low-cost antennas and hardware. There will be significant demand for consumer-level GNSS receivers in the next ten years (*EUSPA EO and GNSS Market Report*, 2022). GNSS signals are confronted with frequent interruptions and reflections in challenging environments. Thus, incoming GNSS signals are often weak and dynamic, which causes the traditional GNSS baseband architecture to be very fragile (Luo et al., 2018b; Ren & Petovello, 2017; Yan et al., 2017).

The classic GNSS baseband architecture consists of code/carrier numerically controlled oscillators (NCOs), correlators, code/carrier phase discriminators, and code/carrier loop filters. A phase/frequency lock loop (PLL/FLL) and a delay lock loop (DLL) separately process the code and carrier components of GNSS signals, respectively (Dierendonck, 1996). To enhance the sensitivity of weak signal processing, many previous studies have applied vector tracking and ultra-tight integration techniques to the traditional baseband architecture (Lashley et al., 2009; Pany & Eissfeller, 2006; Petovello et al., 2008). These techniques have recently been upgraded in the GNSS baseband as a reaction to challenging modern environments (Maier & Pany, 2021; Watts et al., 2021). However, these methods are based on the classic baseband architecture which has not changed in the past 40 years. At present, some practical issues need to be addressed, for example:

First, there is a need to work with weak and dynamic signals in a low-cost front-end application and to determine an appropriate integration time and a loop filter bandwidth, factors that are especially challenging for achieving high accuracy;

Second, this will require a relatively accurate initial code phase and Doppler frequency to guarantee a reliable lock at the beginning of tracking; and

Third, at the start of the tracking process, it takes some time for loop filters to converge.

Some existing research studies have focused on these problems. While ignoring signal dynamics, the long coherent integration (LCI) technique, in which the coherent integration interval is usually from 20 to 200 ms, is an ideal method to increase the signal-to-noise ratio (SNR) (Ren & Petovello, 2017). The more extended integration naturally has a more vital capacity to separate the line-of-sight (LOS) and non-line-of-sight (NLOS) signals in a multipath channel (Gowdayyanadoddi et al., 2015; Luo et al., 2021; Wang & Morton, 2020; Xie & Petovello, 2015). Signal processing benefits from the LCI, particularly in challenging environments. The more real signal power is exploited, the more accurate the local replicas of the LOS signal will be (Groves et al., 2020; Petovello et al., 2008).

As an upgraded version, a super-long coherent integration (S-LCI) algorithm in which the integration interval is greater than 200 ms, has recently attracted increasing attention among those focused on low-cost receiver design (Faragher et al., 2018, 2019). This can be regarded as evidence that LCI/S-LCI is urgently needed. However, the navigation data modulated in older GNSS signals can generate a bit-sign transition that hinders power accumulation in the LCI/S-LCI implementation for the traditional GNSS signals. The bit sign uncertainty can be eliminated using the recorded data bits stream. For example, the GFZ German Research Centre for Geosciences built a service that offers navigation data information for GPS L1 C/A (Beyerle et al., 2009). This state-of-the-art technique uses a combined S-LCI method and has been proposed to mitigate strong interference in the GNSS signal acquisition (Yang et al., 2022). While this work leverages a phase rotation method to compensate for the carrier phase migration generated by the frequency error in the S-LCI, it cannot correct the phase error due to the Doppler rate.

The GNSS modernization has shown its superiority with respect to S-LCI implementation. Unlike the old GPS L1 C/A signal with only a single data channel, the modernized GPS L5 signal is composed of two channels, including an in-phase data channel and a quadrature data-less pilot channel (Macchi-Gernot et al., 2010; Tran, 2004). A known secondary code sequence is modulated in the pilot channel. In other words, the LCI/S-LCI can be realized with modernized signals without bit-sign estimation. Moreover, mainstream modernized GNSS signals are becoming available worldwide; 17 out of the 31 GPS satellites transmitting L5 signals were operational as of December 2021 (*GPS Constellation Status*, 2022). BDS B2a and Galileo E5a signals have modulations and structures that are very similar to those of the GPS L5 signal, with 26 and 24 operational satellites, respectively (*BDS Constellation Status*, 2021; *Galileo Constellation Information*, 2021). More satellites transmitting these signals will be launched in the future. Collectively, these findings indicate that modernized signals have the potential to be widely applied for users’ navigation and positioning in commercial markets (Givhan et al., 2020; Ng et al., 2021).

By applying the LCI/S-LCI to the modernized signal, the GNSS receivers can be improved. However, the baseband becomes more vulnerable to dynamic conditions. As the Doppler rate of the incoming signal increases, the correlation peak after the LCI/S-LCI decreases more easily. In the field of signal modeling with a higher order of dynamics, a fractional Fourier transform (FrFT) is commonly used, such as in radar signal processing (Filip-Dhaubhadel & Shutin, 2020). Recent research shows that FrFT is superior to traditional Fourier transform (FT) for suppressing noise power (Ozturk et al., 2021). A digital FrFT as a practical implementation was also proposed earlier (Ozaktas et al., 1996). As FrFT surpasses FT in the baseband signal processing, the authors have previously performed experiments designed to verify its effectiveness in GNSS. For example, these experiments first verified that the digital FrFT exhibited a shorter mean acquisition time than the fast Fourier transform (FFT) and could efficiently improve the tracking performance in a dynamic situation (Luo et al., 2018a; Luo et al., 2018b). Next, the accuracy of Doppler rate estimation could also be optimized without raising the computational load by the closed-form model of the correlation peak based on the digital FrFT (Luo et al., 2019b). Results from the most recent research study revealed that S-LCI based on the FrFT improved the time of arrival (TOA) estimation in a multipath channel (Luo et al., 2021).

Based on these findings, a baseband maximum likelihood estimator (MLE) using the S-LCI correlation from the FrFT output is proposed for a snapshot GNSS receiver that will overcome the drawbacks of the classic GNSS baseband architecture in harsh environments. This algorithm can jointly estimate the Doppler frequency and TOA of the incoming GNSS signal. By exploring the inner connections of code/carrier frequency and phases through the S-LCI correlation, this joint MLE is expected to increase the estimation accuracy of the Doppler effect (or carrier Doppler frequency) and TOA (or code phase) compared to the traditional FLL/PLL/DLL architectures (Progri et al., 2005). In summary, the main contributions of this work are as follows:

To improve the GNSS signal estimation accuracy in a weak environment, the S-LCI technique was used in the correlating process to increase the SNR, and the FrFT was used to estimate and compensate for the Doppler rate in the incoming signal;

The traditional FLL/PLL/DLL architectures require a relatively accurate Doppler frequency through acquisition to initialize the continuous tracking process to guarantee the accuracy and reliability of baseband signal processing despite frequent signal interruptions. Therefore, this work leverages the abundant geometry distribution of S-LCI correlations in the code domain (CD), frequency/Fourier domain (FD), and fractional Fourier domain (FrFD) to estimate the snapshot initial Doppler frequency error and initial carrier phase error followed by a relaxation of the initial guess of Doppler frequency;

After the initialization, the traditional FLL/PLL/DLL must take some time to converge. This can be unfriendly to a reliable and accurate navigating system, especially in challenging environments. To address this concern, this work adopts the gradient descent in the joint Doppler effect/TOA MLE to accelerate the convergence and obtain the super-resolution measurements (SRMs) of GNSS.

The paper will be organized as follows: Section 2 investigates the baseband signal processing algorithms for snapshot GNSS receivers as shown in Figure 1, where the main contributions are highlighted in red; Section 3 derives the Cramer-Rao lower bound (CRLB) of the MLE proposed in Section 2; numerical simulations and discussions are provided in Section 4; Section 5 realizes a prototype of the proposed snapshot receiver and assesses it with real-world GPS L5Q signals; and Section 6 concludes this paper and states future works.

## 2 BASEBAND SIGNAL PROCESSING FOR SNAPSHOT GNSS RECEIVERS

### 2.1 Dynamic GNSS Signal Models

A dynamic GNSS signal model is the prerequisite of this research. The dynamic signal model for the intermediate frequency (IF) GPS L1 C/A signal is modeled as shown in Equation (1):

1

with

where *Δa* is the amplitude error of the baseband signal; *t* is the time variable; *D*(*t*) is the navigation data; *C*(*t*) is the spreading code; *f _{I}* is the IF;

*ϕ*is the initial IF carrier phase;

_{I}*τ*is the propagation delay about the distance between the satellite and the user at the time of signal emission;

_{PRN}*τ*(

_{dyn}*t*) and

*ϕ*(

_{dyn}*t*) are the respective propagation delay residual and carrier phase residual determined by LOS dynamics and signal propagation time;

*v*is the dynamic term related to the velocity of the changing distance while

*β*is the counterpart related to acceleration;

*c*is the speed of light;

*f*is the initial Doppler frequency shift;

*f*is the radio frequency;

_{r}*μ*is the Doppler rate;

*ϕ*

_{0}is the initial carrier phase directly related to the user’s position; and

*ϕ*is the carrier phase caused by the frequency-mixing process.

_{LO}For the traditional GNSS baseband signal processing, the higher-order dynamic term β related to the Doppler rate in *τ _{dyn}* (

*t*) and

*ϕ*(

_{dyn}*t*) is usually not modeled. Therefore, the local signal cannot be adequately replicated as the real signal. Consequently, the longer the coherent integration interval, the more apparent the discrepancy between the real signal and the local replica will be.

### 2.2 Digital Partially matched filter (PMF)-FrFT technique

This work uses FrFT to compensate for the Doppler rate in the S-LCI correlator outputs. The practical approach to implementing fast FT is FFT; correspondingly, the counterpart for FrFT is digital FrFT (Ozaktas et al., 1996). This part will first introduce the FrFT. Then, a means to realize the digital FrFT in the GNSS baseband will be discussed.

Assuming that the signal *χ*(*t*) satisfies Dirichlet conditions, the FrFT for *χ*(*t*) can be defined as an integral transform as shown in Equation (2) (Almeida, 1994; Ozaktas & Kutay, 2001):

2

with

where *K _{p}*(

*t, u*) represents the transformation kernel of FrFT; ;

*p*is the FrFT order; and

*α*denotes a rotation angle from the time-frequency plane in terms of the ordinary FT to an extended plane and satisfies . The FrFT is a generalized form of FT and it is equivalent to the FT algorithm when

*p*is equal to 1 (Almeida, 1994). The diagram that compares FrFT and the ordinary FT for dynamic signals, also known as Chirp signals here, is shown in Figure 2.

While an ordinary FT transfers the dynamic signal from the time to the frequency domain, the signal power cannot be concentrated in this normal time-frequency plane. After FrFT, i.e., a rotation with an angle *α* for this plane, the dynamic signal power can be concentrated in the FrFD.

Digital FrFT for the GNSS baseband signal processing will be discussed. However, the traditional baseband algorithm is first explored as a comparison. As shown, *T _{t}* and

*N*denote the respective S-LCI interval and incoming IF samples over

_{t}*T*in the following description.

_{t}The work process of the correlators based on the digital FrFT in the GNSS baseband is illustrated in Figure 3. The partially matched filter (PMF) technique is used for the digital FrFT to reduce the computational burden, similar to the PMF-FFT technique (Stirling-Gallacher et al., 1996). The input for digital PMF-FrFT is the slow-time correlation sequence (Figure 3; larger unfilled purple circles).

As a result, three types of correlator outputs are formed.

First, a sequence of fast-time correlations (Figure 3; smaller unfilled purple circles) is formed in the time domain (TD) after integrating and dumping (I&D) of the incoming IF samples;

Second, a sequence of slow-time correlations is formed in the TD when the fast-time correlation sequence passes through the

*N*-point summation operators. This is a down-sampling process. One summation process represents one PMF. As usual,_{S}*N*is much smaller than_{S}*N*and they satisfy the condition_{t}*N*=_{t}*N*where_{S}P*P*and*N*are the IF samples within one PMF and the number of PMFs, respectively;_{S}Third, a sequence of S-LCI correlations is formed in the FrFD after the PMF-FrFT processes the slow-time correlator outputs. Thus, compared to the original digital FrFT with the computational load of

*O*(*N*_{t}log_{2}*N*) (Ozaktas et al., 1996), this calculation has been reduced to_{t}*O*(*N*_{S}log_{2}*N*) in the digital PMF-FrFT._{S}

In summary, different from the PMF-FFT, the correlation peak from the digital PMF-FrFT can be compensated by both the estimated frequency rate and shift errors (a two-dimensional [2D] match filtering process). By contrast, PMF-FFT can only compensate for the frequency shift error (i.e., a one-dimension [1D] match filtering process).

### 2.3 Closed-form models of the S-LCI correlations in the FrFD

This section considers the method that utilizes the abundant geometry distribution of the S-LCI correlator outputs in the CD, FD, and FrFD to relax the initial guess of Doppler frequency. The closed-form models for the S-LCI correlations in the first two domains are included in Appendix A; the models in the FrFD will be subsequently derived.

The expression for the sampled and digitalized *s(t*) in Equation (1) are given as shown in Equation (3):

3

where *n _{F}* is the index of incoming IF samples and

*T*is the IF sampling interval and

_{F}*s*[

*n*] will be used to indicate the sequence of

_{F}*s*(

*n*). Letting be the notation of the sequence of local replicas and

_{F}T_{F}*χ*

^{+}[

*n*,

_{S}*m*] be the sequence of correlator outputs in the time-code-frequency domain where

_{c}, m_{f}*n*is the index of the outputs, the correlator outputs are obtained by summing the outputs of the I&D operator between

_{S}*s*[

*n*] and ; this expression is known as the cross-ambiguity function (CAF) (Borio, 2008; Borio et al., 2009). Also,

_{F}*m*is the index of local code replica sequence about code chip offsets and

_{c}*m*is the index of local carrier replica sequence about Doppler frequency offsets. In this work, the outputs of I&D are called fast-time correlation sequences;

_{f}*χ*

^{+}[

*n*,

_{S}*m*] is nominated as a slow-time correlation sequence. These calculations are derived as shown in Appendix A and labeled in Figure 3.

_{c}, m_{f}According to previous research (Ozaktas et al., 1996; Stirling-Gallacher et al., 1996), the digital PMF-FrFT, i.e., implementing the digital FrFT to *χ*^{+}[*n _{S}*,

*m*] as shown in Equation (4) and Equation (5) results in:

_{c}, m_{f}4

for

with

5

where *m _{fp}* and

*m*represent the indexes of the match filters for frequency shift and FrFT order (related to the Doppler rate) in the FrFD, respectively;

_{p}*Δp*is the discrete search step for the FrFT order;

*M*,

_{c}*M*,

_{f}*M*, and

_{fp}*M*are the numbers of match filters in the CD, FD, FrFD (matching frequency shift), and FrFD (matching frequency rate).

_{p}The matched frequency shift and rate in the FrFD are then computed as shown in Equation (6) and Equation (7):

6

7

where is the matched rotation angle.

According to Moscardini et al. (2015), by applying the batches approximation, Equation (4) can be expressed analytically as shown in Equations (8–13):

8

with

9

10

11

12

13

where and are the respective estimated code delay in seconds, Doppler frequency in Hz, and carrier phase in radians from the NCO algorithms; *Δτ _{b}* and

*Δϕ*are the expressions of initial code phase error and initial carrier phase, respectively;

*Δf*[

*m*] is the expression of initial frequency error about the NCO frequency;

_{f}*f*is the code frequency; and

_{c}*f*(Δ

_{nC}*f*[

*m*],

_{f}*μ*) and

*ϕ*(Δ

_{n}*ϕ*, Δ

*f*[

*m*],

_{f}*μ*) are the closed-form expressions for the averaging frequency error and carrier phase over

*T*(Luo et al., 2018b; Ozaktas et al., 1996). The closed-form correlation amplitude shown in Equation (11) has been analyzed in detail in the authors’ previous publication (Luo et al., 2019b) and will not be considered further here.

_{t}In summary, implementing digital PMF-FrFT to the slow-time correlation sequence results in S-LCI correlation sequence variations with both matched frequency shift and rate in the FrFD. Thus, a 2D set of the S-LCI correlations in the FrFD is formed and can be written in matrix form as follows:

By utilizing the S-LCI correlation sequence in the CD, FD (see Appendix A), FrFD, and related closed-form models, an MLE will be subsequently investigated.

### 2.4 The Joint Doppler Effect/TOA MLE

This section explores the joint Doppler effect/TOA MLE based on the closed-form S-LCI correlations for GNSS baseband signal processing.

A diagram of the proposed baseband is shown in Figure 4 which includes the following work process:

First, 1000 slow-time correlation samples are generated with the unchanged carrier and code numerically-controlled oscillator (NCO) frequencies. By contrast, counterparts in the traditional baseband are derived from 1000 different NCO frequencies updated once at each tracking interval.

Second, five channels are used to process one satellite signal. In addition to the traditional early-, prompt-, and late-code delay channels, two extra channels (i.e., low- and high-frequency channels) are used.

Third, a pre-processing algorithm takes the place of the code and carrier discriminators to process the correlator outputs from the five channels.

Finally, a joint Doppler effect/TOA MLE is used instead of carrier/code loop filters to produce the final Doppler frequency error, carrier phase error, and code phase error as sources of GNSS measurements.

#### 2.4.1 Pre-processing algorithms

A diagram of the pre-processing algorithms in the baseband is shown in Figure 5. The pre-processing algorithm starts from the S-LCI correlation sequence in the prompt channel. Decision-making is then implemented to determine a set of frequency shift and rate indexes corresponding to the S-LCI correlation peak.

To implement decision-making, the search step for matching the code phase, frequency shift, and rate must first be determined. The measurements within a more optimistic decorrelated effect can help to generate non-singular and more accurate estimates (Li & Pahlavan, 2004). Thus, a “rule-of-thumb” early-late code and low-high frequency spacings have been adopted as 0.5 chips and 2(3*N _{S}T_{S}*)

^{−1}Hz, respectively; use of these search steps guarantees the lowest correlating power peak that can be tolerated. Similarly, the small-large rotation angle related to Doppler rate estimation is

*πΔ p*where

_{o}*Δp*is an optimum search step of an FrFT order (Luo et al., 2019b). This is shown in Equation (14):

_{o}Then, for , and , the decision-making process in the CD, FD, and FrFD is given by

14

where *DF ^{p}* {·} is the digital operator of the FrFT algorithm with respect to an FRFT order

*p*; and is the measured noisy

*χ*[

^{+}*n*,

_{S}*m*] in the baseband. The matching process is illustrated in Figure 6.

_{c}, m_{f}Next, and are conveyed to the early, late, low, and high channels. , are the measured S-LCI correlations from the prompt, early, late, low, and high channels, respectively. Their normalized amplitudes are computed as shown in Equations (15–19):

15

16

17

18

19

where *N _{F}* is the IF sample number over

*T*.

_{S}The S-LCI correlations are also related to small and large (rotation angles) channels generated by the outputs of the digital PMF-FrFT in the prompt channel. The related amplitudes are calculated as shown in Equation (20) and Equation (21):

20

21

Next, the carrier phase measurements for the respective normal, low, and high (frequency) channels are calculated through a four-quadrant arctangent/pure PLL discriminator as indicated in Equations (22–24):

22

23

24

where and represent the initial frequency error and Doppler rate estimates. A measured carrier-to-noise ratio density related to the incoming signal amplitude is also used in the MLE. In practice, the measured C/N_{0} can be obtained from previous-epoch measurements and predictions of empirical models (Groves et al., 2020), among other sources.

#### 2.4.2 System design of the proposed joint MLE

This section describes the system design of the proposed joint MLE that takes advantage of outputs from the pre-processing algorithms for GNSS baseband signal processing. A diagram of the proposed MLE algorithm is shown in Figure 7.

A linearized system model was built. The state vector is given by Equation (25):

25

with

where the first three states represent initial carrier phase error in cycle, code phase error in chip, and Doppler frequency error in Hz, respectively; *μ* is the Doppler rate state in Hz/s and *Δa* is the state of the normalized amplitude error as explained in Equation (1).

The closed-form vector representing the measurements of the MLE is defined by Equation (26):

26

with

27

where , is the noise variance of the in-phase/quadrature summations of the measured S-LCI correlation in a noise channel. As the S-LCI correlation sequence follows a Rayleigh distribution when there is no signal present in the signal processing channel (Kaplan & Hegarty, 2017), a noise channel is used to estimate the noise variance in this instance (Luo et al., 2019a).

As the unknown state vector is nonlinear with ** f**(

**), a linearized form should be derived to permit estimates to be updated iteratively. This is given by Equation (28):**

*θ*

with

28

where denotes the initial guess (or the previous-epoch estimation) of the state vector. The design matrix then satisfies Equation (29):

29

The measurement vector is then given by:

30

with

#### 2.4.3 Gradient descent optimization

Let represent an estimate of the state vector *θ* with the measurement vector . A loss function in terms of the S-LCI correlation amplitude (related to *Δτ _{b,c}*,

*Δf*,

*μ*), the carrier phase (related to

*Δϕ*), and the C/N

_{c}_{0}measurements (related to

*Δa*) is proposed as shown in Equation (31):

31

where * C_{s}* is the

*a priori*covariance matrix.

Assuming that the optimization problem is unbiased and measurement noise is additive white Gaussian noise (AWGN), the minimization of the loss function presented in Equation (31) can be made equivalent to maximize the linearized joint likelihood function as shown in Equation (32):

Maximize:

32

where det[·] is the determinant function; , and ** H** is defined by Equation (29).

Maximizing Equation (32) by gives the MLE as shown in Equation (33) (Kay, 1993):

33

where the expression *δ θ* is the vector that represents the state residuals of

*θ*as in Equation (28). The iterative estimation of the state vector is then given by Equation (34) as

34

#### 2.4.4 Carrier/code NCO frequencies and carrier phase/TOA updates

Two upgraded algorithms are proposed to synthesize the local carrier/code replicas in the proposed snapshot receiver.

First, both carrier and code NCO frequencies are compensated by the estimated Doppler rate and computed as shown in Equation (35) and Equation (36):

35

36

with

where subscript *n _{C}* is the epoch index for a working snapshot receiver; , and are the estimated Doppler frequency, the estimated frequency error, and the estimated Doppler rate, respectively and

*γ*is the coefficient for the Doppler rate estimation (introduced via an empirical value in this work).

Second, both carrier and code phases of the local replicas are also directly adjusted by the estimated carrier and code phase errors as follows in Equation (37) and Equation (38):

37

38

where is the estimated carrier phase error and and are the updated carrier phase and TOA estimates, respectively.

Of note, , and were all computed from Equation (34).

## 3 CRLBS OF THE PROPOSED MLE

This section investigates the CRLBs of the estimates for the first four states of ** θ**. The S-LCI correlation amplitude in the FrFD is observed in AWGN as indicated in Equation (39):

39

where is the matched correlation amplitude in the FrFD, while the unknown parameter about true amplitude and *w _{b,FrFD}* denotes the AWGN. The subscript

*FrFD*will be omitted in the subsequent derivation for simplicity. Based on the energy-preserving property, in the given condition, the noise variance will embrace the smallest value as shown in Equation (40):

40

subject to

Equation (40) is the prerequisite; then, the CRLB of the MLE can be approximated. At first, the Fisher information matrix (FIM) for ** x** is as shown in Equation (41) (Kay, 1993):

41

Let represent an estimate of the state vector *θ* with the measurement vector ** x**. Then, the mean squared error matrix of satisfies

42

Based on previous discussions, the minimum variances of the estimates of , and can be computed from the FIM as shown in Equation (43):

43

where the minimum covariance matrix of is derived with Equation (40). This is described further in Appendix B.

## 4 NUMERICAL SIMULATIONS

The numerical simulations that were carried out are described in this section. In an effort to progress toward the application of snapshot receivers, this section will include an assessment and discussion of the use of the proposed algorithm to overcome the flaws recognized in the traditional GNSS baseband. The simulating parameters for Monte Carlo experiments are listed in Table 1 where *g* denotes the standard gravity which is simulated as 9.80665 m / s^{2} in this study.

### 4.1 Improvements of the S-LCI Correlation in Dynamic Situations

To generate an intuitive illustration of the merits of FrFT for the S-LCI correlation in a dynamic situation, theoretical SNR curves based on discrete FT and FrFT values were plotted as shown for a pure carrier simulation (Figure 8). The received carrier signal power was set at −160 dBW and the environmental temperature at 290 °K. After I&D and the correlating process, the SNR curves in the FrFD are seldom affected by the incoming Doppler rates and tend to increase with the enlarging integration interval. However, this results in an FD decrease when the coherent integration interval exceeds a specific threshold.

To show how the proposed MLE improves the baseband signal estimation in a dynamic situation, the following two Monte Carlo experiments were implemented.

The findings shown in Figure 9 offer 2D-estimate plots of the Doppler rate and frequency error estimates as the input Doppler rates increase. The digital PMF-FrFT and the proposed MLE results present outcomes that are close to true values. By contrast, the traditional PMF-FFT gradually deviates from the truth as input dynamics increase.

The root-mean-square error (RMSE) (2σ) curves of frequency error estimates from the Monte Carlo experiment were then computed (Figure 10). The accuracies based on the PMF-FrFT and the proposed MLE do not drop as input Doppler rates increase. However, the accuracy of the traditional PMF-FFT is dramatically reduced when the Doppler rate increases. Thus, the results indicate that the traditional PMF-FFT fails to match the input Doppler rate.

The results of these two simulations verify that the PMF-FrFT technique efficiently improves the S-LCI correlation performance in a dynamic situation compared to the traditional PMF-FFT technique used in an ordinary receiver.

### 4.2 Relaxing Initial Frequency Error and Fast Convergence

The proposed algorithm allows the GNSS baseband to tolerate a more significant initial frequency error than the traditional PLL/DLL when implementing a fast convergence. As a comparison, the traditional tracking algorithm is simulated to estimate a 1.024-s digital GPS L1 C/A IF sequence.

As shown in Figure 11, PLL and I/Q correlator outputs are compared when different initial frequency errors and Doppler rates are added to the simulations. When the values are small (up to 1 Hz and 0 g), the tracking loop converges within a short time and provides clean demodulated I/Q correlations over almost the full 1.024-s processing time. Nevertheless, as the initial frequency error and Doppler rate increase, the converging time also increases. Finally, the signal parameters are difficult to estimate within the same processing interval at the initial frequency error of 60 Hz and the Doppler rate of 0.5 g.

To show that the proposed MLE can tolerate a larger initial frequency error with fast convergence, the RMSEs (2*σ*) of Doppler rate and frequency error estimates based on 100-trial Monte Carlo experiments are calculated as initial input frequency errors increase. The input Doppler rate is simulated as 0.5 g.

The results illustrated in Figure 12 demonstrate that the estimation accuracies of the proposed algorithm are maintained at a high level within the maximum frequency error of approximately 200 Hz with RMSEs at approximately 0.04 Hz and 0.05 Hz/s for the frequency and Doppler rate estimates, respectively. However, for traditional tracking, even if the initial frequency error and Doppler rates are below 1 Hz and 0 g, respectively, the standard deviation (STD) (i.e., the 1σ uncertainty that can be read directly and estimated from the PLL curves shown in the first panel of Figure 11) of the Doppler error estimation as high as 0.83 Hz. This is 20 times worse than what might be achieved with the proposed method. Even worse, the traditional tracking loop fails to converge at a much lower frequency of 60 Hz.

In summary, tolerance of the initial frequency error estimate and the convergence time in the baseband signal processing have been significantly improved by the proposed MLE algorithm.

### 4.3 Super-Resolution Performance of the Doppler Effect and TOA Estimation

To examine the super-resolution performance of the S-LCI-based algorithms for the proposed snapshot receiver, the RMSEs of the Doppler rate as well as the frequency, carrier phase, and code phase error estimates from the Monte Carlo experiment were compared with the respective CRLBs in the three-dimensional (3D) plots shown in Figure 13. Related discussions are provided as follows.

First, compared with the ordinary estimation accuracy of the Doppler rate (denoted as PMF-FrFT), the frequency error (denoted as PMF-FrFT), and the carrier phase error (see Equation (22) and denoted as PMF-FrFT+pure PLL disc), the proposed PMF-FrFT MLE results exhibit superior performance.

It is also worth noting that a sharp drop occurs in the MLE accuracy below 20-dB-Hz. Shown in Figure 14 is a plot of the theoretical detection probability curves for the digital PMF-FrFT process where *T _{t}* = 1.024 s,

*T*= 0.25 ms and “Pfa” denotes the false alarm rate. The curves include the power loss behind the frequency shift and Doppler rate match filter resolution. For example, a simulation in which the Doppler rate error and the worst frequency error are ignored can be approximated as (4

_{S}*T*)

_{t}^{−1}≈ 0.244 Hz in this case. In their previous work, the authors proved that the theoretical results were in accordance with the Monte Carlo simulations (Luo et al., 2018b). As is apparent, the curves show that the detection probability drops from 1 when C / N

_{0}becomes lower than 15 dB-Hz (nb: the existing Doppler rate error increases this value from a practical standpoint). In other words, random noise power gradually exceeds the correlation power peak as C / N

_{0}decreases. Considering the measurements in the proposed MLE from Equations (15)–(24), the estimate accuracy may be highly constrained. Moreover, there is no compensation for the long correlation power loss resulting from code grid migration error at the current stage of this work (Yang et al., 2022) or the code Doppler rate, which are two factors that could affect the sensitivity. These factors will be considered in our future work.

When the input C/N_{0} is below 30 but above 19 dB-Hz, the estimated results based on the proposed method reach the CRLB.

Furthermore, as the C/N_{0} increases beyond 30 dB-Hz, the discrepancy between the estimated accuracy and the CRLB becomes larger. This is because the bandwidth of the PMF-FrFT (other than the thermal noise bandwidth) gradually dominates the estimation accuracy. More specifically, the measured long correlation power would be restricted by the slow-time interval of the PMF-FrFT due to the residual code phase error estimation. Furthermore, as mentioned earlier, there is no compensation for code migration or the code Doppler rate in this work. These factors will induce extra noise that is not modeled in Equation (39). As the input C/N_{0} increases, the unmodeled errors cause the proposed MLE to be slightly biased. This issue will be solved in our future work.

Before discussing code phase errors, it is important to recognize that the traditional code error discriminator, i.e., the early minus late (E-L) envelope discriminator (Kaplan & Hegarty, 2017), is used for S-LCI correlations to discriminate the code phase error as a comparison in simulations, i.e.,

44

where the result from this traditional code discriminator is denoted as PMF-FrFT+ EL disc.

For the code error estimates, the accuracy of the proposed PMF-FrFT MLE performs slightly worse than the traditional PMF-FrFT+ EL disc. This may be because:

the code has a much narrower bandwidth than the carrier. For example, the carrier frequency of the GPS L1 C/A signal is 1540 times higher than the code frequency while the ratio is 115 times higher than the GPS L5 signal. Therefore, dynamic stress has much less influence on the code than the carrier. The proposed PMF-FrFT MLE improves the code error estimation much less effectively than the carrier. Of even greater concern, error propagations from the other states in the MLE reduce the code accuracy.

the results of the traditional PMF-FrFT+ EL disc method are almost identical to the CRLB (even at a low C / N

_{0}of 15 dB-Hz). The traditional E-L discriminator produces a true tracking error within its linear region (Kaplan & Hegarty, 2017); thus, the result is consistent with the theoretical analysis.

Next, to explain how the super-resolution performance will be achieved using the proposed S-LCI technique, the code tracking error and sensitivity of the traditional PLL/DLL are first simulated for comparison purposes and cases with coherent integration times of 1 ms and 20 ms were both tested. Code error accuracy curves based on both the proposed MLE and the traditional tracking algorithms were plotted as shown in Figure 15.

At first, the traditional tracking loop starts to be out of lock at 25 dB-Hz and 18 dB-Hz based on the respective 1-ms and 20-ms coherent integration times. Also, the two traditional methods exhibit substantial reductions in code accuracy compared to the three S-LCI-based algorithms.

The code accuracy of S-LCI-based E-L discriminators, i.e., PMF-FRFT+EL disc and PMF-FFT+EL disc, remains consistent with the CRLB curve over all the simulated C / N_{0} cases. For example, when the C / N_{0} is low at 20 dB-Hz compared with the traditional 20-ms algorithm (21 m), the two S-LCI-based E-L discriminators (11 m) and the proposed S-LCI-based MLE (17 m) exhibit improved code accuracies of 47.6% and 19.0%, respectively.

In summary, the accuracy and sensitivity of code phase (or TOA) estimation can be significantly improved in snapshot receivers by using S-LCI signal processing.

### 4.4 Comparisons of the Computational Load

In this part, we compare the computational loads of the proposed algorithms with those of the traditional tracking loop to provide an understanding of algorithm efficiency.

Then, Table 2 lists the terms used to compare the two algorithms. The characteristics related to the proposed baseband architecture refer to Figure 5; *N _{F}* is the IF sample number over

*T*and

_{S}*N*is the number of points of the digital PMF-FrFT as discussed earlier.

_{S}The input Doppler rate is 0 g and the RMSEs of the S-LCI-based algorithms were computed via Monte Carlo 100-trial experiments with 2*σ* uncertainty. Results related to the traditional algorithm were obtained from the tracking over one processing time. Simulation parameters for the proposed MLE include C/N_{0} is 45 dB-Hz; *T _{t}* is 1.024 s;

*T*is 0.25 ms; and

_{S}*Δp*is 0.00055. Simulation parameters for the traditional algorithm including tracking times of 1.024 s and 1.020 s when coherent integration intervals were 1 ms and 20 ms, respectively, and bandwidths of the 3rd-order PLL and 2nd-order DLL that were 18 Hz and 0.1 Hz respectively. The receiver used early-late envelope code and pure PLL discriminators; early-late correlator spacing was 0.5 chips.

_{o}The complexity related to one-time digital FrFT *O*(*N _{S}* log

_{2}

*N*

_{S}) was presented in an earlier publication (Ozaktas et al., 1996). A more detailed version of this point was also presented in the authors’ previous work (Luo et al., 2018b).

## 5 PERFORMANCE OF THE SNAPSHOT RECEIVER WITH WEAK GPS L5Q SIGNALS IN A REAL-WORLD EXPERIMENT

The static real-world GPS L5 IF data are collected with LabSat 3 Wideband in both open-sky (the peak of Mount High West) and indoor (Fortune Metropolis in Hung Hom) areas of Hong Kong. The IF sampling rate is 58 MHz and the local clock is based on the high-performance oven-controlled crystal oscillator (OCXO). This work focuses specifically on theoretical analysis.

Because it is much more difficult to model clock error behaviors, they are not considered in the proposed algorithm. Given this circumstance, we use the high-quality oscillator (i.e., OCXO) instead of the low-cost temperature compensated crystal oscillator (TCXO) to avoid the clock-vibration-induced inconsistency as much as possible between the theoretical analysis and the real-world experimental results. The applications based on the TCXO will be completed in our future work.

Figure 16 shows the experimental set-up. Views of the surrounding environment and the sky plots are shown in Figure 17 with the displayed C / N_{0} values obtained from the U-Blox F9P receiver. The F9P device used the same antenna as the LabSat 3 Wideband through a splitter. The pilot channel of the collected IF L5 signals (L5Q) was used in these experiments as the known secondary codes facilitated the implementation of the S-LCI without external aiding (Hegarty, 2006). The parameter settings for the digital PMF-FrFT used in the proposed snapshot receiver are listed in Table 3.

Before the experiment, initial Doppler frequency and code phase were acquired based on the parallel-code search (PCS) acquisition algorithm (Hegarty et al., 2003). The initial Doppler frequencies obtained for SV18 and SV6 of the open sky were 2366 Hz and −2081 Hz, respectively, and 1033 Hz for SV32 for the indoor environment.

To show how the FrFT technique improves the baseband S-LCI correlation in practice, 3D views of the normalized correlation amplitudes of CAFs were compared using ordinary FD and FrFD, as shown in Figure 18. The S-LCI correlation peak shown in the FrFD panel clearly has a higher resolution. As shown, its power peak is higher and more concentrated than in the case of ordinary FD; it is also more capable of mitigating side lobe interference. In this research, the GNSS measurements of the snapshot receiver using the S-LCI technique were nominated as the SRMs.

The GNSS measurements, excluding the outliers, were assessed quantitatively via computations of RMSEs, STDs, and mean values of the signal parameter estimates from the proposed snapshot receiver as listed in Table 5 as well as RMSEs using a traditional receiver with a 20-ms coherent integration time as listed in Table 4. In the analysis, the Doppler rate is assumed to be unchanged over the full 200-s processing time with reference to the true Doppler rates obtained from the digital PMF-FrFT algorithm over a 10-s S-LCI interval. As shown in Figure 19, the Doppler rate references are matched as −0.1407 Hz/s and 0.0704 Hz/s for SV6 and SV18, respectively. The results of the Doppler rate estimates are plotted in Figure 20. Clearly, the estimated frequencies are significantly de-noised by the proposed snapshot receiver.

Several conclusions can be drawn based on the statistical results. At first, as expected, the accuracy of signal estimation increases typically as the S-LCI interval gains for SV6. The exception is that the RMSE of code phase error estimation decreases slightly when the S-LCI interval is up to 5s. The code phase migration due to the code rate error between the received and locally replicated signals in the long coherent correlation may be the source of this unexpected phenomenon. More specifically, the code migration enables the correlation power to deviate gradually from the linear region in the integration process as the coherent processing interval increases. Likewise, the final S-LCI correlation cannot be represented precisely by the form as modeled in Equation (8). Collectively, this will result in a more biased code estimation in the 5-s-based approach.

Next, a worse result based on the 5-s S-LCI approach occurs in the estimation process of SV18, where the carrier phase accuracy is also very low. It can be explained that the real Doppler rate in this case is closer to the non-linear region of the numerical model of the digital FrFT match filter. The proposed MLE is built based on the linearized closed-form expressions of ** f**(

**) as shown in Equation (29) and leads to degraded Doppler rate estimation performance. Again, the code and carrier frequency/phases are estimated based on the Doppler rate estimate according to Equations (35) to (38). As a result, this will lead to abnormal accuracy.**

*θ*Finally, compared with the results from the traditional receiver, the proposed snapshot receiver further accelerates the TOA and Doppler effect-dependent estimation accuracy. For instance, based on the RMSEs from SV6, the maximum improvements in the code and frequency error accuracy are 48.3% and 96.8%, respectively.

The results shown in Figure 15 document that the simulated code phase errors based on the early-late amplitude discriminator reach the CRLB model provided. In addition, it was also critical to determine whether the theoretical model could describe real-world signal behavior.

First, we know that the evaluated real-world signals came from an open-sky environment; this excludes contributions from multipath effects and the other unintentional interference factors in the statistical analysis. Second, we used the embedded high-quality OCXO to collect the IF data so that most of the irregular dynamic errors caused by the oscillator vibration are removed. Overall, we have tried to provide statistical results from the purest possible real-world environment to satisfy the condition for the theoretical model.

In this case, we compare the CRLB curves to the RMSEs of the L5Q (BPSK[10]) code error in Figure 21. As shown, the code accuracies from the 0.5-s and 1-s S-LCI methods comply nearly completely with the CRLB curve except for the one at 5-s. The reason for the exception of the 5-s method was explained earlier. Therefore, our findings verify the effectiveness of the proposed snapshot receiver for real-world signals.

Next, the findings shown in Figure 22 are a comparison of the results of the response of the proposed 1-s S-LCI snapshot baseband MLE with the traditional tracking algorithms to the large frequency error. The tested real-world IF data came from the channel of SV6. The top panel proves that the proposed algorithm initiates a prompt response to the 51-Hz error and that it converges to the original estimate level (without manually-added frequency error) within one snapshot epoch. By contrast, the upper bound of tolerance of frequency error of the traditional tracking loop was tested as well. As indicated by the results shown at the bottom of the Figure, the traditional algorithm can converge within a maximum of 20 Hz while getting out of lock when a value of 1-Hz or larger value is added. Compared to the traditional algorithm, the robustness of the proposed algorithm was significantly enhanced as the proposed baseband MLE tolerates an incoming frequency error that is 2.5- fold higher than that tolerated by conventional PLL/DLL. Meanwhile, the proposed algorithm produced more accurate frequency error estimate results in the experiment.

We also assessed the indoor environment where the signal strength from the commercial U-Blox F9P receiver was as shown in Figure 23. The Doppler frequency estimates are provided in Figure 24. The start of the processing time was not strictly aligned with the GPS time in these two figures, as several-second clock biases occur between these two examples. However, Figure 23 contains the entire processing time of Figure 24. A data set that highlights the challenges associated with indoor signal reception would be sufficient in this case.

As shown in Figure 23, the signal strength becomes lower than 30 dB-Hz during the first 100 epochs and drops to approximately 20 dB-Hz for several seconds. During these periods, the commercial receiver is unable to process normal output range measurements.

By contrast, the estimates from the proposed receiver algorithm are much more precise than those generated by the traditional method. However, one noticeable difference from the open-sky result is that the fitted Doppler frequency does not generate a highly accurate reflection of the Doppler estimates. This can be explained by the fact that the severe indoor multipath effect can generate dramatic changes in the Doppler rates. Thus, the matched Doppler rate from the first 10-s signals (it is estimated as 0.2287 Hz/s) does not properly represent the values for the entire evaluation time.

Finally, we compute the RMSE for the code error estimated from the 1-s S-LCI approach as 0.679 m. The related DLL output accuracy from the traditional receiver algorithm (see the configurations provided by the caption of Table 4) is 1.099 m. This improves the code accuracy by 38.2%. Thus, as expected, the proposed snapshot receiver also performs substantially better than the traditional receiver in indoor situations.

## 6 CONCLUSIONS AND FUTURE WORK

This paper proposes the creation of a GNSS baseband MLE using the S-LCI correlation in the FrFD to generate a snapshot receiver. The MLE can jointly process the Doppler effect and the TOA in the baseband. CRLBs were also derived for the proposed MLE.

Compared with the traditional baseband processing algorithms, numerical simulations verified the superiority of the proposed snapshot receiver as follows:

At first, the proposed receiver produces more accurate signal estimates than the traditional algorithms in weak and dynamic conditions which were close to the CRLB;

The proposed snapshot baseband tolerated a 200 Hz initial frequency error and converged within a short 1-s-level coherent processing interval compared to the maximum of 60 Hz of the traditional tracking loop.

However, to evaluate the consistency of these practical results with the theoretical analysis as well as their feasibility in commercial markets, the proposed snapshot receiver prototype was tested in both real-world open-sky and indoor environments (i.e., experiments for the satellites with low elevation angles). The experimental results were consistent with the simulations, including:

The signal estimation accuracy was significantly improved in both open-sky and indoor environments in which the code accuracy was increased by a maximum of 48.3% and 38.2% for SV6 and SV32, respectively;

The experiments verified that the code accuracy was close to the CRLB in the weak open-sky real-world environment;

The proposed snapshot receiver not only promptly converged from an incoming frequency error of 51 Hz (which is 2.5-fold better than the traditional algorithm) but also produced a more accurate frequency error estimate.

In the future, an improved version that considers the multipath signal model will be integrated into the proposed PMF-FrFT MLE.

The proposed technique has a broad scope of application as it overcomes many of the challenges faced by existing GNSS receivers. For example, it achieves high accuracy in weak and dynamic situations, is less likely to suffer from signal interruption, and features signal processing processes that converge rapidly. Besides, a looser initial guess was proved to achieve high-accuracy signal estimation even in challenging environments. Thus, this advance may provide crucial support for future commercial markets, including autonomous driving, urban/indoor navigation and positioning, and the Internet of Things, among other ventures.

## HOW TO CITE THIS ARTICLE

Luo, Y., Hsu, L., & El-Sheimy, N. (2023). A baseband MLE for snapshot GNSS receiver using super-long-coherent correlation in a fractional fourier domain. *NAVIGATION, 70*(3). https://doi.org/10.33012/navi.588

## ACKNOWLEDGMENTS

This research is supported by the University Grants Committee of Hong Kong under the scheme Research Impact Fund for project R5009-21 “Reliable Multiagent Collaborative Global Navigation Satellite System Positioning for Intelligent Transportation Systems.” This research was also supported by funding provided to Professor Naser El-Sheimy from NSERC CREATE and Canada Research Chairs programs. The authors are grateful to Mr. Chin Lok Tsang at The Hong Kong Polytechnic University for collecting the GPS IF data. The authors would also like to thank Dr. Pai Wang from Shanghai Jiao Tong University for careful reading, professional discussions, and valuable suggestions for the manuscript and the anonymous reviewers for their constructive comments to improve the manuscript.

## APPENDIX A THE S-LCI CORRELATIONS IN THE CD AND FD

According to Equation (3), the slow-time correlator output is computed from Equation (45):

45

with

where this process is commonly nominated as I&D to the GNSS baseband processing; when *P* is equal to 1, *χ ^{+}* [

*n*,

_{S}*m*] becomes the fast-time correlation sequence;

_{c}, m_{f}*d*is the spacing of two adjacent match filters in the CD in chips; and

_{c}*Δf*is the corresponding spacing in the FD. Then, its closed-form expression is given by Equation (46):

_{nco}46

with

where *R _{cc}* [·] is the operator of the correlation function of two spreading code sequences;

*Δτ*[

*n*] is the code phase error sequence in the slow-time domain;

_{S}*χ*[

*n*] is the slow-time correlation sequence model between the incoming and local replicated carrier signals; and

_{S}, m_{f}*ϕ*[

*n*] is the carrier phase sequence. The closed-form expression for

_{S}*χ*[

^{+}*n*,

_{S}*m*] has been discussed in the authors’ previous publication (Luo et al., 2019b), and will not be discussed further in this manuscript.

_{c}, m_{f}Next, the S-LCI correlation in the traditional baseband is computed from Equation (47):

47

where *n _{F}* = 0 and subscript

*n*is the index of the S-LCI correlation sequence.

_{C}To summarize, is the set of the S-LCI correlations in the CD; is the counterpart in the FD; is the set of slow-time correlations in the TD.

## APPENDIX B DERIVATIONS OF THE MINIMUM COVARIANCE MATRIX

Firstly, the minimum variance for the amplitude measurement in the FrFD will be introduced. The likelihood function of the amplitude measurements from Equations (15) to (19) leads to Equation (48):

48

Then, the negative of the second derivative of Equation (48) gives the Fisher information of the estimation as Equation (49):

49

The minimum variance of the estimation satisfies Equation (50) (Kay, 1993):

50

where, according to Equation (27), this can be computed from

51

with *T _{t}* =

*N*and

_{S}T_{S}*a*= 2(1 +

*Δa*), where

*a*is the amplitude of the IF signal. Thus, the minimum variance model is as per Equation (52):

52

According to Equation (9), the minimum variance of can be computed by the respective minimum variances of and .

Then, based on the closed-form Equation (8) of , linearizing about *Δf*[*m _{f}*] and

*μ*gives the approximation as

53

When the observation number increases, the PDF of becomes more concentrated around the mean value of its truth. Thus, the observed values of lie in a small interval satisfies . Over small intervals, the nonlinear transformation can be approximately linear as noted above. Then, based on this model, the closed-form *μ* and *Δf*[*m _{f}*] over the small interval can be approximated as and , and their derivations referring to the equation above are given as follows:

where is the constant of true amplitude here. Using a transformation of the parameters, the CRLB of Doppler rate estimation can be derived as:

Similarly, the CRLB of frequency error estimation can be derived as:

The closed-form expression of the carrier phase error is a function of the Doppler rate variable *μ* as follows:

Thus, the CRLB of carrier phase error estimation is derived as:

The variance of the matched carrier phase estimates in the FrFD is then calculated as:

Finally, in this work, the minimum value of the *a priori* covariance matrix for the MLE satisfies the following equation:

where ** I_{i}** is an

*i*th-dimension unit matrix;

*is the zeros matrix with*

**O**_{m×n}*m*th row and

*n*th column.

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.