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

Abstract

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.

1. 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 ƒ_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 $\langle {\widehat{f}}_d^a,{\widehat{\tau }}^a\rangle$ of ƒ_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 estimate ${\widehat{f}}_d^a$ 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 ƒ_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.

2. 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.

2.1. 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

$$\text{y}=\left[y\left(0\right)y\left(1\right)\cdot \cdot \cdot y\left(L-1\right)\right]$$

(1)

where L is the total number of the samples, and

$$y\left(n\right)=r\left(n\right)+W\left(n\right)$$

(2)

where r(n) is a signal that contains a vector of unknown parameters a = [α₁α₂α₃ ⋯ α_K], W(n) is a zero-mean White Gaussian Noise (WGN) random process with variance σ², 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

$$\begin{array}{cc}\hfill p\left(\text{y};\overline{\text{a}}\right)& =\frac{1}{{\left(2\pi {\sigma }^2\right)}^{\frac{L}{2}}}exp\left[-\frac{1}{2{\sigma }^2}\sum _{n=0}^{L-1}{\left(y\left(n\right)-\overline{r}\left(n\right)\right)}^2\right]\hfill \\ & =\frac{1}{{\left(2\pi {\sigma }^2\right)}^{\frac{L}{2}}}exp\left[-\frac{1}{2{\sigma }^2}\sum _{n=0}^{L-1}\left(y^2\left(n\right)+{\overline{r}}^2\left(n\right)\right)+\frac{1}{2{\sigma }^2}\sum _{n=0}^{L-1}y\left(n\right)\overline{r}\left(n\right)\right]\hfill \end{array}$$

(3)

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 $\text{∑}_{n=0}^{L-1}{\overline{r}}^2(n)$ in Equation (3)) does not depend on a̅, then it is possible to show [17] that the corresponding ML estimate â_ML of a is

$${\widehat{\text{a}}}_{\mathrm{ml}}=arg\underset{\overline{\text{a}}\in {\mathit{\text{ran}}}_{\text{a}}}{max}\sum _{n=0}^{L-1}y\left(n\right)\overline{r}\left(n\right)$$

(4)

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

$$R\left(\overline{\text{a}}\right)=\sum _{n=0}^{L-1}y(n)\overline{r}(n)$$

(5)

â_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\left(n{T}_s\right)=\sum _{m=1}^{N_{\upsilon }}{y}_m\left(n{T}_s\right)$$

(6)

where N_υ is the number of satellites in view, and

$$y_m(n{T}_s)=A_m{C}_m(n{T}_s-{\tau }_m)d_m(n{T}_s-{\tau }_m)cos(2\pi (f_{IF}+f_{d,m})n{T}_s+{\phi }_m)$$

(7)

where A_m is the amplitude of the signal, C_m(nT_s −τ_m) = c_m (nT_s −τ_m) sb(nT_s −τ_m) is the product of the satellite spreading code c_m(nT_s −τ_m) and subcarrier s_b(nT_s −τ_m) used in the new GNSS systems [18], such as in Galileo (if no subcarrier is present, then s_b(nT_s −τ_m)=1), τ_m is the code delay, d_m(nT_s −τ_m) is the navigation data, f_IF is the intermediate frequency, f_d,m is the Doppler frequency shift, φ_m is the phase of the carrier, and T_s is the sampling interval (the inverse of the sampling frequency f_s).

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

$$\overline{r}(n)=C(n{T}_s-\overline{\tau })e^{2\pi \left(f_{IF}+{\overline{f}}_d\right)n{T}_s}$$

(8)

where the parameter vector becomes $\overline{\text{a}}=\left[\overline{\tau },{\overline{f}}_d\right]$, and the energy of r̅(n) is not related to a̅. So the accumulation process in acquisition can be expressed mathematically as

$$R(\overline{\tau },{\overline{f}}_d)=\frac{1}{L}\sum _{n=0}^{L-1}y(n)C(n{T}_s-\overline{\tau })e^{2\pi \left(f_{IF}+{\overline{f}}_d\right)n{T}_s}$$

(9)

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

$${\widehat{\text{a}}}_{\mathrm{ml}}=arg\underset{{\mathit{\text{ran}}}_{\text{a}}}{max}\left|R(\overline{\tau },{\overline{f}}_d)\right|$$

(10)

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 $\left[-\frac{f_{sp}}{2},\frac{f_{sp}}{2}\right]$, 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.

2.2. 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

$$S_i(\overline{\tau },{\overline{f}}_d)=\sum _{m=1}^{N_{\upsilon }}{R}_{m,i}(\overline{\tau },{\overline{f}}_d)$$

(11)

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]

$$\begin{array}{cc}\hfill R_{m,i}(\overline{\tau },{\overline{f}}_d)& =\frac{1}{T_d}\underset{0}{\overset{T_d}{\int }}{y}_m(t)\cdot \{C_i(t-\overline{\tau })e^{j2\pi {(f_{IF}+{\overline{f}}_d)}^t}\}dt\hfill \\ & =\frac{A_m}{2T_d}{e}^{-j\phi m}\{F\{P_{T_d}(t)\}*F\{b_{m,i}(t)\}\}f=\mathrm{\Delta }f\hfill \\ & =B_m\sum _{k=-\infty }^{\infty }{a}_{k(m,i)}sin\text{c}\left(\left(\mathrm{\Delta }f-\frac{k}{T_p}\right)\pi T_d\right)e^{-j(\mathrm{\Delta }f-\frac{k}{T_p})\pi T_d}\hfill \end{array}$$

(12)

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

$$\mathrm{\Delta }f=f_{d,m}-{\overline{f}}_d$$

(13)

b_m,_i(t) is the product of two spreading codes, that is

$$b_{m,i}\left(t\right)=C_m(t-{\tau }_m)C_i(t-\overline{\tau })$$

(14)

and P_Td is a window function defined as

$$P_{T_d}=\{\begin{array}{cc}1,& 0<t<T_d\\ 0,& \text{otherwise}\end{array}$$

(15)

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

$$a_{k(m,i)}=\frac{1}{T_p}{\int }_{T_p}{b}_{m,i}(t)e^{-jk(2\pi /T_p)}dt$$

(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 −r̅ = 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

$$\sum _{k=-\infty }^{\infty }{a}_{k(i,i)}\left|\text{Sinc}\left(\left(\mathrm{\Delta }f-\frac{k}{T_p}\right)\pi T_d\right)\right|$$

(17)

to demonstrate the relationship among different Sinc functions. As Figure 3 shows, if a₀ is the peak component in Equation (12) (the subscript i,i is omitted to simplify the notations) only the components strictly adjacent to a₀(i.e., a_-2a_-1a₁a₂ …) 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_-2a_-1a₁a₂ …) are far smaller than α₀, the mathematical expression of the CAF in the main lobe (subscript “ml”) can be written as

$$\begin{array}{cc}\hfill S_i{(\overline{\tau },{\overline{f}}_d)}_{\text{ml}}& =R_{i,i}{(\overline{\tau },{\overline{f}}_d)}_{\mathrm{ml}}+\sum _{m\ne i}{R}_{m,i}{(\overline{\tau },{\overline{f}}_d)}_{\mathrm{ml}}\hfill \\ & \approx \frac{A_i}{2}{e}^{-{\phi }_i}{a}_0\text{sinc}\left(\mathrm{\Delta }f\pi {\text{T}}_d\right)e^{-\text{j}\mathrm{\Delta }f\pi T_d}\hfill \end{array}$$

(18)

where α₀ = α_0(i,i) = R(Δτ), and Δτ = τ_m − τ̅. As a conclusion the approximate expression of the CAF in its main lobe is

$$S_i{(\overline{\tau },{\overline{f}}_d)}_{\text{ml}}\approx \frac{A_i}{2}{e}^{-j{\phi }_i}R(\mathrm{\Delta }\tau )\text{Sinc}(\mathrm{\Delta }f\pi T_d)e^{-j\mathrm{\Delta }f\pi T_d}$$

(19)

The validity of this approximation can be improved in two ways:

• 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.

• By decreasing the values of the adjacent coefficients (a₋₂a₋₁a₁a₂ …). 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.

3. 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 f̅_d ∈ (−f_dmax, f_dmax) (Y-axis range). 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

$$\overline{S}{(\overline{\tau },{\overline{f}}_d)}_{\text{ml}}=\left|S_i{(\overline{\tau },{\overline{f}}_d)}_{\text{ml}}\right|$$

(20)

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.

The cells A, B and C are characterized by the triplete and (τ̅_A,f̅_AS̅_a),(τ̅_B,f̅_BS̅_B),and(τ̅_C,f̅_CS̅_C), 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̅_dX) _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.

3.1. 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]

$$\begin{array}{cc}{\overline{S}}_A& =\frac{A}{2}\left|R(\mathrm{\Delta }\tau )\right|\left|\text{Sinc}\left(\pi {\text{LT}}_{\text{s}}\right({\text{f}}_{\text{d}}-{\overline{f}}_A\left)\right)\right|\\ {\overline{S}}_B& =\frac{A}{2}\left|R(\mathrm{\Delta }\tau )\right|\left|\text{Sinc}\left(\pi L{T}_s\right(f_d-{\overline{f}}_B\left)\right)\right|\\ {\overline{S}}_C& =\frac{A}{2}\left|R(\mathrm{\Delta }\tau )\right|\left|\text{Sinc}\left(\pi L{T}_s\right(f_d-{\overline{f}}_C\left)\right)\right|\end{array}$$

(21)

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) we write f̅ = f_d − Δf, so

$$\begin{array}{ll}\overline{f}\overline{S}\hfill & =(f_d-\mathrm{\Delta }f)\overline{S}\hfill \\ & =f_d\overline{S}-\frac{A\left|R(\mathrm{\Delta }\tau )\right|sin(\pi L{T}_s\mathrm{\Delta }f)}{2\pi L{T}_s}\hfill \end{array}$$

(22)

Based on Equations (21) and (22), we can write

$${\overline{f}}_A{\overline{S}}_A+{\overline{f}}_B{\overline{S}}_B+{\overline{f}}_C{\overline{S}}_C={\overline{f}}_d{\overline{S}}_A+{\overline{f}}_d{\overline{S}}_B+{\overline{f}}_d{\overline{S}}_C-S_{A,B,C}$$

(23)

where

$$S_{A,B,C}=\frac{A\left|R(\mathrm{\Delta }\tau )\right|}{2\pi L{T}_s}(sin(\pi L{T}_s\mathrm{\Delta }{f}_A)+(sin(\pi L{T}_s\mathrm{\Delta }{f}_B)+(sin(\pi L{T}_s\mathrm{\Delta }{f}_C))$$

(24)

Considering f_sp = 2/(3LT_S):

$$\begin{array}{lll}{\overline{f}}_B\hfill & ={\overline{f}}_A+\frac{2}{3L{T}_s}\hfill & \Rightarrow \mathrm{\Delta }{f}_B=\mathrm{\Delta }{f}_A+\frac{2}{3L{T}_s}\hfill \\ {\overline{f}}_C\hfill & ={\overline{f}}_A+\frac{2}{3L{T}_s}\hfill & \Rightarrow \mathrm{\Delta }{f}_C=\mathrm{\Delta }{f}_A+\frac{2}{3L{T}_s}\hfill \end{array}$$

(25)

from which

$$S_{A,B,C}=sin\left(\pi L{T}_s\mathrm{\Delta }{f}_A+\frac{2\pi }{3}\right)+sin\left(\pi L{T}_s\mathrm{\Delta }{f}_A-\frac{2\pi }{3}\right)=0$$

(26)

By substituting Equation (26) for S_A,B,C in Equation (23) we obtain f̅_AS̅_A + f̅_BS̅_B + f̅_CS̅_C = f_d(S̅_A + S̅_B + S̅_C), from which we can write the true Doppler frequency as

$$f_d=\frac{{\overline{f}}_A{\overline{S}}_A+{\overline{f}}_B{\overline{S}}_B+{\overline{f}}_C{\overline{S}}_C}{{\overline{S}}_A+{\overline{S}}_B+{\overline{S}}_C}$$

(27)

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

$$f_{sp}=\frac{2}{\overline{n}L{T}_s}$$

(28)

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

$$\sum _{k=0}^{\overline{n}-1}sin\left(\pi L{T}_s\mathrm{\Delta }{f}_1+\frac{2k}{\overline{L}{T}_s}\right)=0$$

(29)

where Δf1=f_d−f̅₁ and f̅₁ 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

$$f_d=\frac{\text{∑}_{k=1}^{\overline{n}}{\overline{f}}_k{\overline{S}}_k}{\sum _{k=1}^{\overline{n}}{\overline{S}}_k}$$

(30)

which can be adopted as an estimator of the Doppler frequency in the presence of noise.

For simplicity, hereafter we refer to solution Equation (27) as algorithm R-3 and solution Equation (30) as algorithm R-n̅.

3.2. 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̅(r̅,f̅_d)_ml in the complex form given in Equation (19). In fact it is possible to show that

$$\begin{array}{ll}{f}_{d,\mathit{\text{complex}}}\hfill & =\text{Real}\left\{\frac{{\overline{f}}_C{S}_C+{\overline{f}}_A{S}_A{e}^{-j\pi f_{sp}{T}_d}+{\overline{f}}_B{S}_B{e}^{-j\pi f_{sp}{T}_d}}{S_C+S_A{e}^{-j\pi f_{sp}{T}_d}+S_B{e}^{-j\pi f_{sp}{T}_d}}\right\}\hfill \\ & =\text{Real}\left\{\frac{\left({\overline{f}}_C\left|S_C\right|+{\overline{f}}_A\left|S_A\right|+{\overline{f}}_B\left|S_B\right|\right)e^{-j(\pi \mathrm{\Delta }{f}_C{T}_d+{\phi }_i)}}{(\left|S_C\right|+\left|S_A\right|+\left|S_B\right|)e^{-j(\pi \mathrm{\Delta }{f}_C{T}_d+{\phi }_i)}}\right\}\hfill \\ & =\frac{{\overline{f}}_C\left|S_C\right|+{\overline{f}}_A\left|S_A\right|+{\overline{f}}_B\left|S_B\right|}{\left|S_C\right|+\left|S_A\right|+\left|S_B\right|}\hfill \\ & =f_d\hfill \end{array}$$

(31)

At the same time, similarly to Equation (30), we can generalize this solution as

$$f_{d,\mathit{\text{complex}}}=\text{Real}\left\{\frac{\text{∑}_{k=1}^{\overline{n}}{\overline{f}}_k{S}_k{e}^{-j(k-1)\pi f_{sp}{T}_d}}{\text{∑}_{k=1}^{\overline{n}}{S}_k{e}^{-j(k-1)\pi f_{sp}{T}_d}}\right\}$$

(32)

where f̅_k+1=f̅_k+f_sp and ƒ₁ 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̅.

3.3. 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

$$f_{sp}=\frac{2}{3L{T}_s}=167{\text{H}}_{\text{Z}}$$

(33)

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.

4. 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

$$\begin{array}{ll}\stackrel{\sim }{y}(n{T}_s)\hfill & =y(n{T}_s)+\eta (n{T}_s)\hfill \\ & =Ac(n{T}_s-\tau )s_b(n{T}_s-\tau )d(n{T}_s-\tau )\hfill \\ & cos(2\pi (f_{IF}+f_d)n{T}_s+\phi )+\eta (n{T}_s)\hfill \end{array}$$

(34)

where η(nT_s) represents the white Gaussian noise normally distributed with zero mean and variance ${\sigma }_{IF}^2$, related to the power spectral density S_N(f) = N₀/2 of the analogue noise by the well-known formula

$${\sigma }_{IF}^2=E\{{\eta }^2(n{T}_s)\}=\frac{N_0{f}_s}{2}$$

(35)

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

$$\begin{array}{ll}I(\overline{\tau },{\overline{f}}_d)\hfill & =\frac{1}{N}\sum _{n=0}^{N-1}(y(n{T}_s)+\eta (n{T}_s))C(n{T}_s-\overline{\tau })cos(2\pi (f_{IF}+{\overline{f}}_d))\hfill \\ Q(\overline{\tau },{\overline{f}}_d)\hfill & =\frac{1}{N}\sum _{n=0}^{N-1}(y(n{T}_s)+\eta (n{T}_s))C(n{T}_s-\overline{\tau })sin(2\pi (f_{IF}+{\overline{f}}_d))\hfill \end{array}$$

(36)

which, based on Equation (19) (and omitting the subscript i), becomes

$$\begin{array}{ll}I(\overline{\tau }{\overline{f}}_d)\hfill & =\text{Real}\left\{S(\overline{\tau },{\overline{f}}_d)\right\}+N_I\hfill \\ Q(\overline{\tau }{\overline{f}}_d)\hfill & =\text{Img}\left\{S\left\{S(\overline{\tau },{\overline{f}}_d)\right\}+N_Q\right\}\hfill \end{array}$$

(37)

where Real{ } and Img{ } mean respectively real and imaginary part of a complex value, and

$$\begin{array}{ll}{N}_I\hfill & =\frac{1}{N}\sum _{n=0}^{N-1}\eta (n{T}_s)C(n{T}_s-\overline{\tau })cos(2\pi (f_{IF}+{\overline{f}}_d))\hfill \\ N_Q\hfill & =\frac{1}{N}\sum _{n=0}^{N-1}\eta (n{T}_s)C(n{T}_s-\overline{\tau })sin(2\pi (f_{IF}+{\overline{f}}_d))\hfill \end{array}$$

(38)

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)={\sigma }_{IF}^2/(2N)$, and the envelope of S̃(τ̅, f̅d) = I + jQ is

$$\begin{array}{ll}\stackrel{\sim }{\overline{S}}{(\overline{\tau },{\overline{f}}_d)}^2\hfill & =I{(\overline{\tau },{\overline{f}}_d)}^2+Q{(\overline{\tau },{\overline{f}}_d)}^2\hfill \\ & =\overline{S}{(\overline{\tau },{\overline{f}}_d)}^2+2\text{Real}\{S(\overline{\tau },{\overline{f}}_d)\}{N}_I+2\text{Img}\{S(\overline{\tau },{\overline{f}}_d)\}{N}_Q+N_I^2+N_Q^2\hfill \\ & =\overline{S}{(\overline{\tau },{\overline{f}}_d)}^2+{\omega }_I+{\omega }_Q\hfill \end{array}$$

(39)

where the term w_I+w_q = w is a random process with mean $E\{w\}=E\{N_1^2+N_Q^2\}={\sigma }_{IF}^2/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

$${\widehat{f}}_d=\frac{{\overline{f}}_A{\stackrel{\sim }{\overline{S}}}_A+{\overline{f}}_B{\stackrel{\sim }{\overline{S}}}_B+{\overline{f}}_C{\stackrel{\sim }{\overline{S}}}_C}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B+{\stackrel{\sim }{\overline{S}}}_C}$$

(40)

and

$${\widehat{f}}_{d,\mathit{\text{complex}}}=\text{Real}\left\{\frac{{\overline{f}}_C{\stackrel{\sim }{S}}_C+{\overline{f}}_A{\stackrel{\sim }{S}}_A{e}^{-j\pi f_{sp}{T}_d}+{\overline{f}}_B{\stackrel{\sim }{S}}_B{e}^{-j\pi f_{sp}{T}_d}}{{\stackrel{\sim }{S}}_C+{\stackrel{\sim }{S}}_A{e}^{-j\pi f_{sp}{T}_d}+{\stackrel{\sim }{S}}_B{e}^{-j\pi f_{sp}{T}_d}}\right\}$$

(41)

(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.

4.1. Least-Square Method

Another approach to exploit the CAF points ${\stackrel{\sim }{\overline{S}}}_A,{\stackrel{\sim }{\overline{S}}}_B,{\stackrel{\sim }{\overline{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

$$Q_{LS}(\overline{\text{a}},f_d)={\left({\stackrel{\sim }{\overline{S}}}_A-{\overline{S}}_A\right)}^2+{\left({\stackrel{\sim }{\overline{S}}}_B-{\overline{S}}_B\right)}^2+{\left({\stackrel{\sim }{\overline{S}}}_C-{\overline{S}}_C\right)}^2$$

(42)

is minimized with respect to the unknown parameters, where A̅ = A/(2R(Δ_τ)). This leads to the equations

$$\begin{array}{c}\frac{\partial Q_{LS}(\overline{A},f_d)}{\partial \overline{\text{a}}}=0\\ \frac{\partial Q_{LS}(\overline{\text{a}},f_d)}{\partial f_d}=0\end{array}$$

(43)

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.

4.2. 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

$$\begin{array}{ll}\stackrel{\sim }{y}(n{T}_S)\hfill & =y(n{T}_S)+\eta (n{T}_S)\hfill \\ & =Acos(2\pi (f_{IF}+f_d)n{T}_s+\phi )+\eta (n{T}_S)\hfill \end{array}$$

(44)

which contains the unknown parameter vector θ = [A, f_d, φ]. The corresponding Fisher information matrix I for an observation of N samples has elements [22]

$${[\mathbf{I}(\mathit{\theta })]}_{ij}=\frac{1}{{\sigma }^2}\sum _{n=0}^{N-1}\frac{\partial y(n;\mathit{\theta })}{\partial {\mathit{\theta }}_i}\frac{\partial y(n;\mathit{\theta })}{\partial {\mathit{\theta }}_j}$$

(45)

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

$$\mathbf{I}(\mathit{\theta })=\frac{1}{{\sigma }_{IF}^2}\left[\begin{array}{ccc}\frac{N}{2}& 0& 0\\ 0& 2{\pi }^2{A}^2{T}_s^2\sum _{n=0}^{N-1}{n}^2& \pi A^2{T}_s\sum _{n=0}^{N-1}n\\ 0& \pi A^2{T}_s\sum _{n=0}^{N-1}n& \frac{N{A}^2}{2}\end{array}\right]$$

(46)

The CRLB is found as the [i, i] element of the inverse of the Fisher matrix:

$$var({\mathit{\theta }}_i)\ge {[{\mathbf{\text{I}}}^{-1}(\mathit{\theta })]}_{ii}$$

(47)

Therefore the CRLB for the Doppler frequency estimate is

$$\begin{array}{ll}var({\widehat{f}}_d)\hfill & \ge {[{\mathbf{\text{I}}}^{-1}(\mathit{\theta })]}_{22}\hfill \\ & \ge \frac{12}{{(2\pi )}^2\rho T_s^{2}N(N^2-1)}\hfill \end{array}$$

(48)

where $\rho =A^2/\left(2{\sigma }_{IF}^2\right).$ As N = f_sT_d ≫ 1, Equation (48) can be written as

$$\begin{array}{ll}var({\widehat{f}}_d)\hfill & \ge {[{\mathbf{\text{I}}}^{-1}(\mathit{\theta })]}_{22}\hfill \\ & \ge \frac{12}{{(2\pi )}^2\rho f_s{T}_d^3}\hfill \end{array}$$

(49)

which is the benchmark for our estimation algorithms.

4.3. 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:

$$\mathit{\text{SNR}}=\frac{C}{N_{0}B}$$

(50)

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:

$$f_{\mathit{\text{RMSE}}}=\sqrt{E\left\{{({\widehat{f}}_d-f_d)}^2\right\}}$$

(51)

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_rmse 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 $\sqrt{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.

4.4. 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:

• 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)

• Use the parameter 〈f̅_d, τ̅〉 evaluated in the first step to calculate the new columns (as shown in Figure 10) of the next (M − 1) periods of data.

• Calculate the average of the M columns into one single mean column.

• Pick out the top n̅ cells in the mean column and use Equation (40) to calculate the final Doppler frequency.

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.

5. 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 f = [f̅_C, f̅_A, f̅_B] and S̅ = [S_C, S_A, S_B] used in Equations (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 vector S̅, 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.

5.1. 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

$${\mathrm{\Delta }}_p=\left|\frac{\mathrm{\Delta }fq}{f_{sp}}\right|$$

(52)

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 as

$$\begin{array}{ll}{\widehat{f}}_d-f_d={\delta }_f\hfill & =\frac{{\overline{f}}_A{\stackrel{\sim }{\overline{S}}}_A+{\overline{f}}_B{\stackrel{\sim }{\overline{S}}}_B+{\overline{f}}_C{\stackrel{\sim }{\overline{S}}}_C}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B+{\stackrel{\sim }{\overline{S}}}_C}-f_d\hfill \\ & =\frac{(f_d+\mathrm{\Delta }{f}_A){\stackrel{\sim }{\overline{S}}}_A+(f_d+\mathrm{\Delta }{f}_A+f_{sp}){\stackrel{\sim }{\overline{S}}}_B+(f_d+\mathrm{\Delta }{f}_A-f_{sp}){\stackrel{\sim }{\overline{S}}}_C}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B+{\stackrel{\sim }{\overline{S}}}_C}-f_d\hfill \\ & =\frac{\mathrm{\Delta }{f}_A{\stackrel{\sim }{\overline{S}}}_A+(\mathrm{\Delta }{f}_A+f_{sp}){\stackrel{\sim }{\overline{S}}}_B+(\mathrm{\Delta }{f}_A-f_{sp}){\stackrel{\sim }{\overline{S}}}_C}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B+{\stackrel{\sim }{\overline{S}}}_C}\hfill \\ & =\mathrm{\Delta }{f}_A+\frac{f_{sp}({\stackrel{\sim }{\overline{S}}}_B|-{\stackrel{\sim }{\overline{S}}}_C)}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B+{\stackrel{\sim }{\overline{S}}}_C}\hfill \end{array}$$

(53)

From Equation (53), and using the definition Equation (52), we can calculate the expected value of δ_f as

$$E\{{\delta }_f\}=f_{sp}E\left\{{\mathrm{\Delta }}_p+\frac{{\stackrel{\sim }{\overline{S}}}_B-{\stackrel{\sim }{\overline{S}}}_C}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B+{\stackrel{\sim }{\overline{S}}}_C}\right\}$$

(54)

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

$${\widehat{f}}_{d,\mathit{\text{rel}}}=\frac{{\overline{f}}_A{\stackrel{\sim }{\overline{S}}}_A+{\overline{f}}_B{\stackrel{\sim }{\overline{S}}}_B}{{\stackrel{\sim }{\overline{S}}}_A+{\stackrel{\sim }{\overline{S}}}_B}$$

(55)

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̅_Cis nearly zero and ${\stackrel{\sim }{\overline{S}}}_C$ mainly contains noise.

To implement this method, a threshold control has to be added, as drawn in Figure 13. Here we use the empirical criterion ${\stackrel{\sim }{\overline{S}}}_B/{\stackrel{\sim }{\overline{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.

5.2. 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 $\mathit{\text{SNR}}\stackrel{\sim }{\rho }={\frac{A^2\left|R(\mathrm{\Delta }{\tau }_q)\right|}{2{\sigma }_{IF}^2}}^2,$ Equation (49) becomes

$$var(\widehat{f_d})\ge \frac{12}{{(2\pi )}^2\stackrel{\sim }{\rho }{f}_s{T}_d^3}$$

(56)

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.

6. 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.