A Novel Residual Frequency Estimation Method for GNSS Receivers

In Global Navigation Satellite System (GNSS) receivers, residual frequency estimation methods are traditionally applied in the synchronization block to reduce the transient time from acquisition to tracking, or they are used within the frequency estimator to improve its accuracy in open-loop architectures. There are several disadvantages in the current estimation methods, including sensitivity to noise and wide search space size. This paper proposes a new residual frequency estimation method depending on differential processing. Although the complexity of the proposed method is higher than the one of traditional methods, it can lead to more accurate estimates, without increasing the size of the search space.


Introduction
In Global Navigation Satellite System (GNSS) receivers, the synchronization block is one of the most important elements. It synchronizes the received signal and demodulates the navigation message, which is necessary to compute the Position, Velocity, and Time (PVT). Traditionally, this block is divided into two parts: acquisition and tracking. The aim of the acquisition stage is to determine which are the satellites in view, and to estimate their code delay θ, and Doppler frequency f d . These raw parameters are then used to initialize the next stage, in which the satellite signals are tracked.
The acquisition can be done computing the correlation between the received signal and a series of several tentative signals generated by the receiver. The tentative signal with the highest correlation value, provided that this is higher than a predefined threshold, is assumed to be equal to the received signal. In the tracking stage, a Delay Lock Loop (DLL) is used to refine the code delays, whereas the frequencies can be refined by a Frequency Lock Loop (FLL) and/or a Phase Lock Loop (PLL). The lock loops include a discriminator to estimate the residual in the frequency and/or phase of the tentative signal w.r.t that of the received signal, and a filter is then used to reduce the noise in the residual. Then, this residual is fed back to the numerical controlled oscillator (NCO) [1]. Although PLL brings a higher accuracy than FLL, it can work well only if the residual of the initial frequency is small enough. However, the accuracy of the frequency estimation provided by the acquisition stage is such that the residual is in the range of a few hundred Hertz, which is too large an input value for the proper behavior of the PLL. Therefore, it is necessary to adopt a frequency estimation refinement method before PLL using an FLL [1], or refinements based on the results of the acquisition stage [2].
In PLL and FLL, the output measurements are fed back to re-configure the lock loops; therefore, they are considered as closed-loops. The most significant advantage of such closed-loops is their memory and their computational resource minimization [3]. However, closed-loops are vulnerable to fading effects and cycle slips [4,5]. Moreover, the possibility to extend integration time to improve the

Signal Model
The digitalized GNSS signal r[n] available in a receiver, after the Analog-to-Digital Converter (ADC) block, can be represented as: where P is the power of the signal; d[n] denotes the navigation data symbol; c[n] is the PRN code at a time instance n; n w [n] is additive white Gaussian noise (AWGN) with zero mean (µ = 0), and variance σ 2 n ; f IF and f d denote the Intermediate Frequency (Hz) and the Doppler shift (Hz), respectively; T s is the sampling period (s), which is the inverse of the sampling frequency F s ; θ is the received code delay (sample); and ϕ is the initial phase of the received signal (radiant).
The correlation value between the received signal r[n] and the tentative signalr[n] can be expressed as [2]: where m stands for the index of the coherent integration interval [(m − 1)N, mN], N = T coh F s denotes the coherent integration time T coh (s) in samples; s m , w m are the signal and the noise components, respectively. Assuming that the frequency is constant in a small interval, it can be written: where (τ = θ − θ) and ( f d = f d − f d ) are the average difference between actual and tentative code delays and Doppler shifts during an interval; and R[τ] is the autocorrelation function [17]. w m can be considered as a normal distribution random variable with zero mean and variance σ 2 w = Nσ 2 n (w m ∼ N (0, Nσ 2 n )). As pointed out in [2], there are three factors that cause loss of measured correlation value: the bit transition, the residual code in the autocorrelation function R[τ], and the residual frequency. The solution of residual code and data transition can be found in [18,19]. Moreover, some modern GNSS signals introduce pilot channels, which do not contain data bits. Although the pilot channel still has a secondary code, it can be synchronized by multi-layer code based acquisition [20]. Therefore, the data transition problem can be ignored, and, for simplicity, in the following sections, code and data are considered to be wiped off (τ = 0 and d[n + θ] = 1). The residual frequency can be eliminated by using Double Differentially (DDF) Coherent method [21], or estimated by the frequency estimators as introduced in Section 1. The DDF method works effectively for the signal detection; however, the frequency estimation is still needed for the tracking process. In this paper, only the problem of estimating f d is considered. In reality, the Doppler frequency is not constant but changes continuously overtime. This causes a drift component in the residual frequencies as analyzed in [22]. The Doppler drift also affects the correlation value and the accuracy of the frequency estimate. However, in the case of not high dynamic receiver, the Doppler drift can be ignored [22]. Moreover, the results in [7] showed that, in the high dynamic case, the convention frequency estimator can follow the Doppler changes with a proper architecture of tracking loop. Therefore, the Doppler drift is not considered in this paper.

Frequency Estimation Methods
From Equations (2) and (3), it can be seen that, in the free-noise case, the phase difference between two consecutive correlation intervals is equal to 2π f d T coh : Based on Equation (4), in order to estimate the residual frequency, several methods have been proposed in the literature and are hereafter shortly described.

Classical Kay's Method
The method, which was introduced in [8], estimates the residual frequency from the phase difference between two consecutive correlation values (see Figure 1). The estimator has the form: where the total observation time is T 0 = M · T coh (seconds), and M is the number of observed interval. This method is based on the phase wrapping algorithm. It has a low computational complexity, but it only works effectively with high SNR because it is sensitive to noise [9].

Conventional Differential Combination Method
In [16], a method to estimate the frequency that is the residual in the closest tentative signal was proposed. This method, referred to Conventional Differential Combination (CDC) method, was proved to be suitable for signals with a very low SNR. The estimated residual frequency is derived as: The difference between the CDC method and Kay's method is the order of argument function and integration (see Figures 1 and 2b). It should be noted that the argument function is a nonlinear function, while the integration acts as a linear noise filter. Therefore, if the argument function is performed before the integration, the effect of white noise in the correlation values will change, and the filtering effectiveness of the integration will be decreased.

Modified Generalized Differential Combination Method
In [2], a new method based on post-correlation combination is proposed and used to improve the sensitivity of acquisition stage. It is also called the Modified Generalized Differential Combination (MGDC) method. In this method, the correlator outputs are differently combined in different spans, as depicted in Figure 2a,c. A span-i, which is denoted as A i , is defined as the differential combination of the correlator outputs evenly spaced i intervals: Differential processing for frequency estimating. (a) differential operations; (b) conventional differential combination; (c) modified generalized differential combination; (d) modified generalized differential combination.
Considering the free-noise case, substituting the contribution of the useful signal in Equation (3) into Equation (7), it can be rewritten: Therefore, the residual frequency can be estimated by: where w i is chosen depends on the number of combinations in the respective span. The number of usable spans for the MGDC method is limited as presented in the next section. However, the results in [2] shows that: if the number of usable span is more than one, the more number of valid spans is used, the better performance of the estimation is because the more information is explored.

Wrapped Phase Problem
In all the methods above, the argument function is used to estimate the residual phase component in the correlation values. It should be noted that the argument function returns a wrapped phase, which means: Thus, if the actual phase is unwrapped (k = 0), the estimation phase error can be multiple times of 2π. To avoid the problem, the phase at the input of argument function has to be in the range [−π, π].
Before considering wrapped phase problem for each method above, it is important to note that, if the search space is defined as the parallel code phase search acquisition [23], we have: where ∆ f Step is the Doppler step size of the search space. Inequality (11) presents the variation interval of the residual frequency before the estimation. For Kay's method and the CDC method, from Equations (4)-(6), the wrapped phase constraint can be written as: Since Inequality (11) is achieved before the frequency estimation, Inequality (12) is always satisfied. For the MGDC method, the wrapped phase constrain can be written as: Consequently, the maximum number of spans in Equation (9) that can be used is constrained. For examples, if T coh = 1 ms, and f d = 125 Hz, from Inequality (13), i ≤ 4, and the combination A 5 must not be used to estimate the residual frequency. The wrapped and unwrapped phase of vector A 5 are different as shown in Figure 3. The function arg(A 5 ) does not return the unwrapped phase value (dashed angle in Figure 3), and instead it provides the wrapped phase (solid angle in Figure 3), so that including vector A 5 in Equation (9) leads to a consistent error in the estimated phase. From Inequalities (11) and (13), in the worst case-when the Doppler step is maximum-the maximum number of usable spans is 2. However, it is not recommended to satisfy sharp relation (13) since noise has unexpected influence on the computation of the argument function because, in such conditions, the noise could lead to consistent errors in the phase determination. It follows that only span-1 should be used on this case that, obviously, is equivalent to the CDC method. The only ways to increase the number of usable spans is to reduce the Doppler step in the search space. However, this leads to an increase in complexity, especially when a long integration time is used.

Proposed Method
Recall that, in the MGDC method, higher order spans returning a wrapped phase different from the unwrapped one (as it is the case of span-5 in Figure 3) can not be used for the phase estimation according to Equation (9). To overcome this problem, note that, in the ideal case, the phase difference between two consecutive span-i vectors is equal to 2π f d T coh regardless of the index of the spans: (14) where . The index of used span no longer appears in Equation (14). This allows for estimating the residual phase using A 1,i and without being affected by the unwrapped phase problem. A i can be considered as the first order differential of correlations, and A 1,i can be considered as a second order differential. From Equation (8), for a specified number of correlation values, the higher span order would use less number of correlation values, which leads to smaller magnitude of A i and A i . In other words, the magnitude G i represents the amount of information in each combination, and thus can be considered as a weighted value in a new frequency estimator (see Figure 2a,d), which follows: where 2 ≤ K ≤ M − 1 is the number of combinations. A large value of M is not recommended because it causes a great delay in the tracking process and reduces the accuracy instantaneously in the case of high dynamic users. Therefore, in this paper, M is chosen as equal to or less than 200.
Since the unwrapped phase problem is avoided, the number of combinations in this estimation method is not limited as it is in the MGDC method. Hence, this method is more effective than the others in case of a large residual phase, and does not require an increasing of the frequency domain size in the search space. Especially in weak signal scenarios, where adopting the longer T coh is a convenient approach to increase the sensitivity of the synchronization block, smaller Doppler step leads to a noticeable increase in the acquisition/batch processor complexity, which will be shown in the next section.
Another important point about the derivation of relation (15) is the order in which argument and sum functions are computed that can affect obtained results, as it has been already discussed in Section 2.2.2. The authors have first implemented the sum that, acting as a linear filter, reduces the effect of the noise before the nonlinear argument function is computed.

Theoretical Analysis
This section presents the probability density function of estimated frequency of CDC and MGDC methods. To the best of the author's knowledge, the theoretical analysis for the Kay's method and new MGDC method are not easy to derive. However, their performance can be assessed through the Monte Carlo simulation, as introduced in Section 5.
According to [24], a complex random variable X with the real part (X) and image part (X) are independent and follow the normal distribution: , has joint probability density function (pdf): By using the polar coordinates, (16) is transformed into: where: r = |X|, φ = arg(X), and φ ∈ [−π, π]. Therefore, by integrating (17) over r, we obtain the pdf of φ as follows [16]: where: Substituting φ = 2πT coh f into Equation (18), we have: Equation (22) can be applied to compute the probability of frequency estimated from the complex normal random variable X. However, applying Equation (22) for Kay's method is inappropriate because R m R * m−1 is a product of two normal random variables, which does not follow the normal distribution. The pdf of R m R * m−1 is quite complex as computed in [25], and it is not easy to derive p f ( f ) from that. For CDC and MGDC methods, if the number of combination M is high enough (M ≥ 20), applying the Centre Limit Theorem (CLT), A i in Equation (7) can be considered as following the normal distribution (for i = 1, MGDC is equivalent to CDC). Therefore, Equation (22) can be used to obtain the estimated frequency distribution of CDC and MGDC methods.
Applying Equation (22), the expected values and variances of the real and the image part of A i are needed. From Equations (3) and (7), we have: From Equation (3), s k and s k−1 are deterministic, so: As proved in [22]: A i,2 can be represented as: The components A i,21 and A i,23 can be analysed by the same way, as follows: Then: The in-phase component of A i, 22 can be computed as: Similarly, the quadrature component of A i,22 is: Consequently, the variance of the real part and image part of A i can be written as: From Equation (24), the mean values of the real part and image part of A i are: For the CDC method, applying Equations (31) and (32) with i = 1 into Equation (22), the pdf of the frequency estimated by CDC method is obtained. Then, the mean of the f CDC can be computed as: and the variance of the estimated frequency is: Figure 4 shows the theoretical line of p f ,CDC ( f ), and the normalized histogram of f CDC obtained by Monte Carlo simulation with 500,000 simulation runs. The input residual f = 0 Hz in Figure 4a, and f = 250 Hz in Figure 4b. The matching between pdfs and the histograms validates the theoretical analysis above.  For the MGDC method, because the combinations A i are correlated, it is not easy to compute the pdf and variance of estimated frequency. However, the pdf of the estimated frequency at each span can be calculated by applying Equations (31) and (32) into Equation (22), with i being the index of the considered span. Then, the mean value of f MGDC can be calculated as:

Performance Analysis
In this section, the accuracy and the complexity of the proposed method is assessed and compared with the methods introduced in Section 2.2 using Monte Carlo simulation for 50,000 test points. Regarding the accuracy of the evaluated estimators, two signal conditions are considered. The first one is the nominal range of C/N 0 , which is 35 dB-Hz and above [18,26], and the second one is the weak signals, in which C/N 0 is in the range of [20 ÷ 35] dB-Hz. In the nominal conditions, the integration time equals 1 ms and the number of combinations M = 20 are adequate to detect and track the signal [7]. While the signal is weak, two possible solutions can be adopted: extending the integration time, and/or increasing the number of combinations. In this section, to guarantee the proper complexity, the used number of combinations M does not exceed 200, the number of combinations for the second order differential is K = 20, and the integration time is extended up to 10 ms. The residual frequencies are simulated with considering Inequality (11). Figure 5 shows the simulation results of mean errors and standard deviation estimations of the methods when the residual frequency f d = 100 Hz and T coh = 1 ms. From Inequality (13), i ≤ 5, or the number of usable spans in MGDC method is up to 5. However, the simulation results in both mean errors and standard deviation estimations of MGDC span-1, 2, 3, 4 show a better performance than MGDC span-1, 2, 3, 4, 5. This is because when the equal sign occurs in Inequality (13), the idea vector A 5 is matched with the boundary between −π and π. This leads to an ambiguity in the estimated phase of vector A 5 , which can jump from −π to π under noise effect, and causes errors in the estimation. The same result can be found in Figure 6, which shows the results of worst case, when f d = 1 4T coh = 250 Hz. Although the theoretical analysis in Section 2.2 indicates that the maximum usable spans for the case is 2, the simulation results of MGDC span-1, 2 show a worst performance. As a result, the equal sign in Inequality (13) is not recommended. In Figures 5 and 6, the new MGDC method gives the best performance in both mean error and standard deviation. This can be explained because the number of combinations in the new MGDC is higher than those of the other methods. In Figure 6, since the residual frequency is equal to 1 4T coh , the number of usable spans in MGDC method should be 1 as analyzed in Section 2.2.3. Thus, MGDC and CDC are congruent. While the MGDC method cannot improve the performance of estimation in comparison with CDC, the new MGDC reaches this goal with the best mean error and standard deviation.

The Accuracy of Estimations
In both Figures 5 and 6, the performance of Kay's method reduces rapidly when the C/N 0 < 40 dB-Hz. This is compatible with the analysis in Section 2.2 about the nonlinear effect of argument function on Kay's method. Thus, in the rest of this part, which focuses on weak signals, the method will not be considered.
In the case of weak signals, the solutions to improve the performance of the receiver are increasing the number of combinations and/or using a longer integration time. The first solution is preferred in terms of effective implementation. In this scenario, the number of combinations for the first order differential is chosen to be 50, 100, 200. The number of combinations for the second order differential is K = 20, this is chosen to guarantee a proper complexity. In Figure 7, because the residual is 100 Hz, even span-4 can be used for the MGDC method. It can be seen in Figure 7 that the estimations provided by the proposed method are always more accurate than those of the CDC method. Similar to Figure 5, for several low C/N 0 signals, although the traditional MGDC gives smaller estimation variances, its mean errors are larger than those of the proposed method. When the number of the second order differential is up to 200, the new MGDC gives the best performance.
In Figure 8, when the residual is 250 Hz, only span-1, which makes MGDC coincide with CDC, is usable. The results have the same evolution with Figure 6, where the performance of the proposed method is much better than the others. It can be seen from Figures 7 and 8 that the proposed method can provide an accurate estimation even in the case of low C/N 0 . The estimation error is about tens of Hertz, which is suitable to pass into a PLL. In general, with the same number of used correlation values, the new MGDC is much more accurate than CDC. In comparison with MGDC, the improvement of the new MGDC depends on the number of combinations and the power of signal. However, it should be noted that, in order to apply MGDC, a greater frequency search space is required. It leads to the increase of computational load. This issue will be discussed more in Section 5.2. For the second solution, which uses longer integration time, from Inequality (11), the amount of frequency error at the input of estimation will be much smaller. In Figure 9, the integration time T coh = 10 ms; then, the residual frequency must be smaller than 25 Hz. To avoid the equal sign in Inequality (13), and possibly to compare the performance with the MGDC method, we chose f d = 20 Hz. The estimation results were shown in Figure 9. Since f d = 20 Hz, the maximum number of used span is 2. Thus, MGDC span−1,2 gives a better performance than CDC does, while the MGDC span−1,2,3 gives the highest error (more than 10 Hz, which is half of input error). At the same time, the new MGDC standard deviation of estimation error is the best and smaller than 1 Hz for all of the investigated signals.

The Computation Complexity
This part considers the computation complexity of the above methods. For Kay's method and CDC method, their complexity can be found easily from (5) and (6). The summary of computation complexity is shown in the Table 1. The calculation for the complexity of the proposed method is described below. Since the number of correlation values is equal to M, the numbers of multiplications (α) and summations (β) are: The complexity of traditional MGDC can be derived in a similar manner.
In MGDC and the new MGDC method, since K ≤ M − 1, the number of operations is maximum when K = M − 1. Figure 10a shows the maximum number of multiplications of the proposed method, when M is from 1 to 200, and the number of multiplications of the CDC and Kay's method. Obviously, the complexity of the proposed method is the highest and it increases when the number of combinations grows up. However, it allows decreasing the search space in the batch processor of open-loops, and acquisition stage of closed-loops architecture, therefore reducing the overall computation burden. In more detail, from Inequality (11), the maximum frequency step is: In the case of high dynamic users, the Doppler frequency is around ±10 kHz. Thus, the minimum size of frequency search space is: Obviously, the longer the integration time, the greater the size of frequency search space. The load of calculations for the correlators using Fast Fourier Transform (FFT) at each interval is given by [27]: where N = T coh F s is the number of samples for one integration interval. From Equation (40), the change of computational load of acquisition stage (or batch processor) versus the Doppler step is shown in Figure 10b. Comparing with the maximum number of operations of the new MGDC method in Figure 10a, it can be seen that the computational burden of acquisition stage is much higher than that of the frequency estimators. Therefore, to save on the computational cost, reducing the search space is more effective than choosing a low cost frequency estimator.
To apply the MGDC method, the size of search space must be increased depending on the number of spans we want to use. For example, using up to MGDC span-2 requires the double density of search space. Similarly, MGDC span-3, 4, 5 requires the frequency grid to be three, four, or five times denser than CDC and the new MGDC do. In more detail, assume that the sampling frequency is 4 MHz and the integration time T coh = 1 ms; then, N = 4000 (samples) and the minimum number of searched frequencies N F = 41. At each interval, the number of multiplications and summations are 4.8042 × 10 6 and 7.5082 × 10 6 , respectively. For each added frequency in the search space, the number of calculations increases 5.2383 × 10 5 . This number will be much greater if a longer integration time and/or a greater sampling frequency are used. Considering Equations (36)  Another considerable example is the particular cases in Figures 7 and 8, in which the standard deviation of estimation error ( f ) is about 50 Hz at C/N 0 = 26 dB-Hz. The methods and their parameters are summarized in Table 2. It can be observed that, while the difference between the complexity of CDC and MGDC is not much, the number of multiplications and summations of new MGDC is so high. However, the key difference here is the frequency error ∆ f d at the input of each method. In normal cases, if the input frequency error is ∆ f d , from Equation (39), the maximum frequency step in search space is: ∆ f Step = |2∆ f d |. Therefore, the minimum value of N F must be: Hence, in the considering condition (C/N 0 = 26 dB-Hz, var( f ) = 50 Hz), the acquisition stage requires a minimum frequency space size of 101 for CDC and MGDC span-1,2,3,4 , and 41 for the new MGDC. Thus, the load of calculation in the acquisition stage can be reduced to less than a half if the new MGDC is then applied.

Conclusions
In this paper, a new method to estimate the residual frequency for GNSS signal processing is proposed. This new method gives a higher accuracy than the traditional method, such as: Kay's, CDC, and the MGDC method. Moreover, although based on the idea of MGDC, it is not constrained by the amount of residual frequency input as the MGDC method. The complexity of the proposed estimation method is higher than the others; however, it reduces considerably the complexity of the acquisition stage. The proposed method is suitable to be adopted to the frequency estimator in open-loop tracking architecture, or applied in closed-loop to reduce the transient time from acquisition to tracking.
In future works, the effect of data transition and Doppler drift will be considered. A full open-loop architecture, which uses the new MGDC method to estimate frequency, is going to analyze in both theoretical scenarios and simulation.