Optimized Carrier Tracking Loop Design for Real-Time High-Dynamics GNSS Receivers

Carrier phase estimation in real-time Global Navigation Satellite System (GNSS) receivers is usually performed by tracking loops due to their very low computational complexity. We show that a careful design of these loops allows them to operate properly in high-dynamics environments, that is, accelerations up to 40 g or more. Their phase and frequency discriminators and loop filter are derived considering the digital nature of the loop inputs. Based on these ideas, we propose a new loop structure named Unambiguous Frequency-Aided Phase-Locked Loop (UFA-PLL). In terms of tracking capacity and noise resistance UFA-PLL has the same advantages of frequently used coupled-loop schemes, but it is simpler to design and to implement. Moreover, it can keep phase lock in situations where other loops cannot. The loop design is completed selecting the correlation time and loop bandwidth that minimize the pull-out probability, without relying on typical rules of thumb. Optimal and efficient ways to smooth the phase estimates are also presented. Hence, high-quality phase measurements—usually exploited in offline and quasistatic applications— become practical for real-time and high-dynamics receivers. Experiments with fixed-point implementations of the proposed loops and actual radio signals are also shown.


Introduction
A fundamental task of every Global Navigation Satellite System receiver is to synchronize with the visible satellite signals. Since Direct Sequence Spread Spectrum (DS-SS) signals are utilized, code and carrier synchronization is required, but a correlation stage is necessary to despread the signals before the synchronization algorithms can be applied. In real-time receivers the required economy of operations usually precludes the use of complex estimation schemes and tracking loops are preferred. Due to the correlation process these loops are necessarily discrete. The typical trade-off in tracking loop design is bandwidth versus dynamic performance: output noise increases with a larger loop bandwidth, while dynamic tracking error decreases with it [1]. Thus, the loop design becomes particularly challenging when the receivers are subject to high dynamics. To overcome this limitation other receiver structures have been proposed in [1], claiming tracking capability up to 150 g of acceleration, in contrast with the 5 g regularly assigned to tracking loops. However, the required computational burden is large since several simultaneous correlations and Fast Fourier Transform (FFT) computations are needed. In this paper we show a careful design of the digital loops that can expand their tracking ability to acceleration steps up to 40 g or even more, keeping a low computational load and reasonable tracking threshold values at the same time.
The loop structure known as FLL-assisted PLL [2] is very often adopted for GNSS receivers. It consists of a Phase-Locked Loop (PLL) and a Frequency-Locked Loop (FLL) in a coupled mode, with the advantage of reducing locking times and avoiding false locks. This solution is also a legacy of analog loops since the FLL or Automatic Frequency Control (AFC) has been used to reduce frequency errors as a previous stage to phase lock for analog PLL [3]. The advantages of adding the FLL to track spread spectrum signals in dynamic environments were already studied in [4]. For high-dynamics GNSS receivers, the focus is on carrier loops because the carrier shares the same dynamics as the code. Then, the estimation of the carrier frequency can be used to aid the estimation of the code frequency, and a first-order code loop is enough [5]. Usually, implementations of FLL-assisted PLL are not based on optimal digital loop solutions, with each loop designed separately, leaving the analysis of their interactions and possible modifications to the simulation stage [2,5,6]. Moreover, schemes adopted to discriminate phase or frequency errors are often justified because of their similarity with well-known analog solutions rather than with an optimality versus implementation complexity criterion. We will show that digital implementations of optimal discriminators are not necessarily more complex and allow designing the FLL-assisted PLL in a coupled way.
Nevertheless, the FLL-assisted PLL leads to a more complex design and a computationally more expensive implementation than a single PLL. Moreover, when coupledloops lose phase lock for a moment, they present cycle slips introducing a phase ambiguity. We will show how to use the same frequency information as that of an FLL to build a nonambiguous phase detector, the Unambiguous Frequency-Aided (UFA) phase discriminator. A PLL with this new phase discriminator, that is, a UFA-PLL, keeps the desirable properties of an FLL without demanding an extra loop and avoiding cycle slips. Other nonambiguous phase discriminators are known for analog PLLs, that is, with analog loop filter, such as the sequential discriminators built with flip-flops presented in [7,8] or the nonsequential discriminator of [9]. While their goals are quite similar to ours, they increase the PLL implementation complexity, demanding some digital circuitry and a digital-to-analog converter to get the analog phase error. On the contrary, the UFA phase discriminator is easily implemented and naturally suited for a software-based PLL, leading to a less complex implementation than a FLL-assisted PLL. Section 2 introduces the UFA-PLL structure for GNSS tracking loops.
The optimum loop filter structure for analog PLLs was introduced in [10], solving the mentioned bandwidth tradeoff by minimizing a quadratic functional. A widespread technique for designing digital loops is discretizing an analog loop with a sample rate 1/T at least ten times faster than loop bandwidth B N [5,11,12]. As B N T increases above the rule-of-thumb value of 0.1, the resulting loop deviates from optimal and may become unstable [5], especially when accounting for the delays of a digital implementation. This limit imposed to the loop bandwidth is not fundamental and an attempt to avoid it has been presented in [13]. They introduced a digital loop design based on pole placement that allows somewhat larger B N T values. However, the pole location is assigned with standard second-order analogsystem rules. Our approach is to consider a completely digital loop model and pose the bandwidth trade-off directly in the digital domain, building upon the early and often overlooked work of [14] for hybrid loops. We include two delays in the loop to consider the effect of the correlation stage, similar to the inclusion of an accumulator before the loop error discriminator for signals without spreading codes [15]. Our method [16] allows the design of stable loops with B N T > 0.1, a particularly useful feature for high-dynamics receivers. Specifically, we will focus on dynamics modeled as acceleration steps, that is, unbounded jerk, as in the case of launching vehicles when the engine turns on or off. In Section 3 we first derive the optimal loop filter for arbitrary phase inputs and then for the case of acceleration steps that produce quadratic ramps of input phase or a linear ramp of the input frequency. Simulations comparing the different loop structures are also shown.
Optimization gives the structure for the loop filter, leaving the choice of T and B N unsolved. Usually, these parameters are selected based on some rule of thumb [2,5], and the ultimate loop performance, as measured by the pull-out probability and/or tracking threshold, is obtained later by simulation. An optimal choice of these fundamental parameters demands an analysis of the nonlinear aspects of the tracking loop with noise. This is quite difficult, although some results are known for analog loops with relatively simple loop filters, by solving a Fokker-Planck equation [17]. They can be extended to digital loops when an analog approximation is valid [18]. Our approach is to get a reasonable approximation for the pull-out probability and its relationship to the loop parameters. This new approach introduced in [19] allows us considering dynamics modeled as acceleration steps and digital loop filters with zero stationary error response to these inputs. Previous analyses are based on stationary loop responses or sinusoidal acceleration profiles [2,5]. For these cases, we derive approximate expressions for the probability of starting a nonlinear behavior of the mentioned loops. These expressions quantify the role of B N and T and let us choose them in order to obtain lower tracking thresholds for different dynamic scenarios, as presented in Section 4.
Our optimized digital carrier tracking loops also allow smoothing of the phase estimates incorporating more measurements, at the expense of some delay. In general, an output delay of a few samples should not be a limitation since the navigation task in a GNSS receiver is usually slow compared with the loop sample rate. This update of the phase estimates can significantly reduce the noise variance and the transient responses in high-dynamics environments. This strategy is suitable for real-time receivers because it can be efficiently calculated. Hence, some of the precise positioning techniques would be applicable in real-time and for highdynamics receivers. Consider, for instance, smoothing of code delay measurements with carrier phase estimates in stand-alone receivers [12,20], or differential positioning applications [12,21], or even attitude estimation with GNSS signals [22]. In all these cases, an improvement in the phase estimation has a direct impact on the positioning performance. The expressions for optimal smoothing filters are derived in Section 5, and their efficient implementation is also discussed there. In addition, we present experimental results obtained with actual RF signals and a fixed point implementation of our loops tracking acceleration steps of up to 40 g. Finally, the conclusions of this work are given in Section 6.

Digital Loops Models
Correlations of the received signal with the locally generated replicas for each visible satellite are the inputs to the discriminator of the carrier tracking loops in a GNSS receiver. The complex correlation for a given satellite with International Journal of Navigation and Observation carrier power to noise power spectral density C/N 0 and for the ith correlation interval of duration T can be written as [12] where Δτ i = τ i − τ i is the code delay estimation error, Δ f i = f i − f i the frequency estimation error, and Δθ i = θ i − θ i the phase estimation error, all assumed constant during the integration time. The sequence n i is a complex white Gaussian noise process with unit variance, R(·) is the code correlation function, and sin c(x) = sin(πx)/(πx). It is also assumed that the signal has binary data bits D i = ±1 and that correlations are computed within a single bit period. This Binary Phase Shift Keying (BPSK) modulation is present in many GNSS signals like the GPS civil signals or in the data components of composite modernized GNSS signals [23]. In tracking conditions (i.e., after the acquisition process has been completed [12]), estimation errors are small and then the functions sin c(·) and R(·) can be approximated by 1. In this case the expression (1) reduces to where I i and Q i are the so-called in-phase and in-quadrature correlations respectively, n Ii = e{n i }, n Qi = m{n i }, and we have defined These sequences allow to model the carrier tracking loop as a purely digital single-input single-output (SISO) system. When the frequency is changing according to a constant acceleration error of a m/s 2 , we verifiedby numerical integration-that expression (2) is a good approximation if aT 2 /λ 1, where λ is the wavelength of the signal. For L1 GPS λ = 0.19 m and with T = 5 ms, this implies that a 7600 ≡ 775 g. In this case, the terms Δ f i and Δφ i have to be reinterpreted as the average frequency error and average phase error during the correlation interval, respectively.
In the following we briefly review the basic concepts of PLL and FLL-assisted PLL from our digital point of view, and later we introduce the UFA-PLL.

PLL Model.
The phase estimation error is typically obtained using one of several possible discriminators [5], which give the desired phase modified by different memoryless nonlinearities. The optimal one-maximum likelihood estimator-is given by where the notation [·] π indicates that its argument is kept within the interval (−π/2, π/2] by adding or subtracting π as many times as needed. The zero-mean noise term n φi has a rather complicated probability distribution [24], but in high C/N 0 it can be approximated by a Gaussian distribution with zero mean and variance σ 2 φi ≈ 1/(2TC/N 0 ). A four-quadrant tan −1 (·) is not appropriate if there is BPSK data modulation because the discriminator becomes sensitive to the data phase changes. On the contrary, for signals without data the range of the discriminator can be doubled with a four-quadrant tan −1 (·). We chose this discriminator because it is not amplitude dependent and the calculation of tan −1 (·) can be easily implemented with a lookup table, since in practice I i and Q i are frequently quantized to a few bits.
In order to close the loop in our model, it is of crucial importance to consider the delays present in a real implementation. Failure to account for a delay may turn unstable an optimal loop design. Since ours loops are digital, a single sample delay is expected but in fact there are two. One of them is due to the time spent in I i and Q i calculations.
The other delay appears because the estimated values used in the present correlations have to be known before the calculations begin. That is, the value φ i is obtained with the loop filter output of the (i − 1)th correlation interval, which in turn is calculated with the estimation errors of φ i−2 . Then, with these considerations, the model of a PLL using the classical loop filter structure of type 3, that is, with three accumulators, is shown in Figure 1.

FLL-Assisted PLL Model.
To add an FLL to our previous PLL, a frequency discriminator is needed. In a digital loop a frequency error estimate may be obtained as the difference of two successive phase errors, and in fact this is often correct. A problem appears when the discontinuities caused by [·] π make that the difference to be wrong in ±π. However, our discrete system cannot distinguish frequencies greater than half of the sample rate, that is, phase changes of π between consecutive samples, and so the measured frequency errors must be bounded. In fact, if the phase discriminator is insensitive to BPSK data, the phase changes caused by 4 International Journal of Navigation and Observation Figure 2: Block diagram of the typical FLL-assisted PLL structure.
frequency errors must lie in the interval (−π/2, π/2] [25]. Thus, the difference of two consecutive outputs of the phase discriminator can be corrected just using the operation [·] π . Therefore, the frequency discriminator for the FLL can be obtained by (4) Figure 2 shows a diagram of the FLL-assisted PLL presented in [2], where the second-order loop filter of the FLL shares the same cascade of accumulators used by the PLL filter. In the locked condition e i = Δφ i and e fi = Δφ i − Δφ i−1 are small enough to justify a linear analysis of the loop. The complete loop is seen as an equivalent PLL with filter coefficients p 3 , p 2 + f 2 , and p 1 + f 1 , instead of p 3 , p 2 , and p 1 . Thus, the FLL is inserted into the model of the PLL at a design stage. This eliminates the constrain of using a narrow bandwidth FLL to not significantly perturb the PLL behavior, as in [2,6]. A wide bandwidth FLL allows the loop to have two regions of operation: "phase-locked" as it was described before, and "frequency-locked" when the dynamics unlocks the PLL but the FLL keeps the frequency error within the linear range of its discriminator. In the latter region the loop is governed by the FLL (coefficients f 1 and f 2 ) and the phase error input acts like a zero-mean perturbation [25]. As soon as the dynamics let the loop reduce its frequency error close to zero, the phase lock can be restored.

The UFA-PLL Model.
As we have seen so far, due to the cyclic nature of phase a memory-less discriminator is unable to distinguish changes of an integer number of cycles-or half cycles if there is BPSK data-, that is, its output is ambiguous. However, it is possible to obtain a frequency error estimate from these ambiguous phase error estimates correcting their difference with the nonlinear operation [·] π . This is the reason why an FLL can cope with carrier tracking in situations when a single PLL cannot. Assume that there is BPSK data and the PLL phase error is rising and crosses the value π/2. The output of the phase discriminator abruptly changes to a value close to −π/2, reversing the evolution of the PLL phase. Hence, the phase error will increase since the PLL is now moving in the wrong direction. We should instruct the phase discriminator with information of the phase derivative to keep moving in the right direction, that is, we should feed it with proper frequency information available at the FLL. Therefore, the idea of the Unambiguous Frequency-Aided (UFA) phase discriminator is to use the same frequency information used by the FLL to get a better phase discriminator. It works correcting the ambiguous values of e i by adding or subtracting an integer number of π so that the difference of successive values of the corrected phase error, u i , gives the right frequency error. Then, the equations that define our new phase error estimate, for i ∈ N, are where we define I π (x) = x − [x] π , an operation similar to the function integer part but with steps at the multiples of π. A practical formula to compute k i can be derived noting that and then k i π = −I π (e i − u i−1 ). Substituting this in (5), we can recursively calculate the UFA phase error from the ambiguous e i : with starting value u 0 = e 0 . Then, the PLL structure in Figure 1 transforms into a UFA-PLL just adding a block that implements (7) immediately after the phase discriminator output, as shown in Figure 3.
It is interesting to note that the UFA scheme acts like the phase unwrapping algorithm proposed in [18] for correcting cycle slips in the phase estimates of feed-forward synchronizers. In this case, the phase correction does not affect the phase estimation process since it is done once the estimation stage is finished. On the contrary, the UFA phase discriminator modifies the behavior of our feedback estimator, the PLL, changing its nonlinear characteristics. As a result, cycle slips and the rather complex transient responses induced by them are avoided as long as the frequency error is compatible with the loop sample rate.

Equivalence between UFA-PLL and FLL.
We saw that the frequency error estimate can be obtained as the difference of two successive phase errors if the result is kept in range by adding or subtracting an integer number of π. Thus, the frequency discriminator for the FLL can be obtained International Journal of Navigation and Observation Figure 3: Block diagram of the new UFA-PLL structure.
correcting the difference of two consecutive outputs of the phase discriminator just using the operation [·] π . Figure 4(a) shows a block diagram of a digital FLL, with loop filter transfer function F(z). Notice that the two delays and the accumulator that convert the frequency estimation to phase before the feedback are not included in F(z).
An alternative way to obtain the same frequency error discriminator is to use the UFA algorithm previously described. Indeed, the output phase sequence u i is built in such a way that the difference of consecutive values produces the right frequency error, as seen from (4) and (5). Therefore, the schemes of Figures 4(a) and 4(b) are equivalent. The interesting fact in Figure 4(b) is that most linear blocks are adjacent. Thus, the differentiator cancels with the accumulator without changes in the dynamic loop response, except for the mean value of the phase error, leading to the equivalent UFA-PLL model of Figure 4(c). In fact, this zero-pole cancellation shows why the FLL is insensitive to constant phase errors whereas the equivalent UFA-PLL is not. More importantly, the equivalence reveals that the nonlinear behavior of the UFA-PLL is equal to that of a FLL with the same F(z), and then their tracking capacity and noise resistance are the same.

Optimal Loop Filter Design
We propose to design the digital loop filter minimizing a twoterm quadratic functional to handle the bandwidth trade-off mentioned in the Introduction. The input signal is assumed to have a part related to phase evolution φ i plus additive, zero mean,and noise n φi . The functional to be optimized is where α 2 is a weighting factor that controls the tradeoff between noise and transient response, that is, the loop bandwidth, σ 2 N is the noise variance at the loop output, and E T (φ i ) represents the energy of the tracking error Δφ i transient response. Since the functional uses the energy of the transient response, the optimum filter must produce a zero stationary response for the given input.
Suppose F(z) is the loop-filter transfer function to be found, and consider that the linear model hypothesis holds for a PLL or FLL. The closed loop transfer function including the delays is where Y (z) must be a causal and stable rational transfer function. The minimization of the functional J written in terms of Y (z) is shown in the Appendix A. The optimum transfer function is given by where Ψ(z) and X(z) can be obtained from the spectral factorization of (A.6) and from the partial fraction expansion of (A.8), respectively. We repeat them here for completeness, where Φ(z) is the z-transform of φ i . The relation between minimizing (8) to Wiener filtering [26] is sketched in Appendix B. Observe that whereas the former is a mixed criterion with a term depending on the stochastic part and a term depending on the deterministic part of the phase signal; the latter criterion stems from a purely stochastic formulation. The connection between both approaches arises when modeling the input phase as white noise passing through a rational transfer function. The Wiener filtering approach offers other possibilities such as keeping the optimality for a wide range of admissible transfer functions via a robust approach as in [27] or considering continuous models for the phase as in [28]. Optimum loop filters for an input phase step, frequency step, and frequency ramp were derived in [16]. In the following, only the last result is presented for the sake of brevity. Analog loop filters optimized for these kind of inputs are the origin of the classical methods of filter design for type one, two, and three loops, respectively. As it will be seen, our purely discrete design for each case has one extra pole, due to the loop delays. This additional pole does not appear when discretizing analog designs, but it has a decisive influence on the stability or the range of achievable product B N T.

Optimum Filter for a Frequency Ramp. The ramp is modeled as
whereΔω is the rate of frequency change. Denoting ν = Δω 2 T 4 γ 2 , from (11) it is necessary to solve (z − 1) 6 − νz 3 = 0. The six roots of this polynomial are obtained using the fact that three of them are the inverses of the other three. This allows us to express the following equations: that determine the values of z 1 , z 2 , complex conjugates and a real z 3 . Using these values and (11), we get and replacing in (12) Then, the corresponding X(z) has only three poles in z = 1, and the closed-loop transfer function of (10) is where Then, the optimum loop filter with four poles, three of the input and the extra one, is . For the purpose of implementation, it is desirable to use a cascade of accumulators. Then, (18) can be rewritten as where The closed loop noise equivalent bandwidth is shown in Figure 5. This curve allows choosing the appropriate value of ν for a given normalized noise bandwidth. Observe also from this figure that the larger is the step or the parameter ν, the more emphasis is given to the transient energy, causing the normalized noise bandwidth to increase. The product B N T levels out to a value of approximately 54.5. This part is not included in the figure since it is of minor importance for most designs of practical interest.

Design Example: Loops for Launching Vehicles.
We simulate and compare the loop models presented in Section 2 taking as an example a GPS carrier tracking loop for launching vehicles [29]. In this case the dynamic input can be modeled as an acceleration step, which becomes a quadratic ramp in terms of phase and a linear one, in terms of frequency. For these inputs the optimal loop filter for a PLL was obtained in Section 3.1. The case of an FLL-assisted PLL, taking results from [29], leads to a type 2 FLL and a type 3 PLL. The FLL-assisted PLL in [29] was designed to operate in "phase-locked" mode with steps up to 10 g and in "frequency-locked" mode up to 20 g of acceleration. These requirements were too demanding for the commonly adopted correlation time of 10 ms, and then it was lowered to 5 ms (at the cost of almost doubling the processor load and an increase in the tracking threshold). A typical rule of thumb for keeping a reasonable distance from the pull-out values of the loop is that the peak of the error transient has a maximum value given by half the linear range of the phase discriminator, an eighth of cycle [2]. As it will be shown in the simulations this condition was obtained with a value of ν = 0.00025 for the PLL. However, as we will explain Figure 6: Block diagram of the FLL-assisted PLL in Section 3. Figure 7: Block diagram of the UFA-PLL designed in Section 3. in Section 4, this rule is not that useful. The resulting filter transfer function for the PLL was where Since the FLL design does not affect the previous results, it was designed wider than strictly necessary in order to facilitate the posterior implementation. The selected transfer function is where D = 0.6 and E = 0.5, resulting in that the extra pole needed for the FLL and the PLL is the same. Then, f 1 = E = 0.5 and f 2 = D − E = 0.1. This implis p 1 = 0 and p 2 ≈ 0. With these simplifications the complete loop design reduces to the diagram showed in Figure 6. This FLL loop can track steps up to 40 g with transient error peaks smaller than 25 Hz, half of the linear range of the frequency discriminator, with an equivalent noise bandwidth of B N = 61.3 Hz. We will use the previous loop filter as a basis for the comparison of different loop configurations. We consider a PLL and a UFA-PLL with the same loop filter as before, that is p 1 = C = 0.5, p 2 = 0.105, p 3 = 0.0123. This structure is showed in Figure 7. The phase error response for a step of 10 g of acceleration, common to the three loops as expected, is depicted in Figure 8. In this case the phase error detected by the discriminators is equal to the actual phase error since its magnitude is always less than a quarter of cycle.
In Figure 9 the discriminated phase error during a 40 g step in the three loops is illustrated. It can be seen that In the case of the FLL-assisted PLL the loop lost phase lock for a moment, but could still track the dynamics because the FLL remained locked. This nonlinear behavior could correspond to a cycle slip. To verify this analysis the actual phase error for each configuration is shown in Figure 10.
Clearly, there is one cycle slip in the tracked phase of the FLLassisted PLL, whereas the UFA-PLL is able to track this step of acceleration without any cycle slip. The limit of the tracking capability of the FLL-assisted PLL and the UFA-PLL is the frequency error-the phase change between samples-, since both are frequency-aided loops. If this error becomes greater than 50 Hz in magnitude, the frequency estimation will be ambiguous, in the same way as the phase of the PLL does, and the loops will lose their lock. This can be caused by excessively high dynamics or noise power or a combination of both. Notice that the noise power considered in this case is twice the input phase noise power due to the differencing. In fact, we can use the UFA algorithm applied to the frequency discrimination to further extend the dynamics resistance of the loops. However, it will be of little practical importance due to the noise power increase caused by a new differentiation. The frequency error of both loops during a 40 g step is shown in Figure 11. The peak error is 25 Hz-half the limit-and thus, using the same rule of thumb that we used for the phase error of the PLL, it can be argued that 40 g is the level of acceleration steps that can be tracked with a reasonable safety margin for noise effects. In Section 4 we will give a totally different approach for the consideration of the noise in the UFA-PLL. Another shared feature of the FLL-assisted PLL and the UFA-PLL is the resistance to false locks. For the sake of brevity this analysis is not included here but it can be seen in [16].

Pull-Out Probability Analysis
In this section we will compute an approximation to the pullout probability for a loop in a given operating condition.
Since the pull-out is necessarily a consequence of nonlinear behavior, a simple way to bound the pull-out probability is with the probability of entering nonlinear behavior, that is, the limits of the tan −1 (·) range. We are interested in tracking acceleration steps. Then, it is clear that for a given noise level the instants when the loop is closer to these boundaries approximately correspond to the peaks of the loop transient response. Therefore, we focus on calculating the probability of entering nonlinear behavior at the instant of the transient response peak, given that the loop behavior has been linear up to that time.

PLL Analysis.
The PLL enters nonlinear behavior when the phase error becomes larger than π/2. This is equivalent to a sign reversal of the in-phase component, with respect to the sign of the data bit. Then, the probability of nonlinear behavior at the transient peak i = p is [19] P P = P Δφ p + n φp > π 2 = P cos φ p − φ p + n Ip < 0 , where n Ip = D p n I p / TC/N 0 has variance σ 2 = 1/(2TC/N 0 ). Assuming that the PLL has had a linear behavior up to the analyzed instant, φ p , can be thought of as a deterministic value plus output noise n φp . Subtracting this deterministic value from φ p we find the peak value of the loop transient response to an acceleration step, which is shown in Figure 12 as a function of ν for an integration time of 10 ms. Note that with acceleration we refer to the Doppler rate the loop has to track, scaled in units of g = 9.8 m/s 2 instead of Hz/s to keep an easy physical interpretation. This peak value is proportional to the square of the integration time because an acceleration step is a parabolic ramp of phase. The values of the normalized closed-loop equivalent noise bandwidth are also plotted in Figure 12 for completeness. The noise n φp depends on the filtered past input noisefrom instants p − 2, p − 3, and so on-and is statistically independent of n Ip . Then, we can write φ p − φ p = K p aT 2 + n φp , where a is the acceleration amplitude in g's. In (22), we obtain P P = P cos K p aT 2 + n φp + n Ip < 0 .
To further simplify (23), we use a "low noise" approximation cos(n φp ) ≈ 1 and sin(n φp ) ≈ n φp and obtain P P ≈ P cos K p aT 2 − sin K p aT 2 n φp + n Ip < 0 .
Using central limit theorem arguments, it is reasonable to approximate the distribution of n φp as a zero mean Gaussian with variance approximated by 2B N Tσ 2 φi ≈ B N /(C/N 0 ), using the high C/N 0 variance expression for the input noise. In this case, both random terms in (24) become a single Gaussian random variable with zero mean and variance (1 + 2B N Tsin 2 (K s aT 2 ))σ 2 . Therefore, where Q(x) is the cumulative Gaussian distribution from x to ∞. Since the function Q( such that the larger is f PLL , the smaller is P P for a given C/N 0 , and then a low tracking threshold is attained. For example, f PLL (ν, T, a) for a 5 g acceleration step is plotted in Figure 13. It can clearly be seen that the larger values of f PLL are found in the region 5 ms < T < 7 ms and 10 −3 < ν < 10 0 , and the maximum is approximately at ν = 0.05 and T = 6 ms. The value of f PLL in this region is about 0.009. Larger values of ν are not preferred because they lead to loops that produce larger phase estimation error and only a slightly lower P P . However, it must be emphasized that for values of ν > 10 −3 the loop bandwidth is B N T > 0.5. These values of bandwidth show that a filter loop design based on discretization of analog solutions, only valid for B N T < 0.1, is not appropriate to design loops with better tolerance to nonlinear behavior. If ν = 0.02, then B N T = 1.5 and, with T = 5 ms, leads to an optimum loop bandwidth B N = 300 Hz, which is too large from the point of view of output phase error variance. This shows that the main cause of pull-out in PLLs with narrow bandwidths is the transient error response, due to the input phase and the input noise, rather than output noise. In other words, 5 g of acceleration is very demanding for a single PLL and then a large bandwidth is required to track them with small pull-out probability.

UFA-PLL
Analysis. An equivalent description of the UFA algorithm presented in Section 2 is to consider it as a modified tan −1 (·) function that produces output values in the range (−π/2 + u i−1 , π/2 + u i−1 ] instead of (−π/2, π/2]. Hence, we conclude that nonlinear behavior will occur if the actual phase error differs from the previously discriminated one by more than π/2. Assuming that the behavior before the analysis time has been linear, u i−1 = Δφ i−1 + n φi−1 and then Writing the phase noise terms as a function of the corresponding in-phase and in-quadrature components, it can be shown that this condition is equivalent to a sign reversal of International Journal of Navigation and Observation The quadratic terms in (28) have zero mean, and σ 4 variance and are uncorrelated between them and with the linear ones. For practical values of σ 2 their variance is much smaller than the variance of the linear terms. Even in this case, they cause the probability distribution of the sum in (28) to differ considerably from Gaussian. A more accurate calculation will require numerical computations of the actual distribution. On the contrary, our aim is to get reasonable and easy-to-handle approximation and, then, we will discard them, but being aware thatthe Gaussian assumption is only a coarse approximation. Then, P U ≈ P cos K u aT 2 + Δn φd + n Ii + n Ii−1 < 0 (30) where Δn φd = n φd − n φd−1 . The input noise terms are independent of each other and of the output ones because of both loop delays. Therefore, comparing (30) with (23) we replicate the reasoning for the PLL, but doubling the input noise contribution, and changing K p by K u and B N T by B N T. Then, using the same approximations made for (25), we get and then f UFA (ν, T, a) = Tcos 2 K u aT 2 1 + B N (ν)T sin 2 (K u aT 2 ) (32) is the function to analyze which values of ν and T are better for the design of low tracking threshold loops. For accelerations of 20 g f UFA (ν, T, a) is plotted in Figure 15 as an example. It can be seen that the larger values of f UFA are found in the region of T near 5 ms and 10 −4 < ν < 10 −2 , and the maximum is approximately at ν = 0.0025 and T = 5.5 ms. The value of f UFA in this region is about 0.004. In this case, the optimum loop bandwidth is about B N = 70 Hz, which is a more reasonable value than in the case of the PLL. For the UFA-PLL, the minimum pull-out probability and minimum output variance seem not to be as contradictory criteria as for the PLL. This can be understood noticing that the ability of the UFA-PLL for tracking in high dynamics depends on the smoothness of the transient error response rather than its absolute value, and then the output noise contribution becomes more relevant. Another important fact that must be emphasized is that even for 20 g accelerations it is not advisable to use correlation times T lower than 5 ms. Notice that, as it was explained in Section 2, this probability analysis applies also for an FLL as long as the right F(z) is used in the computations of K u and B N T.
International Journal of Navigation and Observation

Simulations.
In this section we assess the accuracy of the approximations made in the previous analysis. According to them the filter design of Section 3 is almost optimum when used in a UFA-PLL for tracking steps of 20 g, that is, it produces the minimum tracking threshold. The loop evolution with input noise according to a given C/N 0 was simulated until the transient peak instant. Runs entering the nonlinear region before this peak were discarded. A variable number of runs were used in order to reduce simulation time as much as possible but keeping statistical significance of the results. Specifically, 100,000 runs were enough for the lowest C/N 0 values, whereas 100 million had to be used for the highest. A step of 5 g is considered in the simulation of the PLL, which produces a transient peak of K p 5 g (5 ms) 2 = 1 rad. Using (26) we found f PLL = 0.0033. The results of the simulations compared with expression (25) are presented in Figure 16. It can be seen that the approximation is slightly optimistic and that the error is almost constant in the simulated range of C/N 0 . For the UFA-PLL we adopt a step of 20 g that produces a transient peak of K p 20 g (5 ms) 2 = 2 rad and K u 20 g (5 ms) 2 = 0.38 rad of peak difference. Using (25), we get f UFA = 0.004. The results of the simulations compared with expression (31) are shown in Figure 17. In this case the approximation is still acceptable for tracking threshold determination, but now it is pessimistic and the error grows for increasing values of C/N 0 . This behavior is caused by the Gaussian approximation that neglects the quadratic noise terms in the probability expression (28).

Tracking Threshold Analysis.
To illustrate how this analysis can lead to loop designs with lower tracking thresholds we consider the design of a type 3 UFA-PLL for 20 g acceleration steps. The tracking threshold is determined by a given probability of pull-out, and then we will define it by a given level of probability of starting a nonlinear behavior. Therefore, if P is the admissible pull-out probability and C/N 0 | TH is the tracking threshold, then or equivalently, Clearly, minimum tracking threshold values will be achieved when f UFA (ν, T) is near its maximum, that is, using ν ≈ 0.00025 and T ≈ 5.5 ms. Considering ν = 0.0003 and T = 5 ms, for example, the loop parameters are B N = 80 Hz, B N = 51 Hz K p aT 2 = 2, K u aT 2 = 0.38, and f UFA = 0.004. Taking for instance a value of P = 0.001 and replacing in (34), we find that C/N 0 | TH ≈ 34 dB/Hz. Another design, similar to those based on analog prototypes, can use ν = 0.000001 and T = 2 ms and then B N T = 0.1. The loop parameters become B N = 50 Hz, B N = 13 Hz K p aT 2 = 1.4, K u aT 2 = 0.1, and f UFA = 0.002. Hence, for the same P = 0.001 if we replace in (34) the tracking threshold results C/N 0 | TH ≈ 37 dB/Hz. Therefore, the use of the digital design method together with the proposed pull-out probability analysis can lower 3 dB the tracking threshold compared with traditional analog-based designs. An additional advantage is the use of longer values of T, requiring less computational load than analog designs.
It has to be mentioned that, even though the actual probability distribution can be different because of the Gaussian approximation, only the arguments of the Q(·) are used in the comparison of both designs. Therefore, the comparison is not affected by the Gaussian approximation, that is, a modified Q(·) function could be used for a more accurate probability calculation but the 3 dB threshold gain would remain.

Optimal Smoothing of the Phase Estimates
Due to the presence of two delays in the loop, the phase estimate obtained at a given instant is not computed with measurements up to this instant, but with measurements up to two previous instants. Thus, in the notation of [30] the loop phase estimate at the feedback branch is actually φ[i | i−2] = φ i . Naturally, the use of "closer" measurements would produce a smoothing effect, that is, a better estimation. The real-time constraint does not allow taking advantage of this for the loop itself, but it is possible for other purposes as data detection and raw data generation for the navigation processes of the GNSS receivers. In this case, the optimal phase estimate can be obtained in the same way as before, but without forcing the two delays as in (9).

One Sample Smoothing.
In this case the problem is equivalent to obtaining the optimal loop filter with only one delay. Then, if F (z) is the new loop-filter transfer function to be found, the corresponding closed loop transfer function is where Y (z) is the rational and stable transfer function to be found minimizing J(Y (z)) in (8). Following the same optimization process done in the Appendix A, the result is where Ψ(z) is the same minimum phase rational function of (11) and X (z) is the rational and stable transfer function obtained from Noting that G (z) = G(z)/z, it is simple to relate X (z) with X(z) since the only change needed is to extract the possible pole in z = 0 of W(z −1 )/z to obtain W (z −1 ). Hence, For the case we are interested in, which is tracking of acceleration steps, according to (16) G(0) = 0. Then, Therefore, the new optimum close-loop transfer function of (36) is and the corresponding optimum loop filter is now only the three poles of the input, Even more interesting is the expression of (36) in terms of the cascade of accumulators, If it could be possible to implement this loop filter, the feedback of the complete loop would be φ[i | i −1]. Of course T (z) = Y (z)z −1 cannot be implemented with a real-time loop, which has two delays, but Y (z)z −2 can. Indeed, if the loop filter structure of the UFA-PLL is slightly modified as shown in Figure 18, it can be shown that amazingly (2) i . Clearly, the delay on the feedback branch precludes the use of this value for the loop, but not for the rest of the GNSS receiver.

Two Samples Smoothing.
The previous process can be applied again. Now, the problem is equivalent to obtaining the optimal loop filter without delay. Then, if F (z) is the loop-filter transfer function and the corresponding closed loop transfer function is where Y (z) is the rational and stable transfer function to be found minimizing J(Y (z)) in (8) and replicating (36)-(39), we obtain Therefore, the optimum transfer function of (43) is where P 2 = (p 2 − 2p 3 ) and P 1 = (p 1 − 2p 2 − p 3 ). Now T (z) = Y (z) cannot be implemented with a real-time loop, but Y (z)z −2 can. Again, based on the loop structure of Figure 18 it can be shown that

More Samples Smoothing.
As it was previously mentioned, since G(z) in (16) has four zeros at z = 0 for inputs modeled as accelerations steps, the previous smoothing procedure can be done two more times. In this way, if some latency is allowed, the following phase estimates can be obtained based only on the real-time tracking loop of Figure 18-the first three equations are repeated for clarity: International Journal of Navigation and Observation Figure 18: Block diagram of the UFA-PLL model with loop filter structure modified for phase smoothing.
These phase estimates can be interpreted as a prediction in the context of estimation theory [30].
with the signal dynamic model adopted for the phase estimation. If we consider i ] as the state of the loop input phase model, the transition matrix must be Since can be found propagating backwards this state as {A −2 x i−1 } (1) . Notice that using this backward propagation process with the matrix A −1 , all the estimators of (46) can be obtained. Actually, more smoothed estimates can be built. For example, the equation can be used, but it is not optimal. As it will be shown in the simulations it can be considered useful because of its extremely simple implementation. In some way, the quantity of zeros in (16) at z = 0 gives a measure of the backward propagation capacity of the states estimated by the loop filter.

5.4.
Simulations. The simulated loop model is a UFA-PLL as shown in Figure 18 with the same filter coefficients of Section 3. The phase estimation error for an acceleration step of 50 g (starting at i = 5 and without noise) for the different estimators is plotted in Figure 19. In this situation the loop error grows up to almost one cycle and therefore the data detection during this transient will not be possible. However, applying the smoothing process described by (46) the transient error is consistently reduced each time that a new input sample is used for the estimation. The response of the suboptimal estimator of (48) is also shown. It can be seen that its transient response is slightly worse than the obtained with φ[i | i + 2]. The smoothing process also produces a decrease of the estimation noise variance. In Figure 20 the standard deviation of the six previous phase estimators is plotted for three different signal levels. These results were obtained simulating a linearized loop fed with Gaussian noise of variance 1/(2TC/N 0 ). As expected, an increase of 3 dB in the signal corresponds to a reduction of approximately √ 2 in the standard deviation. It is also possible to verify that

RF Test Experiments
The FLL-assisted PLL and the UFA-PLL designed in Section 3 were implemented in a System Developer Kit (SDK) for GPS receivers from SiRF [31]. Due to the real-time nature of this task all calculations for the loop had to be done in a fraction of 5 ms. They were programmed in fixed-point arithmetic, using some scaling and approximating coefficient values by powers of 2. Details of these implementations were given in [32]. To verify the tracking capability of the loops with real signals and without relying on expensive equipment like a GPS signal simulator, we used an RF signal generator to produce a frequency modulated carrier at 1575.42 MHz (the L1 GPS frequency). The signal was not spread with the code of a particular satellite, and thus the code generators of the GPS receiver were turned off during the test. This is not a limitation since the focus of this analysis is on the carrier loop, rather than the code loop. A triangular waveform was used as a frequency modulation to simulate steps in acceleration. The frequency deviation was selected according to the magnitude of the step (an instantaneous frequency 14 International Journal of Navigation and Observation deviation of Δ f corresponds to a λ · Δ f instantaneous velocity, where λ is the L1 wavelength). The selected carrier power was −113 dBm. Taking into account the noise of the 50 Ω output resistance of the generator and a noise figure of approximately 8 dB of the RF front-end gives a C/N 0 = 53 dB/Hz, which is a relatively high value for GNSS receivers. This value was selected to obtain low noise curves for a visual comparison with the simulated responses since the noise performance has been already characterized.
A test of 10 g acceleration steps is shown in Figure 21, depicting the tracked frequency detrended by a linear fit that accounts for the local clock drift. The amplitude of the triangular waveform was increased gradually up to the desired value to avoid large frequency steps. The measured phase error response of the FLL-assisted PLL at the output of the phase discriminator in one of the steps is shown in Figure 22, the response of the UFA-PLL is the same. The simulated response of the loop (as in Figure 8) is also displayed to appreciate that the implemented loop is properly characterized by our model.
The same experiment was performed for acceleration steps of 40 g. The phase error responses of the FLL-assisted PLL and the UFA-PLL are presented in Figures 23 and 24, respectively. In the case of Figure 23, the step presented is negative and therefore the phase response is upside down with respect to the simulation in Figure 9. Again, a fine agreement between the measurements with real laboratory generated signals and the simulations can be appreciated.

Conclusions
A new carrier tracking loop design method for real-time GNSS receivers has been presented, which is completely optimized from the perspective of the digital nature of the correlation measurements. An analysis of the phase and frequency discrimination ideas from this point of view allowed us to choose optimum discriminators often discarded because of the complexity of their analog counterparts. Also, the known structure of FLL-assisted PLL has been considerably improved leading to a carrier loop that operates normally in phase locked condition and in frequency locked condition if the dynamics become severe enough. The effect of coupling the FLL to the PLL is considered at the design stage allowing a fine control of the effective loop bandwidth. Moreover, this approach allowed us to develop the UFA algorithm that corrects the cycle ambiguities of measured phase errors using the frequency information exploited by an FLL. With this algorithm it was possible to conceive a PLL that has the same advantages of an FLL-assisted PLL but avoids cycle slips and yet is easy to implement. Regarding loop filters, their optimization was achieved directly in the digital domain. Our procedure solves the bandwidth trade-off considering the discrete nature of the filter and the unavoidable computational delays in the loops and makes possible to design loops extending beyond the restrictive constraint B N T < 0.1. This limitation is caused by the standard technique of discretizing analog loop filters, which only gives acceptable results for narrow loop designs, not suitable for high-dynamics receivers. Even more important, with our method the designed loop bandwidth is the same as the actually implemented one, since there are no approximations involved. This is of great importance considering the difficulties found in ( [5] pp. 183) with respect to the stability of designed loops. An approximate analysis of the nonlinear behavior of digital loops has been presented. The pull-out probability is approximated as an efficient tool for selecting the correlation time and the loop bandwidth so that the tracking threshold is minimized. Contrary to classical pull-out probability studies, our approach considers nonstationary scenarios as is often found in high-dynamics applications. The UFA-PLL examples presented show that 3 dB of improvement in the tracking threshold can be attained by properly selecting the integration time and loop bandwidth. It is worthwhile to emphasize that even for accelerations of 20 g, a loop sample time shorter than 5 ms is not appropriate.
It was also shown that using some state variables of the loop, smoothed phase estimates can be efficiently built with a latency of only a few samples. This can be very useful in many GNSS applications that use phase measurements, such as code measurement smoothing, differential positioning, and attitude estimation. It can also extend the use of these techniques to real-time and high-dynamics applications, where more complex phase estimation schemes are not practical and the usual tracking loop estimates do not provide enough phase accuracy. The simulations presented show that using the proposed smoothing scheme the transient responses of the loop to an acceleration step of 50 g can be almost eliminated and the estimation noise reduced by half.
Finally, our technique was used to design a loop that can optimally track steps of 20 g with a tracking threshold of C/N 0 ≈ 34 dB/Hz. It was implemented in a GPS receiver using fixed-point arithmetic and tested with frequency modulated RF signals. Simulations and experimental results confirm that our new loop designs can track the input phase in these severe conditions, with the same implementation complexity as the usual loops. This work has not considered adaptive schemes because they tend to increase the computational requirements and transient behavior. An adaptive version of the UFA-PLL could be formulated for situations where time-varying filters are affordable. Our assumptions do not deal with heavily unmodeled dynamics or disturbances on the tracked signals that justify a robust H ∞ approach. Nevertheless, given the Wiener filtering relation established for the present hypotheses, the loops obtained are equivalent to a steadystate version of Kalman-filter based-tracking loops [33]. Including robustness as one of the loop design considerations can also be addressed with the Wiener filtering formulation [27,34], and a balance between robustness and adaptiveness will be pursued elsewhere.

A. Minimizing J(Y )
Expression (8) for J can be minimized applying a standard procedure from variational calculus.
Assuming that the input noise is white with power spectral density η/2, the variance of the output noise can be calculated as where Y (e jω ) is the frequency response of Y (z) and the last integral extends over the unit circle of the complex plane.  (11) is the same as (B.3) and Ψ(z) = γrβ(z)/D(z). Moreover, (12) is similar to (B.5), and thus γrQ 1 (z)/zD(z) = X(z). We can see how both solutions match.