Fast Nearly ML Estimation of Doppler Frequency in GNSS Signal Acquisition Process

It is known that signal acquisition in Global Navigation Satellite System (GNSS) field provides a rough maximum-likelihood (ML) estimate based on a peak search in a two-dimensional grid. In this paper, the theoretical mathematical expression of the cross-ambiguity function (CAF) is exploited to analyze the grid and improve the accuracy of the frequency estimate. Based on the simple equation derived from this mathematical expression of the CAF, a family of novel algorithms is proposed to refine the Doppler frequency estimate with respect to that provided by a conventional acquisition method. In an ideal scenario where there is no noise and other nuisances, the frequency estimation error can be theoretically reduced to zero. On the other hand, in the presence of noise, the new algorithm almost reaches the Cramer-Rao Lower Bound (CRLB) which is derived as benchmark. For comparison, a least-square (LS) method is proposed. It is shown that the proposed solution achieves the same performance of LS, but requires a dramatically reduced computational burden. An averaging method is proposed to mitigate the influence of noise, especially when signal-to-noise ratio (SNR) is low. Finally, the influence of the grid resolution in the search space is analyzed in both time and frequency domains.


Introduction
The main purpose of the acquisition and tracking systems of a Global Navigation Satellite System (GNSS) receiver is to provide an estimate of the Doppler frequency f d , the code delay τ , and the phase of the carrier, φ, of the signal transmitted by each visible satellite. The task of the acquisition system is to detect the visible satellites and to provide, for each detected satellite, a coarse estimate f a d ,τ a of f d and τ . This parameter vector is then passed to the tracking systems, whose task is to refine this estimate. The refinement of Doppler frequency estimate is generally performed by a classic phase lock loop (PLL), which requires an initial estimate much more accurate than the one provided by the acquisition system. Therefore it is necessary to improve the accuracy of the estimatef a d to an acceptable level before starting the operations of the phase tracking loop. A system typically adopted by a GNSS receiver to reach this goal is a frequency lock loop (FLL), which is generally integrated within a PLL. The first refinement is done by a robust FLL operating at wideband, then the loop bandwidth is gradually reduced and finally the system switches to a PLL scheme [1,2]. Other methods [3,4] refine the frequency estimate by exploiting the phase difference between two successive periods of data. An interpolation method is introduced in Reference [5] to estimate the true value of the Doppler frequency, but it is based on an empirical approximation.
In most of the previous methods, usually the estimates of f d and τ are picked in a search plane only considering the peak cell without any usage of the other cells. In the fields of communications, audio, medical, instrumentation, and others [6], the problem of estimating the frequency of a tone contaminated with noise is tackled for example by Quinn [7,8], MacLeod [9,10], and Jacobsen [6,11], by exploiting the idea of refining the final frequency estimate using the peak sample and two neighbors of the discrete Fourier components. At the same time there are other methods, studied in Reference [12,13], which utilize the phase information. These methods cannot be directly applied to the acquisition of a GNSS signal, because of the very low signal-to-noise ratio and the different signal model, but they can inspire us to do some innovation in GNSS frequency estimation.
In this paper, the peak and neighbor points of the cross-ambiguity function (CAF) in the frequency domain are used to derive a simple formula that greatly improves the accuracy of the frequency estimate provided by the acquisition system. The CAF was initially derived in Reference [14] using statistical principles, then Reference [15] presented a new approach of the CAF derivation. Furthermore in this paper the approximation in CAF main lobe is analyzed in details, based on this approximation, and a new family of methods for refining the estimate of the Doppler frequency is proposed, which exploits the cells close to the peak in the search plane. Compared with the traditional methods, these methods significantly improve the accuracy without increasing the computation complexity or using additional received data.
A preliminary version of this work was presented in Reference [16]. With respect to that previous work, here we extend and complete all the mathematical derivations, extend the performance analysis with appropriate comparisons, derive and discuss the Cramer-Rao lower bound (CRLB) for the frequency estimator showing that the proposed approach is close to the CRLB (quasi-ML approach), and include the theoretical analysis of other non-AWGN nuisances. This paper is organized as follows. In Section 2 the signal model is presented and the approximate mathematical expression of the CAF in the main lobe is obtained. In Section 3 a new family of algorithms is derived and proved to work perfectly in the absence of noise. In Section 4 the performance of the proposed algorithms is investigated in the presence of additive noise; both the CRLB and a least-square (LS) solution are derived as benchmark, and the comparison shows that the new algorithms can approach very closely the CRLB. Besides that, a simple averaging approach based on non-coherent sums is proposed to improve the accuracy of the algorithms in low SNR conditions. Furthermore, in Section 5, the effects of other nuisances, uncorrelated with the additive noise, are analyzed, and some countermeasures are proposed. Finally, in Section 6 the conclusion is drawn.

Fundamentals of the New Algorithms
The basic scheme of the acquisition method proposed in this paper is illustrated in Figure 1. The left part of the figure indicates the traditional GNSS acquisition process from which a two-dimensional search grid (marked in green color) is generally obtained, while the right side shows the presence of a new additional block able to refine the Doppler frequency estimate in a simple way. The innovation proposed in this paper refers to the algorithms used by this additional block to refine the frequency estimate.

Acquisition Process
The acquisition system for GNSS application is based on the maximum-likelihood estimation theory, which can be briefly described as follows [17].
The incoming sampled signal can be denoted as a vector where L is the total number of the samples, and where r(n) is a signal that contains a vector of unknown parameters a = [α 1 α 2 α 3 · · · α K ], W (n) is a zero-mean White Gaussian Noise (WGN) random process with variance σ 2 , and 0 ≤ n ≤ L − 1.
The ML estimate of the parameter vector a can be found by maximizing the likelihood function, which depends on the probability density function (PDF), that is where the test signal r(n) has the same structure of r(n), but the unknown parameter vector a is substituted by a vector a, whose elements are variables defined in a range ran a that contains all the possible values of the unknown vector a, that is to say, a ∈ ran a .
If the energy of the test signal r(n) (that is the term (3)) does not depend on a, then it is possible to show [17] that the corresponding ML estimate a ML of a is So, in this case, the ML estimation actually depends on the scalar product R(a) between the test signal and the received signal, defined as a ML can be found by searching the maximum R(a) in the range ran a . In GNSS field, without considering the influence of noise, the received signal, after down-conversion and sampling, can be written as [1] y (nT s ) = Nv m=1 y m (nT s ) (6) where N v is the number of satellites in view, and where A m is the amplitude of the signal, From Equation (7), we can learn that in principle, the satellite signal actually contains four unknown parameters: code delay (τ ), Doppler frequency (f d ), carrier phase (ϕ) and data bit. However in the acquisition process, only two of them are estimated, which are τ and f d .
With respect to the parameter data bit, in the implementation, a non-coherent acquisition scheme is used to solve the problem, so here we assume that there is no data-transition in the accumulation period.
Considering the parameter carrier phase, its influence can be removed by involving two components in the acquisition process that are in-phase component (I) and 90 • phase-shifted quad-phase (Q) component [19]. Therefore, the test signal r(n) can be written as where the parameter vector becomes a = τ , f d , and the energy of r(n) is not related to a. So the accumulation process in acquisition can be expressed mathematically as Equation (9) is known as cross-ambiguity function (CAF). Based on Equations (4) and (5), the ML estimate of [τ, f d ] can be obtained [17], as where |·| means the modulus of a complex value, and the range of a, ran a will be discussed in Section 3.
There are mainly three acquisition schemes [19]: serial search acquisition, parallel frequency space search acquisition and parallel code phase search acquisition. No matter what kind of scheme is used, a two-dimensional search grid ( Figure 4) is always obtained, and the resulting estimated vector is selected as the location of the peak cell, and, at the same time, the other cells in the search grid are abandoned. However, because of the large frequency searching step f sp , the frequency estimate error is located in the range [− fsp 2 , fsp 2 ], so the initial Doppler frequency estimate is usually not accurate enough to pass to the tracking loop directly.
In order to refine the Doppler frequency estimate, a system typically used is the Frequency Lock Loop (FLL) mentioned in Reference [1,2]. An FLL needs additional data and a special structure, which is generally embedded inside the tracking loops. Another typical technique [3,4], which exploits the phase relation of consecutive data (p.150 in Reference [3]), also needs additional data, and, at the same time, encompasses an ambiguity problem in the phase measurement that has to be solved. Actually this technique is essentially similar to an FLL with a particular discriminator.
In this paper, we develop new methods to refine the Doppler frequency estimate, based only on the search grid already evaluated by the acquisition; that is, we do not have to compute new correlations, but we only use the neighbor cells of the CAF peak, already available in the search grid.

Analytical Expression of the CAF
The CAF [15] is used in radar, sonar and other similar systems to estimate the time delay and the Doppler shift of an incoming signal. An accurate estimation of these signal parameters generally requires the evaluation of several CAF samples, at the cost of an increased computational complexity. In this paper we propose a family of methods that exploits the knowledge of an approximate expression of the analytical formula of the CAF, given in Reference [15], to reach a trade-off between accuracy and complexity.
Following the approach presented in Reference [15], the CAF associated to the generic i-th satellite code, locally generated for each trial value of code delay τ and Doppler frequency f d , can be written as where R m,i (τ , f d ) is the contribution to the CAF of the m-th signal y m (nT s ). Its analytical expression is [15] where B m = (A m /2T d )e −jϕm , the subscript i denotes the i-th satellite code generated by the local generator, T d is the integration time, F {} denotes Fourier transform, the symbol * denotes convolution operation, T p is the code period, τ is the code delay estimate introduced in the local code, Sinc(x) = sin(x)/x, f d is the Doppler frequency estimate introduced in the local carrier, Δf is the Doppler frequency estimate error, which can be expressed as and P T d is a window function defined as Since b m,i (t) is a periodic signal with a period equal to the code period, its Fourier transform leads to a line spectrum with coefficients given by (16) and the convolution with the line spectrum leads to the summation in Equation (12). In Figure 2 the distribution of a k(m,i) is shown in the case τ m − τ = 0.2 T ch , where T ch is the chip duration. Thanks to the property of Pseudo Random Noise (PRN) code, as expected, a 0(m=i) predominates over the other a k(m,i) . To better understand the nature of the summation in Equation (12), let us refer to Figure 3 where for simplicity, we represent the quantity to demonstrate the relationship among different Sinc functions. As Figure 3 shows, if a 0 is the peak component in Equation (12) (the subscript i, i is omitted to simplify the notations) only the components strictly adjacent to a 0 (i.e., a −2 a −1 a 1 a 2 ...) affect the shape of the main lobe of Equation (12), while the contribution of faraway components can be ignored. So, if we can guarantee that the adjacent coefficients (a −2 a −1 a 1 a 2 ...) are far smaller than a 0 , the mathematical expression of the CAF in the main lobe (subscript "ml") can be written as where a 0 = a 0(i,i) = R(Δτ ), and Δτ = τ m − τ . As a conclusion the approximate expression of the CAF in its main lobe is  The validity of this approximation can be improved in two ways: 1. By enlarging the integration time T d . In fact as T d increases, the adjacent components will "move" away relatively. In other words, as T d increases the width of the lobe decreases, while the distance between two sinc functions stays constant, as it depends on the code period. 2. By decreasing the values of the adjacent coefficients (a −2 a −1 a 1 a 2 ...). This can be obtained by improving the accuracy of code delay estimate, so as to work close to the maximum of R(Δτ ).
Based on the CAF expression in Equation (19), new algorithms for a better estimation of the Doppler frequency are discussed hereafter, both in ideal (i.e., noiseless, Section 3) and realistic (Section 4) scenarios.

Doppler Frequency Evaluation in the Absence of Noise
The signal acquisition process is basically a two dimensional search in a grid plane (commonly referred to as search space), as shown in Figure 4, where τ ∈ (0, T p )(X-axis range), and The variables under test τ and f d are discretized with a step τ sp for the code delay, and a step f sp for the Doppler frequency. The integration time is T d = LT s (where L is the total number of integrated samples). The number of trial points in the two axes are N τ = T d /τ sp , and N f = 2f dmax /f sp . Therefore the grid plane contains N τ × N f cells, and each cell (marked by yellow color in Figure 4) corresponds to a parameter pair τ , f d . Finally the decision variable for the acquisition is The purpose of a traditional acquisition system is to find the coordinates of the peak cell of the grid plane when the satellite we want to detect is visible. To improve the accuracy of the estimates, the steps τ sp and f sp must be decreased, at the expenses of the computational complexity, since the number of points of the search space increases. The empirical value f sp = 2/(3LT s ) is a typical choice [8] for the Doppler frequency step.
An example of search space is shown in Figure 4. The column (marked in green) crossing the peak cell (marked A) contains cells that share the same code delay. When the i-th satellite is visible, the function in this column is as shown in Figure 5, where the X-axis contains the variable Doppler frequency f d , while the Y-axis represents the absolute value of the acquisition test statistic S(τ , f d ) ml . As indicated in Figure 3, the width of the main lobe is 2/T d = 2/(LT s ). Since f sp = 2/(3LT s ), this guarantees that three adjacent points of the frequency domain CAF (A, B, C) are located in the main lobe, and then we can assume that Equation (19) is always valid in these points.  where the code delay is the same for all the cells in the same column (τ A = τ B = τ C ) and we adopt the notation S X = S(τ X , f d X ) ml , X = (A, B, C) for simplicity. In the case of no noise, these triplets can be used to find the true Doppler frequency f d , as it will be shown hereafter. This will be also the starting point of the estimation method proposed in this paper, when the measurements are affected by noise.

True Solution Based on the Absolute Value of CAF
In Reference [16], which is our initial work on this topic, we can also find the basic idea about the solution based on the absolute value of the CAF. Under the assumption that there is no data transition in the integration interval, based on Equations (19) and (20) we can write the following equations [16] where the unknowns are the amplitude A, the value of the code correlation function R(Δτ ), and the true Doppler frequency f d . We are going to show that the value of f d can be easily computed from the above system of equations [16].
First of all we observe that, when the points are located in the main lobe (like points A, B and C in Figure 5), we can write Sinc(πLT s (f d − f A )) = Sinc(πLT s (f d − f A )). Now, according to Equation (13) Based on Equations (21) and (22), we can write where Considering f sp = 2/(3LT S ): , from which we can write the true Doppler frequency as This expression gives the correct value of the Doppler frequency f d , independently from the code delay error Δτ , in the absence of noise and other nuisances. This equation represents a weighted average of the three measured points, and can be used as a first promising estimator of the Doppler frequency even when the CAF is affected by noise. Notice that Equation (27) is valid only when the Doppler frequency step is f sp = 2/(3LT S ). More in general, if we choose the Doppler frequency step as where n is a positive integer, we can choose the cells as follows: • When n is odd, we take (n − 1)/2 cells at each side of the peak cell, as shown in Figure 6(a).
• When n is even, we take n/2 cells at one side and [(n/2) − 1] cells on the other side of the peak cell, as Figure 6(b) shows.
No matter what value n assumes (even or odd), there will be n cells taken into the final calculation, and Equation (26) will become where Δf 1 = f d − f 1 and f 1 is the trial value of the Doppler frequency in the first cell of the set (i.e., the cell 1 marked in Figure 6). Then Equation (27) becomes which can be adopted as an estimator of the Doppler frequency in the presence of noise.

True Solution Based on Complex Values of CAF
Based on References [7,10], we can obtain another similar solution by using the test statistic S(τ , f d ) ml in the complex form given in Equation (19). In fact it is possible to show that At the same time, similarly to Equation (30), we can generalize this solution as where f k+1 = f k + f sp and f 1 is the Doppler frequency value in cell 1 marked in Figure 6. For simplicity, hereafter we refer to solution Equation (31) as algorithm C-3 and solution Equation (32) as algorithm C-n.

Test of Validity
The formulas in the previous sections show that, in the absence of noise and other nuisances, the proposed equations are able to evaluate the true Doppler frequency, hence reducing the estimation error range from the traditional (−f sp /2, f sp /2) to theoretically zero. A possible residual error can arise due to the fact that the method is based on an approximate formulation of the test statistic Equation (19). Therefore to test its validity, we set up a simulated acquisition campaign, where the Doppler frequency estimation error due to the traditional acquisition method is compared with the residual error introduced by the algorithms R-3, R-n, C-3, and C-n. In the simulations, the GNSS signals are generated by using the signal simulator N-FUELS [20], and several instances of a Galileo E1-b signal are obtained. Firstly we tested the methods in an ideal scenario (i.e., noiseless), with parameters f IF = 4 MHz, f s = 17 MHz, τ = 0.11 ms.
The accumulation time of the acquisition stage is set to the minimum period (4 ms), assuming that no data transition occurs in this period, the Doppler frequency step is and 67 different values of f d are randomly chosen in the range (−4, 500, 4, 500) Hz. This makes the original errors uniformly distributed in (−f sp /2, f sp /2), as also proved in Figure 7, where the cumulative distribution of the Doppler frequency estimation errors is shown. After the refinements are obtained by using R-3 or C-3, the error range decreases to nearly (−0.8, 0.8) Hz, which is a residual numerical error due to the approximations introduced in Equations (18) and (19). Therefore, we can conclude that in the absence of noise, in experiments, the new methods eliminate the error in the evaluation of the Doppler frequency due to the discretization of the search space, just using three cells selected in the main lobe of the frequency-domain CAF. The only constraint is that the Doppler frequency step has to be as given in Equation (28), with n = 3. Moreover the obtained results show that the method is not affected by the code delay error Δτ .
The practical situation in which the noise is unavoidable is discussed in the next section.

Doppler Frequency Estimation in the Presence of Noise
In real scenarios, we have to consider the influence of noise, which can be modeled as an Additive White Gaussian Noise (AWGN), as commonly done in the literature. Then the received signal model becomes where η(nT s ) represents the white Gaussian noise normally distributed with zero mean and variance σ 2 IF , related to the power spectral density S N (f ) = N 0 /2 of the analogue noise by the well-known formula valid when the transfer function of the equivalent front-end filter is assumed flat over the whole digitization bandwidth (−f s /2, f s /2). Without considering the data-transition, the in-phase and quadrature components of the CAF for the local parameters f d and τ can be written as which, based on Equation (19) (and omitting the subscript i), becomes where Real{} and Img{} mean respectively real and imaginary part of a complex value, and According to Reference [21], we know that N I and N Q are still white noise processes with variance var(N I ) = var(N Q ) = σ 2 IF /(2N ), and the envelope of S(τ , f d ) = I + jQ is where the term w I +w Q = w is a random process with mean E{w} = E{N 2 I +N 2 Q } = σ 2 IF /N including all the noise contributions. At this point, working as in Equations (27) or (31), we can develop two new estimators in the presence of noise, that is and (for simplicity the Doppler frequency step is as in Equation (33)). Hereafter, we refer to Equation (40) as algorithm "R n -3" and to Equation (41) as algorithm "C n -3". As before, they can be generalized to "R n -n" and "C n -n".
In the following two subsections, we discuss two terms of comparison worth to be considered for the frequency estimators proposed so far, firstly a least squares solution and secondly the CRLB on the variance of the estimator. Performance comparisons obtained in simulation are presented in Section 4.3.

Least-Square Method
Another approach to exploit the CAF points S A , S B , S C as defined in Equation (39) is to set up a least squares (LS) problem as follows.
Since we know that the main lobe of the CAF is a sinc function, it is possible to use the LS method in which the fitting curve is the sinc function Equation (21). The LS method requires that the sum of the squared residuals is minimized with respect to the unknown parameters, whereĀ = A/(2R(Δτ )). This leads to the equations Since we are interested in the performance comparison with R n -3 or C n -3, we will use a numerical method to solve Equation (43) for f d .

Evaluation of the Cramer-Rao Lower Bound
Usually, a statistical estimator is characterized by its bias (mean error), variance (mean square error), and the threshold SNR (signal-to-noise ratio) [22]. Hence, here the Cramer-Rao lower bound (CRLB) is proposed as benchmark to compare the performance of different algorithms.
Without considering the data-transition problem and ignoring the influence of the code delay estimation error (which will be discussed in Section 5.2), the received samples Equation (34) can be mathematically expressed as which contains the unknown parameter vector θ = [A, f d , ϕ]. The corresponding Fisher information matrix I for an observation of N samples has elements [22] [ where y(n; θ) = A cos(2π(f IF + f d )nT s + ϕ), and N is the number of samples used for the estimation of the unknowns. Therefore the Fisher matrix is The CRLB is found as the [i, i] element of the inverse of the Fisher matrix: Therefore the CRLB for the Doppler frequency estimate is where ρ = A 2 /(2σ 2 IF ).
As N = f s T d 1, Equation (48) can be written as which is the benchmark for our estimation algorithms.

Simulation Experiments for Performance Assessment
To test the algorithms R n -3, C n -3 and LS, we performed several simulation experiments with different values of the carrier-to-noise density ratio (CNR), corresponding to a signal-to-noise ratio ρ given by: where B is the one side front-end bandwidth, assumed ideally flat over the whole digitization bandwidth. The parameters used in the experiments are f IF = 4 MHz, f s = 17 MHz, τ = 0.11 ms, and f sp is set as in Equation (33). We use the root-mean-square error (RMSE) computed for the different algorithms as a metric of performance comparison: where f d is the Doppler frequency estimate, and E{·} is the expected value (estimated as a temporal average along the simulation runs). The metric f RM SE is calculated for each algorithm and compared with the square root of CRLB.
In the first group of experiments, we set T d = 4 ms and we executed 1000 runs for each CNR. Then, according to Equation (51), we calculated the corresponding RMSE. The results can be seen in Figure 8. We can see that in this case the algorithm R n -3 is better than C n -3, as it achieves a lower RMSE closer to the CRLB. At the same time we can see that the algorithm R n -3 is very close to the least-square method, though the latter is slightly better. However, considering the computation complexity of the LS method, the algorithm R n -3 appears really competitive with respect to the LS. Finally, the threshold CNR (below which the RMSE rapidly worsens) is around CNR = 38 dB-Hz in all the proposed methods.   In the second group of experiments, we changed the integration time to T d = 8ms. Similarly, we executed 1,000 runs for each CNR and obtained the results reported in Figure 9, where we can observe that, as expected, the RMSE is lower than in Figure 8, since a longer integration time reduces the effects of noise. Again, R n -3 is the closest to the least-square method, and the threshold CNR decreases to around CNR = 35 dB-Hz in all three proposed methods.  Comparing Figures 8 and 9 we can observe that, first, at the same CNR, the RMSE in Figure 9 is decreased by about a factor of 2.8 with respect to Figure 8. This is in agreement with the theoretical CRLB given in Equation (49). In fact when T d = 4 ms is replaced with T d = 8 ms, the CRLB bound decreases by a factor of √ 8. Furthermore, when CNR is relatively high (like CNR≥ 45dB-Hz in Figure 8), the proposed three algorithms are very close to the CRLB and each other.

Averaging Method
Since noise is dominant in the acquisition process, in order to increase the robustness of the proposed approach in the presence of noise, the performance of a simple averaging method, based on the idea of non-coherent sums, is assessed here. The main steps of the method represented in Figure 10 are: 1. Find the initial code delay using a first period of data, then pick out the column (marked in blue color in Figure 4), and save (2J + 1) points (J points at each side of the peak point as illustrated in Figure 10  In the following simulation, we set CNR = 43 dB-Hz, n = 3, J = 2, M = 4, and we execute 1,000 independent runs, for both averaging and non-averaging strategies. The result is shown in Figure 11, where we can see that the new averaging method decreases the error range from (−25,25) Hz to nearly (−15,15) Hz.

Analysis of Other Non-AWGN Nuisances
The performance of the algorithms presented so far depends not only on the additive noise but also on other nuisances, whose impact is analyzed in this section. In particular we observe that the methods R n -3 and C n -3 are based on the two measured vectors (27) and (31), obtained by reading the peak cell and the adjacent cells in the same column (marked as C, A and B in Figure 4). In particular, the code delay τ p is kept constant in these three cells.
Since the search space is discretized, in general even in the absence of noise f p = f d , and τ p = τ . This can be seen as a quantization error Δf q = f p − f d in the frequency domain, and Δτ q = τ p − τ in the time domain. We know that Equations (27) and (31) completely eliminate the quantization errors in the absence of noise. However, in the presence of noise, we experience an accuracy degradation due to the noise influence on vectorS, as shown in the simulation results of the previous section. The purpose of the following analysis is to state the influence of the frequency and code delay quantization errors on the proposed frequency estimators, in the presence of noise.

The Influence of the Peak Point's Location
In this section we study the influence of the quantization error in the frequency domain, which can be re-elaborated as from which it is evident that 0 ≤ Δ p ≤ 0.5. Using Equation (40) and setting f sp as in Equation (33), the estimation error can be expressed aŝ From Equation (53), and using the definition Equation (52), we can calculate the expected value of δ f as From this result we observe that the residual error depends on both the noise contribution in the vector S and the parameter Δ p . From Figure 12 we can see that Δ p influences the accuracy of the proposed three methods, especially when it comes to the top limit 0.5. This result suggested us to adopt a strategy to mitigate this effect. Recalling Figure 4, in which "A" is the peak point, "B" is the second high point and "C" is the lowest point independently from their relative position on the frequency axis, we can easily observe that, when Δ p is close to the limit 0.5, then S C is close to zero. In this case the residual error defined in Equation (39) introduces a significant error in the estimate Equation (40), since noise dominates in point "C". This is experimentally proven in Figure 12, where the curves of R n -3 and C n -3 show an increasing RMSE as Δ p increases, while the LS appears slightly more robust.
So when Δ p is close to 0.5, to limit the accuracy degradation we changed the algorithm (40) as which will be used whenever the point "C" gets close to zero. The idea is to ignore the "C" term, because when Δ p is close to 0.5, the true value S C is nearly zero and S C mainly contains noise. Figure 12. Influence of the quantization error on the RMSE as a function Δ p . To implement this method, a threshold control has to be added, as drawn in Figure 13. Here we use the empirical criterion S B / S A > 0.92 to decide whether Δ p is critically close to 0.5 or not. In Figure 12 we can observe the result of the method drawn in Figure 13 (continuous line with square markers), which is able to consistently reduce the RMSE when Δ p is close to 0.5.

The Influence of the Code Delay Error
The peak column selected in the search space (marked in Figure 4) also depends on the resolution of the search space in the code delay domain. The mentioned quantization error Δτ q = τ p − τ affects the CAF samples with an amplitude scale factor |R(Δτ q )| as shown in Equation (21), where now Δτ = Δτ q . This factor is expected to affect the performance of the estimators C n -3, and R n -3.
The influence of such a code delay error can be quantified by modifying the expression of the signal-to-noise ratio ρ in the CRLB Equation (49), so as to take into account the term R(Δτ q ). Thus, defining the modified SNRρ = A 2 |R(Δτq)| 2 2σ 2 IF , Equation (49) becomes From Equation (56), because R(Δτ q ≤ 1), we can conclude that the CRLB value increases as the code delay Δτ q increases; this increase will be also experienced by the RMSE of the proposed algorithms. In conclusion the effect of the quantization error Δτ q is an attenuation, which does not modify the results of Figures 8 and 9, except for a scaling factor in the abscissa.

Conclusions
In this paper a new family of algorithms is proposed for the fine estimation of the Doppler frequency based on an approximate analytical expression of the CAF. The proposed methods have been analyzed in both ideal (i.e., noiseless) and realistic (i.e., noisy) scenarios and compared with a similar LS approach. The CRLB has been derived and used as benchmark for performance assessment. The influences of non-AWGN nuisances are also analyzed under a theoretical perspective. In application, from the experiments, we can see that the method R n -3 almost achieves the performance of LS, which is very close to CRLB, but the complexity of R n -3 is notably lower. Moreover, performance can be improved by adopting a simple averaging method.