Modified Costas Loop for Carrier Phase Tracking in GPS Receivers

. The carrier phase received at the receivers of the Global Positioning System (GPS) links is used to detect navigation data and to precisely determine the position, speed and time corresponding to the user's equipment. Therefore, subsystems for carrier phase tracking are crucial parts in all GPS receivers. When the propagation conditions are favorable, the method frequently used for phase tracking is based on Digital Phase-Locked Loop (DPLL)) and implemented through the discrete Costas loop operating under the modulated L1 carrier, in the case of a GPS receiver. This technique is quite simple, well known and very suitable for implementation in low-cost receivers. In this article, we revisit the traditional Costas loop design and point out some issues that affect the phase tracking performance of this loop. In order to overcome these problems, we propose some modifications to the traditional Costas loop. The resulting architecture presents better performance and complexity equivalent to the original loop. Another contribution of this work is the mathematical analysis to evaluate the performance of the new architecture when operating on an Additive White Gaussian Noise (AWGN) channel. Various results from computational simulations carried out with the two architectures, in different operating scenarios, including AWGN, dynamic stress and ionospheric scintillation are presented and discussed. We conclude that the new architecture outperforms the traditional Costas loop in terms of the variance of the estimated phase error, root mean squared error of the estimated phase and robustness to cycle-slip and loss of lock.


Introduction
The satellites of the constellation of the Global Positioning System (GPS) transmit to the receivers, located on or close to the ground, carriers in frequencies in the L band that are processed to determine the geographical position of user's equipment [1].In this process, a crucial operation is the determination of the frequency and phase of the received carrier that are used to detect the navigation data and to solve the position, velocity, and time (PVT) in the GPS receiver [2].
The phase and frequency offset between the oscillators of the transmitter on the satellite and the receiver on the ground, the time delay due to the propagation of the signal and the relative speed between transmitter and receiver cause phase deviation and frequency deviation (Doppler shift) in the received signal which need to be estimated at the receiver.In addition to these effects, there are situations in which other phenomena such as multipath, shadowing, and ionospheric scintillation can strongly affect the phase of the received signal [3].
Usually, frequency and phase synchronization are done in two steps.The first step is called acquisition and in it a rough estimate of the frequency offset is made.The refinement of this estimate and the determination of the phase offset occur in the next step called tracking.Because of the time variability of frequency and phase offsets due to the dynamics of satellite orbits, among other factors, the tracking step needs to be maintained while the receiver is receiving the transmitted signal.
The method traditionally used for tracking the carrier phase is the DPLL (Digital Phase Locked Loop) whose performance is quite satisfactory when the propagation conditions are favorable, i.e., when the ratio of carrier power to noise density (C/N0) (carrier-to-noise density ratio) is high and ionospheric scintillation and other distortions caused by propagation are absent.On the other hand, in scenarios with severe scintillation or low C/N 0 , tracking with classic DPLL becomes inadequate, with the occurrence of cycle-slip or even complete loss of lock [4].Several new and more robust techniques for carrier phase tracking in GPS receivers, such as FLL-Assisted PLL, Kalman filter based PLL and Adaptive PLL, have been proposed in the literature in recent years [2], [3], [5][6][7].
The most robust architectures for carrier phase tracking in GNSS receivers designed to work precisely even in harsh propagation scenarios, with large dynamic stress and severe ionospheric scintillation, are based on extended Kalman filter whose associated state matrix incorporates autoregressive (AR) models for the magnitude and scintil-lation phase [8].Usually, these robust and high precision carrier tracking loops are designed to work from pilot carriers, without modulation.In paper [9], the parameters of AR models are determined and adjusted online, allowing the phase tracking system to work well even with varying scintillation intensity.The phase tracking loop proposed in [6] is based on the operation of two Kalman filters: a filter dedicated to carrier phase tracking and another for estimating AR process parameters.Recently, the authors of [10] presented a new Kalman filter-based phase tracking architecture in which the parameters of AR models are estimated through radial base function networks.One thing that all these architectures have in common is the high computational complexity required.
On the other hand, there are many applications intended for favorable propagation conditions where the utilization of low-cost and low-complexity GPS receivers is the most suitable choice.In these cases, the phase tracking may be done from the legacy L1 modulated carrier of the GPS system and the traditional DPLL implementing the Costas loop is the solution normally adopted [1], [4], [11][12][13].Focusing on this type of application, this work revisits the traditional Costas loop project and points out some limitations of this architecture.Aiming to overcome these limitations, it proposes a new DPLL architecture for carrier phase tracking in a GPS receiver intended for applications in the same scenarios where the traditional Costas loop is recommended.The solution presented consists of some changes in the original Costas loop that result in a significant improvement in the system's performance.On the other hand, the complexity of the proposed architecture is practically equivalent to the traditional Costas loop.
The main contributions of this work are: 1. Discussion on some points of the discrete implementation of Costas loop that reduces its performance; 2. Proposal for a new architecture for tracking the carrier phase; 3. Performance evaluation of the proposed tracking loop via mathematical analysis and computational simulations.
The discrete Costas loop for GPS receivers found in the literature [1], [4], [11][12][13] receives the passband signal from IF (Intermediate Frequency) receiver stage and converts it to baseband through an integrates-and-dump accumulator (I&D) in each arm of the loop.Therefore, the sampling frequency at which the loop filter operates depends on the I&D accumulation period.On the other hand, in the proposed architecture, as will be detailed in Sec. 3, the conversion to baseband and the correlation/accumulation carried out by I&D occur outside the tracking loop and, in addition, the loop sampling frequency is independent of the accumulation period.These simple modifications have a significant impact on system performance.This article is organized as follows: In Sec. 2, to establish a reference, we describe the traditional Costas loop.The model of the proposed architecture, including the details of the signals and functional blocks, as well as math-ematical analysis for evaluating the performance of the new loop are presented in Sec. 3. In Sec. 4, the results of the computational simulations that were carried out with the objective of evaluating the performance of the two architectures for different operating conditions are presented.In Sec. 5, the article is concluded.

Costas Loop Description
The block diagram of the traditional Costas loop is illustrated in Fig. 1 [1], [4], [12][13][14].The received signal, after being converted into an intermediate frequency (IF), is discretized, resulting in the signal x [n].Therefore, x[n] is composed of L1 modulated carrier added with white Gaussian noise and may be described by the following equation: where  [4] corresponding to the receiver thermal noise.
The purpose of the Costas loop is to estimate the carrier phase, shown in Fig. 1   The discriminator used in the Costas loop must be insensitive to phase jumps resulting from carrier modulation.
The following examples are suitable for this purpose: Of these, the DD has the smallest phase error variance when the carrier is distorted with only Gaussian noise because of the absence of quadratic loss [4], [12].
The loop filters commonly used in the Costas loop for carrier phase recovery in GPS receivers are first or second order whose transfer functions are expressed by [3] ( ) ( ) ( ) where , and The constant ω n is called the natural frequency of the loop.Typical values for choosing for constants a, b and c are √2, 1.1 and 2.4, respectively [1].When the loop filter is first order, because of the NCO, the loop will be second order.Second-order loop filter results in third-order loop.
Considering that second-order loop presents non-zero stationary error for phase tracking when there is an acceleration offset (Doppler rate) in the phase of the received signal, it is recommended to use third-order loop when receivers are subjected to high dynamic stress [1].Higher order loops are susceptible to instability.
The NCO is implemented by a discrete integrator whose transfer function is ( ) and is followed by the carrier generator block responsible for generating the signals cos During tracking, the phase error is much smaller than  , which allows linearization of the entire loop, resulting in the following transfer function [4] for the case of the second-order loop: In this equation, G(z) is given by ( ) and models the effect of I&D blocks inserted into the loop.
From (5), it is evident that G(z) contributes with the position of the poles and zeros of H(z).Thus, the operation of the loop is affected by G(z), mainly when 1/T a is comparable to the equivalent noise bandwidth of DPLL [4], denoted by B LD .However, when 1/T a >> B LD the effect of I&D is negligible, and the transfer function given in ( 5) is simplified to The value of B LD is determined by [16] ( ) ( ) When the sampling frequency at which the DPLL operates is high, that is, B LD T a << 1, the value of B LD and the equivalent noise band of the equivalent analog PLL, represented by B LA , practically coincide.However, as the sampling frequency 1/T a is reduced, the noise equivalent bandwidth of the DPLL is increased, that is, we have B LD > B LA .Table 1 illustrates the effects of sampling frequency and I&D employment on DPLL noise bandwidth.Assuming that the equivalent analog PLL noise bandwidth is B LA = 10 Hz, B LD is determined for H(z) expressed by ( 5) and (7).It is observed that a reduction of 1/T a and the use of the I&D block within the loop cause an increase in B LD .
One of the parameters used to define the performance of the Costas loop is the estimated phase error variance (σ φ 2 ).For the Costas loop working with a conventional discriminator (CC), this variance is given by [17] where C is the carrier power and N 0 is the one-sided power density of the Gaussian noise at receiver.The parenthetical term in (9) represents the quadratic loss due to the CC discriminator.It is interesting to note that for a given loop filter, an increase in T a reduces the quadratic loss, but on the other hand, it can increase B LD .
Another important characteristic for phase tracking systems is their robustness against the occurrence of cycleslip.The cycle-slip event for the Costas loop is characterized by jumps in the recovered phase, with magnitude equal to nπ rad.This phenomenon is caused by the nonlinear behavior of the loop when C/N 0 is too low or the received signal is distorted by ionospheric scintillation or other type of fading [4], [18].The less susceptible to cycleslip more reliable is the phase tracking system.

Modified Carrier Phase Tracking
This section presents a new architecture for carrier phase tracking in GPS receivers.The differences between this solution and the traditional Costas loop, presented in Sec. 2, are: 1.It works from the baseband signal; 2. The integrated-and-dump (I&D) blocks were removed from within the DPLL loop.This alteration eliminates losses resulting from the use of I&D within the loop; 3. The sampling frequency at which the DPLL operates is independent of the integration time Ta and the adopted value is such that B LD ≅ B LA .Thus, for a given selected loop filter, the DPLL noise bandwidth B LD is kept close to its minimum value.

Architecture Description
The block diagram of the proposed architecture is shown in Fig. 2. The signal r[n] at the input of the architecture represents the received discrete baseband signal and is at the sampling frequency 1/T s .Assuming the receiver is operating in tracking, when the bit timing and the C/A code are already known, the signal r[n] can be modeled by (10)   where A[n], d[n] and θ [n] are carrier amplitude, binary symbol and carrier phase as defined in (1).The signal z[n] is a discrete, complex, Gaussian, and white noise with zero mean and variance 2N 0 /T s corresponding to the receiver thermal noise.
The integrate-and-dump block accumulates N samples of r[n] as the signal is received and outputs the accumulated value at each period T a = NT s .Therefore, the sampling frequency at the I&D output is equal to 1/T a = 1/(NT s ).As it was considered that the bit synchronization is perfect, the N samples added will always belong to the same bit.The signal of the I&D output can be approximated by where samples of z[n] are independent, it turns out that the noise z D [n] is also complex, Gaussian, and white with zero mean and variance 2NN 0 / T s .
Aiming for the DPLL to operate at a single sampling frequency, independent of the value of T a and high enough for B LD ≅ B LA , the signal r D [i] is filtered by a discrete rectangular filter whose impulsive response is [ ] Consequently, the output signal of the rectangular filter is expressed by and will have a new sampling frequency equal to M/(NT s ).
Then, before entering the DPLL to recover the phase passes through a magnitude normalizer to guarantee the loop performance, regardless of the received signal level.The normalized signal at the DPLL input is given by where φ[k] represents the phase noise resulting from the Gauss- DPLL finally estimates the carrier phase, where  � [k] indicates the estimated phase.The operation of this DPLL is like the traditional Costas loop.However, as the input signal is already in baseband, it is not necessary for the I&D blocks to exist within the loop.In analogy with the Costas loop, possible discriminators for this DPLL can be implemented by the following operations:

Determination of Phase Error Variance
The portion of noise present in signal r F [k] (13) is given by where z D [l] is a discrete, complex, Gaussian, and white stochastic process with zero mean and an autocorrelation function where δ[m] is the discrete unit sample signal.

The autocorrelation function of z
where the operator E{x} indicates the statistical mean of x.From the definition of p[n] in (12) and from ( 15), (16), and (17), we get [ ] ( ) Therefore, the process z F [k] is cyclostationary with period M. To determine the power spectral density of z F [k], a new process derived from it and with a randomized delay [19] is defined by the equation where q is uniformly distributed over the interval [0,M).The autocorrelation function of z FS [k] is [19] [ ] Considering that the power spectral density of z F [k] is equal to the discrete-time Fourier transform of (20), we get ( ) The signal r F [k] can be represented by where z F [k] = z FX [k] + j z FY [k] is a complex Gaussian stochastic process, with power spectral density given by ( 21).The processes z FX [k] and z FY [k] are also Gaussian and independent with spectral density equal to S zF (e jω )/2.
Considering that e jπ = −1 and e j2π = 1, we can represent the symbol received in (22) as d[k] = e jψ [k] and ψ[k] is equal to π or 2π.In this way, we rewrite (22) in the following way: where z F [k] was replaced by z F [k]e j(θ[k]+ψ[k]) , because as both processes are zero-mean Gaussian processes with the same autocorrelation function, they are statistically identical.From (23), the magnitude and phase of r F [k] are determined, given by respectively, where the phase noise part is The signal r f [k] passes through the magnitude normalizer, described by ( 14), deriving the signal In this equation, θ [k] represents the phase of the carrier to be estimated and φ [k] the phase noise caused by the Gaussian thermal noise of the receiver.In high signal-tonoise ratio condition, that is, C/N 0 >> 1, the phase noise is approximated to Therefore, from ( 21), ( 22) and ( 28), and recognizing that , we derive the power spectral density of φ [k]: ( ) Assuming a third-order loop and the use of a DD discriminator, the DPLL linear model shown in Fig. 2 is illustrated in Fig. 3.

Using (3) and knowing that the transfer function of this loop is defined by H
Equation ( 30) shows that DPLL in the tracking phase works as a low-pass filter.Therefore, the variance of the estimated phase error can be determined by the expression [20] ( ) where B LD is the noise bandwidth of H B (z) and T B is the loop sampling period.In this case, we have T B = NT s /M.Setting ω = 0 in (29) and substituting in (31) results It is interesting to note that this expression differs from (9) only by the absence of quadratic loss, which in this case does not really exist because a DD discriminator was used.The result of this analysis is valid when the decision of the symbol used in the DD discriminator is correct.This is true under the condition that C/N 0 >> 1.

Numerical Results and Discussions
In this section, performance evaluation results of the traditional Costas loop and the proposed architecture are presented.Both architectures are evaluated via computational simulation, under various aspects of operation.A series of numerical tests were performed aiming to determine: 1.The variance of the estimated phase error (σ φ 2 ) as a function of C/N 0 .

Phase error (φ
of ionospheric scintillation and dynamic stress. 3. The root mean squared error (RMSE) of the estimated phase (φ RMSE ) in scenarios with ionospheric scintillation and dynamic stress.
4. The estimation of the probability of occurrence of cycle-slip or loss of lock [4] as a function of the scintillation intensity.
Specifically, the tests were performed with third order loops, having B LA = 10 Hz, T a = 20 ms and DD discriminator.The DPLL of the new architecture operates at the sampling frequency 1/T B = 2 kHz.
All models of the simulated systems were implemented in Mathworks Simulink running on a computer (3.7 GHz AMD Ryzern 5 5600X 6-Core processor) with a Windows operating system.
The block diagram of the entire system implemented in Simulink and used for simulations is described in Fig. 4. The signals x[n] and r[n] at the inputs of the loops to be evaluated are given by ( 1) and (10), respectively.The phase of these signals, denoted by θ [n], is modeled by where  The input and output parameters used in the simulations and indicated in Fig. 4 are described in Tab. 2.
The model of traditional Costas loop used in this article for comparison with the proposed architecture is exactly the same as the solution presented in [4], [14].The validation of this model was done via simulations in which results like those presented in figures 7 and 8 of [4] were obtained.Therefore, the performance results of the modified Costas loop will always be referenced to the performance of this discrete Costas loop found in the literature.
The plots of the variance of the estimated phase error (σφ 2 ) as a function of C/N 0 for the two architectures are shown in Fig. 5.It is evident from the results of Fig. 5 that the proposed architecture clearly outperforms the traditional architecture by approximately 3.5 dB.We can easily verify the blue curve points reasonably agree with (32).Ionospheric scintillation causes amplitude and phase fluctuations in trans ionospheric GPS signals [21].The statistical model CSM (Cornell Scintillation Model) synthesizes the effects of equatorial ionospheric scintillation for GPS links through a stochastic process c[i], as indicated in Fig. 4, that multiplies the GPS received baseband signal [22].The scintillation severity in CSM is defined by two parameters: the S4 index, which can vary from 0 to 1 and the decorrelation time (τ 0 ) [22].The phase RMSE curves for the tracking architectures, when subjected to scintillation with intensity {S 4 = 0.5, τ 0 = 0.48 s} and {S 4 = 0.2, τ 0 = 0.48 s} and C/N 0 = 35 dB, are plotted in Fig. 7.The simulation was conducted as follows: From 0 to 50 seconds, the architectures were submitted only to Gaussian noise and Doppler dynamics to have the RMSE convergence.From 50 to 150 s, the systems were also submitted to ionospheric scintillation.It can be observed that the proposed architecture performs better than the traditional Costas loop.
In sequence, Figures 8(a To evaluate the robustness to cycle-slips of the new architecture, in Fig. 10(a), 10(b) and 10(c) we plot the corrected estimated phase of both systems for a carrier-tonoise ratio of 35 dB-Hz; a Doppler shift of 1 Hz; and moderate (S4 = 0.5, τ 0 = 0.48 s), strong (S 4 = 0.7, τ 0 = 0.35 s) and severe (S 4 = 0.9, τ 0 = 0.2 s) ionospheric scintillation conditions, using a simulation time of 500 seconds.From the results shown in Fig. 10(a), in moderate scintillation scenario, we see that both architectures keep lock for 500 s interval, but it occurs a cycle-slip event for traditional Costas loop (red curve).From Fig. 10(b), in strong scintillation condition, we observe cycle-slips for both architectures but only traditional Costas loop (red curve) loss lock permanently.We notice, when the scintillation is severe, as shown by results in Fig. 10(c), although the cycle-slip events be very frequent for the modified Costas loop (blue curve), it remains locked at every interval while the traditional Costas loop fast losses the lock.As in the proposed architecture the equivalent noise band of the loop is constant, its performance in terms of phase RMSE and robustness to cycle-slip when subjected to scintillation is much lower than solutions based on adaptive Kalman filters [8][9][10].Scintillation changes the C/N 0 ratio and the loop would need to be continuously optimized for each C/N 0 value.On the other hand, the loop needs to have a minimum bandwidth to be able to track the phase of the scintillation.This requirement conflicts with the need to reduce bandwidth when C/N 0 is low.Therefore, only the most sophisticated solutions, such as those described in [8], [9,10], are suitable for strong scintillation scenarios.
In summary, in this section it was possible to validate the proposed carrier phase recovery architecture through several computer simulations, considering a wide range of scenarios operations, including Doppler dynamics, additive white Gaussian noise and even ionospheric scintillation.

Conclusion
In this work, we reviewed the traditional Costas loop that is used for carrier phase tracking in GPS receivers.We present the loop design as a whole and detail the possible types of discriminators and loop filters.We argue that the use of I&D as an integral part of the DPLL and the dependence of the sampling frequency on the integration time of the I&D cause degradation in the loop performance in terms of the phase error variance.
As a solution to these traditional Costas loop problems, we propose a new architecture for carrier phase tracking in a GPS receiver.
The proposed solution works with the received signal in baseband, removes the I&D from within the DPLL and increases the sampling frequency after the I&D, which in this case works as a matched filter.Increasing the sampling frequency reduces the noise bandwidth of the DPLL, reducing the estimated phase error variance.
An unprecedented mathematical analysis is presented to evaluate the performance of the proposed architecture.It is noteworthy that the value of the variance of the estimated phase error, determined analytically, was confirmed via computational simulations.
The Costas loop and the proposed architecture were simulated in different operating scenarios with the aim of comparing the performance between them.The proposed solution obtained lower variance of the estimated phase error, lower phase RMSE with Doppler stress and lower RMSE with scintillation.In term of error variance of estimated phase, the gain of proposed solution is about 3.5 dB for BLA = 10 Hz, T a = 20 ms.Regarding the robustness to cycle-slips the probability of cycle-slip or lock loss, we also confirm the superiority of the proposed architecture.
From the results shown, it is clear that the proposed modifications in the Costas loop bring expressive improvements in the performance of the carrier phase tracker in GPS receivers.
As a continuation of this research, we are studying new techniques to improve the robustness of the proposed architecture in operating scenarios with low C/N0 or ionospheric scintillation.A first line would be to replace the loop filter, which is fixed, for adaptive filters with variable bandwidth.Another possibility to be considered is to use DPLL based on the Kalman filter instead of traditional DPLL.
by the signal  � [n].After conversion to baseband, through the multipliers in each loop arm, the converted signals pass through the I&D blocks, which are integrated-and-dump accumulators with an integration period equal to T a .The frequency 1/T a is called the loop pre-detection bandwidth and corresponds to the sampling frequency of the resulting baseband signals I[k] and Q[k].These signals are used by the discriminator to generate the error signal e[k] which is proportional to the difference between the phase of the received carrier and the estimated phase.The error signal, smoothed by the loop filter, is delivered to the NCO (Numerically Controlled Oscillator) whose output is the estimated phase  � [n].Finally, the carrier generator block derives the signals cos(ω c n +  � [n]) and -sin(ω c n +  � [n]), thus closing the feedback loop.The discrete frequency ω c corresponds to the intermediate frequency (IF) of the received carrier.

1 .
AT: e[k] = arctan(Im{e B [k]}/ Re{e B [k]}); 2. CC: e[k] = Re{e B [k]} Im{e B [k]}; 3. DD: e[k] = sign(Re{e B [k]}) Im{e B [k]}.The mathematical operations Re{x} and Im{x} indicate the real and imaginary parts of x, respectively.The loop filter and NCO blocks, in this case, are identical to those used in the traditional Costas loop.The carrier generator delivers the complex signal exp(-j � [k]) to loop feedback.
denoted by Doppler phase resulting from the system dynamics, with θ 0 being the initial phase and f d being the residual Doppler shift.The phase contribution due to scintillation is indicated by θ s [n].The scintillation generator delivers the discrete process c[i],determined according to the CSM model[22], to the signal generator.The phases estimated by the two loops are  � [n] and  � [k] as shown in Fig.4.The measure processing block from the estimated phases and the generated phases determines the phase error (φ[n]), the variance of the phase error (σφ 2 ), the root mean squared error (RMSE) of the estimated phase (φ RMSE ) and an indication of cycle-slip (CS).

Figure 6
Figure 6 shows phase RMSE plots as a function of time for the two loops in the condition of C/N 0 = 35 dB and residual Doppler shift of 1 Hz.From the curves, the superiority of the modified Costas loop is clearly demonstrated.

Fig. 8 .
Fig. 8.Comparison of corrected estimated phase for both loops for weak scintillation.

Fig. 9 .
Fig. 9. Comparison of corrected estimated phase for both loops for moderate scintillation scenario.
), 8(b), 9(a) and 9(b) show the corrected estimated phase, that is estimated phase subtracted from Doppler phase ( � [n] − θ d [n]), of both traditional and modified Costas loop, during events of weak (S 4 = 0.2, τ 0 = 0.48 s) and moderate (S 4 = 0.5, τ 0 = 0.48 s) ionospheric scintillation.It can be seen from the results that both architectures are capable of tracking the scintillation phase.It is also evident that the modified Costas loop tracks phase scintillation better than the traditional Costas loop.

Figure 11 Fig. 10 .
Figure 11  shows the curves that represent the estimation of probability of cycle-slip or loss of lock in function of ionospheric scintillation intensity (S4) for the two studied

Fig. 11 .
Fig. 11.Probability of lock loss or cycle-slip as a function of S 4 .architectures.Simulations of 60 seconds were performed 50 times.Only Gaussian noise, with C/N 0 = 35 dB, is inserted in received signal during the first 10 seconds of each simulation.It is necessary to have the acquisition of lock by the loops initially.The scintillation is inserted at signal only after 10 s.It is evident by the results presented in Fig. 11 that the proposed solution is more robust to cycleslip or loss of lock than the traditional Costas loop.