GPS Pull-In Search Using Reverse Directional Finite Rate of Innovation (FRI)

Abstract

When an incoming Global Positioning System (GPS) signal is acquired, pull-in search performs a finer search of the Doppler frequency of the incoming signal so that phase lock loop can be quickly stabilized and the receiver can produce an accurate pseudo-range measurement. However, increasing the accuracy of the Doppler frequency estimation often involves a higher computational cost for weaker GPS signals, which delays the position fix. In this paper, we show that the Doppler frequency detectable by a long coherent auto-correlation can be accurately estimated using a complex-weighted sum of consecutive short coherent auto-correlation outputs with a different Doppler frequency hypothesis, and by exploiting this we propose a noise resistant, low-cost and highly accurate Doppler frequency and phase estimation technique based on a reverse directional application of the finite rate of innovation (FRI) technique. We provide a performance and computational complexity analysis to show the feasibility of the proposed technique and compare the performance to conventional techniques using numerous Monte Carlo simulations.

1. INTRODUCTION

Successful positioning in harsh Global Positioning System (GPS) environments has been one of the most challenging problems in GPS, and the demand for location-based service (LBS) in harsh environments for GPS, such as dense urban and indoor environments is continuously increasing. To meet the needs for GPS in harsh environments, GPS receivers may need a faster and higher sensitivity acquisition and tracking functions that can produce accurate carrier phase and code phase estimates in various GPS signal conditions. In the literature, a number of acquisition techniques for weak GPS signals have been introduced. For example, a GPS receiver connected to a synchronized cellular communication network (Kransner 1998, van Diggelen 2009) can make use of the downlink signal measurements and downlink information from the cellular base station to quickly acquire GPS signals with high sensitivity. On the other hand, conventional techniques to enhance the acquisition sensitivity and to reduce the mean acquisition time (MAT) such as long coherent integration and non-coherent accumulation techniques (Parkinson et al. 1996) have been continuously studied, and, recently, techniques such as parallel signal search, FFT-based signal search (Borre et al. 2007), and two-dimensional compressed correlator (Kong & Kim 2013) techniques have been developed.

To produce accurate measurements of code phase and carrier phase of incoming GPS signals, GPS receivers employ tracking functions that require initial code phase and Doppler frequency estimates with sufficient accuracy. In practice, the accuracy required for the initial resolutions of the code phase and Doppler frequency are less than a half chip and about ten Hertz (Tsui 2005), respectively. Since most of the acquisition techniques perform Doppler frequency search for incoming GPS signals with a search step size about \(1/(2T)\) Hz, where the coherent integration time \(T\) is usually 1~5 ms, the Doppler frequency estimate may have an error of up to 250 Hz, which can be too coarse for a carrier phase tracking loop. However, an initial code phase error of up to a half chip can be tolerable for most of the code phase tracking functions (Irsigler & Eissfeller 2003). Therefore, a GPS receiver has a pull-in state between the acquisition and tracking states to increase the accuracy of the Doppler frequency estimate to a high enough level to accurately run carrier phase tracking. In the pull-in state, a GPS receiver tries to obtain a fine Doppler frequency, carrier phase, and \(C/N_0\) estimates of an incoming GPS signal, which can take about 1.5s (Kelley et al. 2002). To improve the pull-in search performance, conventional techniques, such as Tsui (2005) and Goiser & Berger (1996), perform a fast Fourier transform (FFT) to find the maximum frequency component from code wiped-off sample data and then utilize the amplitude of the component to obtain a fine Doppler frequency estimate. Consecutive fine Doppler frequency estimates can be averaged to increase accuracy in the presence of noise, but, in general, the conventional techniques work well for high \(C/N_0\). To lessen the noise effect on the Doppler frequency estimation, Sagiraju et al. (2006) suggest increasing \(T\) up to 5 ms, however, the scheme requires a smaller initial Doppler frequency estimation error than 50 Hz, which, in effect increases the complexity of the acquisition state. On the other hand, Zeng & Li (2010) suggest a serial Doppler frequency search using correlation in the time domain, and employ a curve fitting algorithm to find a more accurate Doppler frequency estimate. The curve fitting may be able to provide good Doppler frequency and carrier phase estimates at the cost of increased computational complexity. In general, since a GPS receiver does not yet have the knowledge on the received data bit in the pull-in state, increasing the accuracy of Doppler frequency and carrier phase estimates of weak GPS signals can be a complex problem due to the computational complexity and algorithmic complexity required to mitigate the effect of bit transition. In practice, the carrier phase estimation is performed by the phase locked loop in the tracking state.

Despite the complexity and importance, there has been little attention paid to the pull-in search of GPS receivers in the literature. To reduce the computational complexity and improve the accuracy of the Doppler frequency and carrier phase estimates in the pull-in state, we propose a reverse-directional finite rate of innovation (FRI) technique with iterative Cadzow denoising applied to consecutive coherent auto-correlation function (ACF) outputs to reduce the computational complexity in estimating the Doppler frequency and carrier phase of an incoming GPS signal with high accuracy and to increase noise robustness. The complexity of the proposed technique is analyzed and the performance of the proposed technique is compared to the conventional techniques using numerous Monte Carlo simulations.

The rest of this paper is organized as follows. In Section 2, we introduce a general continuous form expression for a coherent ACF output, and in Section 3, we verify that a long coherent ACF output for a Doppler frequency can be estimated with consecutive short coherent ACF outputs at a close Doppler frequency. Section 4 provides the details and complexity of the proposed technique that exploits the findings in Section 3, and Section 5 shows the analysis of the performance of the proposed technique and a performance comparison to conventional techniques using numerous Monte Carlo simulations. And a conclusion is drawn in Section 6.

Throughout this paper, the following conventions are used for notation. Vectors or matrices are denoted by boldface symbols. Small letters are used for scalars and vectors, and capital letters are used for matrices.

2. COHERENT ACF OUTPUT OF GPS SIGNAL

Let \(r(t)\) be the frequency down-converted incoming \(L1\) GPS coarse acquisition (C/A) code signal to an IF frequency \(f_{IF}\) , then the complex IF signal can be expressed as

\(r(t)=\sum\limits_{q=1}^Q A_qD_q(t-\tau_q)C_q(t-\tau_q)e^{j[2\pi(fIF+fDq)t+\phi q]}+n(t),  \)                                                                             (1)

where \(Q\) is the number of GPS satellite signals being received, \(A_q, \tau_q, |f_D^q|(\le 5~{\rm kHz}),\) and \(\phi_q\) are the slowly varying amplitude, code phase, Doppler frequency, and unknown carrier phase of the \(q\)-th satellite signal, respectively, \(D_q(t)\) and \(C_q(t)\) are the navigation data bit signal and Pseudo-Random Noise code signal of the q-th satellite at time t, respectively, and \(n(t)\) represents an additive whiteGaussian noise with two-sided power spectral density (PSD) \(N_0/2\). Note that the data bit and chip rates are \(R_b\)(=50 bps) and \(R_c\)(= 1.023 Mcps), respectively, and thet the C/A code has code length \(L_c=1023\) chips so that the code period is only \(T_1\)(= 1 ms). When the coherent integration (correlation) length used in the acquisition state is \(T=L_T T_1\), where \(L_T\) is an integer much larger than 1, and code phase search resolution is \(T_c/2\), the detected code phase \(\tau_q\) and Doppler frequency at the acquisition state can have errors upto \(T_c/4\) and \(1/(4T)\), respectively. Note that the navigation data remains unknown in the pull-in state. To detect and track the code phase \(\tau_q\) and Doppler frequency \(f_D^q\) of the \(q\)-th satellite signal using correlation, the receiver generates a prompt replica signal.

\(r_q(t)=C(t-\hat{\tau}_q)e^{j2\pi(fIF+{\hat{f}_D^q)}t},\)                                                                                                 (2)

where \(\hat{\tau}_q\) and \({\hat{f}_D^q}\) are the code phase and carrier frequency estimates made by the receiver’s acquisition function, respectively. The receiver integrates the product of \(r(t)\) and \(r_q^*(t)\) (complex conjugate of \(r_q(t)\)) over the coherent integration interval \(T\) to see the correlation between the two functions. Defining \(\delta\tau=\tau_q-\tau_q\) and \(\delta f = f_D^q - \hat{f}_D^q\), the mathematical expression of the ACF output over the two-dimensional space is readily available from Kong & Kim (2013).

\(R_2 (\delta\tau,\delta f)={1 \over T} \displaystyle\int_{t_0}^{t_0+T} r(t) r_q^*(t)dt\)                                                                                                        

\(= {A_qD_q \over T} \displaystyle\int_{t_0}^{t_0+T} C(t-\tau_q) C(t-\hat{\tau}_q) e^{j(2\pi \delta ft+\phi)} dt + w\)                                                               

\(= {A_qD_q \sin (\pi\delta fT) \over \pi\delta fT} \times {\sin(\pi\delta f(T_c-|\delta\tau|)) \over \sin(\pi\delta fT_c)} \times e^{j(2\pi\delta ft +\pi\delta ft - \pi\delta f |\pi\delta | + \phi) }+ w\)                                            (3)

where \(D_q(t - \tau_q)\) is assumed as constant \(D_q\) for the integration interval \(T\), and

\(w = {1 \over T} \displaystyle\int_{t_0}^{t_0+T} \left [ \sum\limits_{p=1,p≠q}^Q A_p D_p (t-\tau_p) C_p(t-\tau_p) \times e^{j[2\pi(f_{IF}+ f_D^p ) t + \phi_p ]} + n(t) \right ] r_q(t)^* dt\)                                                                    

\(\simeq {1 \over T} \displaystyle\int_{t_0}^{t_0+T_1} n(t) r_q (t)^* dt\)                                                                                                                                                   (4)

is a complex Gaussian noise with two-sided power spectral density \(N_0/2\), since signal powers are negligible when compared to the noise power. When \(T=T_b(=1/R_b)\) and when \(\phi, \delta\tau\), and \(\delta f\) are small enough, \(R_2(\delta \tau, \delta f)≥0\) indicates that \(D_q=+1\), and \(R_2(\delta \tau, \delta f)<0\) indicates that \(D_q=-1\). Note that there are several techniques to remove the effect of the navigation data signal \(D_q(t)\) from \(R_2(\delta\tau, \delta f)\); a number of bit sequences are estimated using Kalman filter (Psiaki & Jung 2002), the binary value of \(D_q(t)\) is estimated by testing different bit combinations (Petovello et al. 2008), or with the assistance information from the cellular base station in the Assisted-GPS system (Kransner 1998, van Diggelen 2009). Since figuring out the data \(D_q(t)\) is one of the primary tasks in the tracking state, we neglect the effect of the unknown data in this paper.

3. FREQUENCY ESTIMATION WITH CONSECUTIVE SHORT COHERENT ACF OUTPUTS

In this section, we break \(R_2(\delta\tau, \delta f)\) obtained with a long coherent interval \(T \gg T_1\) into \(LT(=T/T_1)\) consecutive \(R_2(\delta\tau, \delta f)\)s obtained with short coherent interval \(T_1\) and investigate the possibility of precisely estimating the Doppler frequency of \(r(t)\) with the \(L_T\) \(R_2(\delta\tau, \delta f)\)s obtained. In the following analysis, we drop the subscript \((·)_q\) (unless necessary) for simplicity in algebraic expressions.

Let the input to the pull-in state be denoted by \(R_2(\delta\tau,\delta f) = R_2^l (\delta\tau,\delta f)|l=0,1,…,L_T-1,\) where

\(R_2^l(\delta\tau,\delta f) = {A \over T_1} \displaystyle\int_{lT_1}^{(l+1)T_1} \left[ C(t-\tau) e^{j(2\pi f_D^t+\phi_q)}+n(t) \right ] \times C(t-\hat{\tau})e^{-j2\pi {\hat{f}_D^t}} dt\)                                                                   

\(= {A\sin\pi\delta fT_1) \over \pi\delta fT_1} \cdot {\sin(\pi\delta f(T_c-|\delta\tau|)) \over \sin(\pi\delta fT_c)} \times e^{j(\pi\delta f[2t_l+T_1-|\delta\tau|]+\phi)}+n_1^l,\)                                                                             (5)

is alreadly available from Kong & Kim (2013) and \((·)^l\) represents that the coherent integration from \(t=lT_1=t_l~ {\rm to}~ t=(l+1)T_1=t_l+T_1\), and

\(n_1^l = {1 \over T_1} \displaystyle\int_{t_1}^{t_l+T_1} n(t)r_q^*(t) dt\)                                                                           (6)

is a complex Gaussian process with two-sided PSD \(N_0/2\). Note that in Eq. (5) the expected ACF output for |\(\delta\tau| \geq T_c\) is assumed as zero. When the code phase estimation is successful, i.e., \(\delta\tau=T_c/2\), the expression Eq. (5) can be simplified to a single-dimensional function

\(R_2^l(\delta\tau,\delta f)|_{\delta\tau = 0} = {A\sin(\pi\delta fT_1) \over \pi\delta fT_1} e^{j[\pi(2l+1)\delta fT_1+\phi]}+n_1^l\)                                                                          

\(= R_1^l(\delta f)\)                                                                           (7)

and let

\(R_1(\delta f)= \sum\limits_{l=0}^{L_T-1} {R_1^l (\delta f) \over L_T}\)                                                                          

\(= \sum\limits_{l=0}^{L_T-1} \left [ {A\sin(\pi\delta fT_1) \over \pi\delta fL_T T_1} e^{j[π(2l+1)δfT_1+ϕ]} + {n_1^l \over L_T} \right ].\)                                                                           (8)

when the carrier frequency of \(r_q(t)\) is \({\hat{f}}_D+ \Delta f\), it is found from Eq. (7) that

\(R_1^l(\delta f-\Delta_f) = {A\sin[\pi(\delta f-\Delta f)T_1] \over \pi(\delta f-\Delta_f)T_1} \times e^{j[(2l+1)\pi(\delta f-\Delta_f)T_1+\phi]}+w^l(\Delta_f),\)                                                                           (9)

where

\(w^l(\Delta_f) = \displaystyle\int_{lT_1}^{(l+1)T_1} n(t)C(t-\hat{\tau})e^{-j2\pi({\hat{f}_D}+\Delta_f)t} dt\)                                                                                                

\(\simeq e^{-j(2l+1)\pi\Delta_f T_1} n_1^l\)                                                                                                                           (10)

which is a complex Gaussian process with the same PSD to \(n_1^l\) Eq. (6). Note that Eq. (10) holds only for a very small \(|\Delta f| T_1 = 1\). From Eqs. (7, 9, 10), it is found that

\(R_1(\delta f-\Delta_f)= \sum\limits_{l=0}^{L_T-1} { R_1^l(\delta f-\Delta_f) \over L_T}\)                                                                                                   (11a)

\(\simeq {\beta(\delta f,\Delta_f) \over L_T} \sum\limits_{l=0}^{L_T-1} R_1^l(\delta f)e^{-j(2l+1)\pi\Delta_f T_1},\)                                                                     (11b)

which indicates that a long coherent ACF output \(R_1(\delta f - \Delta_f)\) with a long coherent integration interval \(T\) and a Doppler frequency estimate \(\hat{f}_D + \Delta_f\) can be obtained from the linear sum of \(L_T\) consecutive short coherent ACF outputs \(R_1^l(\delta f)\) from \(L_T\) integration segments \([lT_1, (l+1)T_1]~(l \in \{0,1,…,R_T-1\})\) of \([0, T_1]\) with a neighboring Doppler frequency \(\hat{f}_D\), when \(|\Delta_f| = 1/T_1\) and \(|\delta f+\Delta_f| = 1/T_1\) are satisfied. In other words, Eq. (11b) states that it is possible to precisely estimate a long coherent ACF output \(R_1(\delta f - \Delta_f)\) for a Doppler frequency hypothesis using consecutive short coherent ACF outputs \(R_1^l(\delta f)\) \((l \in \{0,1,…,R_T-1\})\) with a close Doppler frequency hypothesis. This observation is exploited in the proposed technique for an accurate and fast GPS pull-in search. In addition, it should be noted that the carrier frequency \(\delta f\) of \(\mathbf{R}_1(\delta f)=[R_1^0(δf),R_1^1(δf),…,R_1^{L_T-1}(\delta f)]\) can be estimated by compensating the frequency \(\delta f\) with \(\Delta_f\) which maximizes the linear sum in Eq. (11b).

Fig. 1 shows 104 Monte Carlo simulation results demonstrating the estimation accuracy of the expression Eq. (11b) for a random Doppler frequency \(\delta f\) of the incoming signal from 0Hz to 250 Hz, when \(C/N_0\)=40 dB-Hz and \(L_T=25\). Fig. 1a shows that the result of absolute normalized amplitude difference, i.e., \(\left | \left [ R_1(0)-L_T^{-1}\sum\limits_{l=0}^{L_T-1} R_1^l(δf)e^{-j(2l+1)πδfT_1} /R1(0) \right ] \right |\) for \(\delta f\) varying from 0 to 250 Hz. Fig. 1b shows that the standard deviation of the Doppler frequency and the estimated Doppler frequency using Eq. (11b) is small, also. As a result, it is demonstrated that is a set of samples of a noisy complex carrier signal, with a carrier frequency \(\delta f\), amplitude \(1/L_T\), and phase \(\pi\delta f T_1\), samples at \(1/T_1\) [Hz], and that estimating the carrier frequency and phase of \(\mathbf{R}_1(\delta f)\) can provide a mismatch of Doppler frequency and phase estimates between those of the incoming signal and the receiver replica. In addition, when the phase of \(R_1^l(\delta f)\) is denoted by \(\phi\) as in Eq. (7) the accumulation of \(R_1(\delta f - \Delta_f)\) for all \((l\in\{0,1,…,R_T-1\})\) in Eq. (11a), i.e., the summation of LT consecutive short coherent ACF outputs, will have an overall phase \(\phi - \pi\Delta_fT_1\). Therefore, when the frequency error \(\Delta_f\) is estimated, the phase of the incoming signal \(\phi\) can be readily available from

\(\phi=∠R_1(0)=∠R_1(-\Delta_f)+\pi\Delta_f T_1.\)                                                                                           (12)

consequently, we can see thet estimating the frequency error \(\delta f\) is to search for a solution \(\Delta_f\) that maximizes \(│R_1(\delta f-\Delta_f)│\) such that

\({\hat{\Delta}}_f={\rm argmax}\Delta_f|R_1(\delta f-\Delta_f)|,\)                                                                                               (13)

which is equival to the frequency estimation of carrier signal sample \(\mathbf{R}_1(\delta f)=\{R_1^l(\delta f)|l=0,1,…,L_T-1\}.\)

Fig. 1. Validity of the approximation in Eq. (11b).

4. FREQUENCY AND PHASE ESTIMATION USING THE PROPOSED TECHIQUE

Ther are a number of techniques that can provide a solution for the problem in Eq. (13). For example, instead of searching for \(\Delta_f\), trivial techniques can be the Discrete Fourier Transform (DFT) and the FFT that can estimate the frequency of \(\mathbf{R}_1(\delta f)\) by finding the maximum frequency component such that

\(\widetilde{\delta}f=\overset{\rm argmax}{_{\delta f}} |FFT\{R_1(\delta f)\}|,\)                                                                                               (14)

which may require only \(L_T \log_2 L_T\) complex multiplications. The esimate \(\widetilde{\delta}f\) is asymptotically unbiased for a very large \(L_T\) (Hayes 1996), and the maximum error of the FFTbased scheme Eq. (14) is \(1/(2T)\) [Hz]. Since the frequency estimation error with the FFT (or DFT) is uniformly distributed in \([-1/(2T), 1/(2T)]\) [Hz], the absolute mean frequency estimation error \(1/(4T)\) [Hz] is not sufficiently accurate when T is not very large. In this section, however, we propose a reverse directional FRI technique to increase the frequency estimation accuracy when \(T\) is not very large.

4.1 Proposed Reverse Directional FRI Technique

In fact, the observation \(\mathbf{R}_1(\delta f)\) is the vector of samples of a signal with FRI (Vetterli et al. 2002), i.e., the samples contain only a finite number \(K(<∞)\) of innovation of a periodic signal. For an example, a low pass filtered input signal \(y(t)\) composed of \(T_P\)-periodic distinctive \(K\) Dirac delta functions with complex weights \(x_k (k=0,1,...,k-1)\) as shown in Fig. 2a, the FRI technique (Vetterli et al. 2002), is to estimate the delays \(t_k (k=0,1,...,k-1)\) of the Dirac delta functions in the input using the minimum number of samples of \(y(t)\). The estimation of \(t_k (k=0,1,...,k-1)\) is performed in the frequency domain, since the time domain samples of \(y(t)\) are effectively the samples of the linear sum of complex sinusoids. The delays can be accurately estimated from the zeros of the annihilating filter \(H\)[z] that annihilates the samples of the linear sum of complex sinusoids. And the amplitude xk is estimated by solving the Vandermonde system of linear equations (Vetterli et al. 2002). On the contrary, the problem in Eq. (13) is to estimate the frequency of a single \((K=1)\) complex sinusoid that is equivalently a single Dirac delta function in the frequency domain, which means that the problem in Eq. (13) requires reverse directional application of the FRI technique as shown in Fig. 2b and as investigated in Kim et al. (2011), for synthetic aperture radar (SAR) signal processing. In other words, \(\mathbf{R}_1(\delta f)\) can be regarded as a set of noisy time domain samples of a single (i.e., \(K=1\)) Dirac delta function in the frequency domain as

\(R_1^l(\delta f) = \left ( {A_q\sin(\pi\delta f T_1) \over \pi\delta f T_1} e^{j[\pi\delta f T_1+\phi]} \right ) e^{j2\pi\delta fl T_1} + n_1^l,\)                                                                                       (15)

where the sampling is performed at every \(t=lT_1\). Therefore, we can expect that the single Dirac delta function in the frequency domain is \(T_1^{-1}\)-periodic, and that the location and amplitude of the corresponding delta function in the frequency domain can indicate the carrier frequency and phase of the \(\mathbf{R}_1(\delta f)\), respectively.

Fig. 2. Reverse directional application of FRI for GPS.

Exploiting Blu et al. (2008) for the case of noisy sample input, \(\mathbf{R}_1(\delta f)\), when the input is arranged into a Toeplitz matrix

\(\mathbf{A}=\left [ \begin{matrix} R_1^L(δf),&R_1^{L-1}(δf),&…,& R_1^0(δf) \\ R_1^{L+1}(δf),& R_1^L(δf),& …,& R_1^1(δf)\\ ⋮ & ⋮ & ⋱ & ⋮ \\ R_1^{L_T-1}(δf),& R_1^{L_T-2}(δf),&…,& R_1^{L_T-L-1}(δf) \end{matrix} \right ], \)                                                                       (16)

where \(L≥K\), an annihilating filter

\(H(z)=\sum\limits_{k=0}^1 h_kz^{-k}=1-e^{-j2\pi f_1T_1} z^{-1}\)                                                                                       (17)

that minimizes the total least squares

\(\overset{\min}{_{H}}∥AH∥^2\)                                                                                                           (18)

under a constraint \(\Vert H\Vert^2=1\) provides the frequency estimate \(\hat{\delta}f\), i.e., \(f_1=\hat{\delta} f\). The coefficients of the annihilating filter \([h_0,h_1,…,h_k]\) can be found from the eigenvector of \(\mathbf{A}^T\mathbf{A}\) corresponding to the smallest eigenvalue, which is related to Pisarenko’s method (Pisarenko 1973, Tufts & Kumaresan 1982).

To increase the accuracy of the frequency estimate \(\hat{\delta}f\) in the presence of heavy noise in the signal sample input, there can be two ways to increase the noise robustness of the algorithm in Eq. (18). One way is to increase \(L_T\), which is a trivial solution, and the other way is to apply the iterative Cadzow denoising technique used in Blu et al. (2008). The first step of the technique is to perform a singular value decomposition (SVD) of \(\mathbf{A}\) such that

\(\mathbf{A}=\mathbf{USV}^T,\)                                                                                                         (19)

where U, S, and V are unitary, diagonal with elements in descending order, and unitary matrices of size \([(L_T-L)×(L+1)]\), \([(L+1)\times(L+1)]\), and \([(L+1)\times(L+1)]\), respectively. When the \((K+1)\)-th largest diagonal element of S is not sufficiently smaller than the \(K\)-th largest diagonal element, Cadzow’s iterative denoising technique sets the \(L-K+1\) smallest diagonal elements of \(\mathbf{S}\) to zero to build a new diagonal matrix \(\mathbf {S^{(1)}}\) in the first step, where \((·)^{(i)}\) denotes the \(i\)-th iteration of the Cadzow denoising technique. In the second step, the technique constructs \(\mathbf {A^{(1)}}\) using \(\mathbf {S^{(1)}}\) as

\(\mathbf{A^{(1)}}=\mathbf{US^{(1)}V}^T,\)                                                                                                       (20)

The resulting \(\mathbf {A^{(1)}}\) is not Toeplitz anymore, so the technique obtains the best Toeplitz approximation of \(A^{(1)}\) by averaging the diagonals of \(\mathbf {A^{(1)}}\) in the step. The above three steps are usually iterated Ni(>1) times until the \((K+1)\)-th largest diagonal element of \(\mathbf{S}^{(N_i)}\) becomes \(R_{TH}(\gg 1)\) times smaller than the \(K\)-th latgest diagonal element. From experimental observation, \(R_{TH} \geq 50\) is sufficiently high. When it is found that the Cadzow denoising is completed after \(N_i\)-th iteration, \(\mathbf{A}^{(N_i)}\) can be re-arranged to form a matrix \(\mathbf{B}\) of size \([L_T \times 2]\) as

\(\mathbf{B}=\left [ \begin{matrix} R_1^1(\delta f)',&R_1^2(\delta f)',&…,& R_1^{L_T-1}(\delta f)' \\ R_1^0(\delta f)',&R_1^1(\delta f)',&…,& R_1^{L_T-2}(\delta f)' \\ \end{matrix} \right ], \)                                                                             (21)

where \(R_1^l(\delta f)'\) is the denoised \(R_1^l(\delta f)\) by the iterative Cadzow denoising technique, and the eigenvector of \(\mathbf{B}^T\mathbf{B}\) corresponding to the smallest eigenvalue provides the coefficients of the annihilating filter \([h_0, h_1]\). From Eq. (17) and for \(K=1\),

\(\hat{\delta}f=f_k|_{k=1}=-{1\over 2\pi T_1} \angle \left ({-h_1 \over h_0} \right )\)                                                                                             (22)

is the final estimate of the frequency \(\delta f\) buried in \(\mathbf{R}_1(\delta f)\) found by the proposed reverse directional FRI technique.

Once the frequency \(\delta f\) is estimated, from Eq. (12), the relative phase of the input signal to the receiver generated carrier signal can be easily obtained as

\(\hat{\phi}=\angle \left ( \sum\limits_{l=0}^{L_T-1} R_1^l (\delta f)e^{-j(2l+1)\pi\hat{\delta} fT_1} \right )\)                                                                                         (23)

which is the arc-tangent of the frequency compensated and averaged In-phase and Quadrature-phase consecutive short coherent ACF outputs. Note that the proposed technique can estimate both the frequency and phase difference between the incoming signal \(r(t)\) and the receiver replica signal, and that when there is no frequency, i.e., \(\delta f=0\), the phase estimate of the proposed technique is similar to that of the conventional phase estimation technique.

4.2 Realization and Receiver Complexity

A drawback of the proposed technique for the pull-in search maybe the high computational cost for multiple SVD’s in the iterative Cadzow denoising. To reduce the computational cost, computationally efficient SVD algorithms such as Riemannian SVD (R-SVD) and Golub-Reinsch SVD (GR-SVD), introduced in Chan (1982) can be applied as a solution. For an SVD of the matrix A of size \([(L_T-L)×(L+1)]\), R-SVD and GR-SVD require

\(M_1= [ 4(L_T-L)^2+22(L+1)^2 ](L+1)\)                                                                                        

for \(L\)  sufficiently smaller than \(L_T/2,\)                                                                                       (24a)

\(M_2=\left [ 4(L_T-L)^2 + 8(L_T-L)(L+1)+9(L+1)^2 \right ] (L+1)\)                                                                        

for \(L \approx L_T/2,\)                                                                                                                   (24b)

respectively (Golub & van Loan 1996). Assuming \(L=L_T/4\) for Eq. (24a) and \(L=L_T/2\) for Eq. (24b), and \(L\gg 1\) for both cases, we can find that \(M_1 \approx M_2/3\). Therefore, 1   여기 많이 빠짐

\(M_{CD} \approx N_i(18\alpha_l^2-8\alpha_l+4)\alpha_l L_T^3,\)                                                                                             (25)

where \(N_i \leq 5\) is usually observed in experimentations (Blu et al. 2008). Notice that \(M_{CD}\) in Eq, (25) is a monotonously increasing function with respect to \(\alpha\), so that smaller \(\alpha\) results in less computational cost. Fig. 3 shows the results of \(10^4\) Monte Carlo simulations for a fixed \(R_{TH}=50\) and for various \(\alpha\) to estimate the total number of multiplications \(M_{CD}\) in the proposed technique employing iterative Cadzow denoising with respect to GPS signal strength. Figs. 3a and b show the simulation results for \(L_T=16\) and \(L_T=32\), respectively. In both figures, \(1/4\leq \alpha \leq 1/2\) results in a similar total number of multiplications, but \(\alpha=1/8\) results in the largest total number of multiplications, which is because the number of iterations Ni is increased for smaller \(\alpha\). Overall, the computational cost is the lowest or near lowest when \(\alpha=1/3\), so \(L\approx L_T/3\) is a choice leading to a computationally efficient application of the iterative Cadzow denoising algorithm. Therfore, in the following, we use \(L\approx L_T/3\) for the proposed technique.

Fig. 3. Complexity comparison.

5. PERFORMANCE ANALYSIS AND NUMERICAL SIMULATIONS

In this section, we compare the performance of the proposed technique and other techniques including the conventional technique (Tsui 2005).

Exploiting the analysis on the uncertainty of the estimate with the FRI technique for a single Dirac delta function, the variance of the unbiased frequency estimate \(\hat{\delta}f\) Eq. (22) is found in Blu et al. (2008) and is lower bounded by the Cramer-Rao lower bound (CRLB) as

\(\sigma_{\delta f}^2 \geq {3T_p^{-2} \over \pi^2 L_T^2PSNR},\)                                                                                                       (26)

where \(T_p^{-1}=T_1^{-1}\) is the period of the frequency domain spectrum, and \(PSNR\) represents the peak signal-to-noise ratio such that

\(PSNR = {L_T|E[R_1(0)^l]|^2 \over \sigma_n^2}.\)                                                                                                   (27)

Fig. 4. Performamce comparison

Fig. 4a shows the standard deviation of frequency estimation errors of the proposed technique and other techniques with respect to the carrier to noise denisty ratio \((C/N_0)\). The performance of the proposed techique is almost the same to the theoretical bound CRLB when \(C/N_0\) is not too weak. Note that when the FFT-based frequency resolution governs the estimation accuracy so that the 64-point FFT (i.e., \(L_T=64\)) has twice the constant resolution error of the 123 -point FFT as shown in Fig. 4a, where ‘\(N_F\)-pt FFT (\(N_F\) ms)’ in the legend denote the \(N_F\)-point FFT from \(N_F\) ms of data. Similarly, the performance of the serial Doppler frequency search scheme with coherent correlation length \(T_{CO}=N_{CO}T_1\) ms if is lowerbounded by the Doppler frequency search step size, where \(N_{CO}=20\). In the serial Doppler frequency search simulations, we chose the search step size of \(10^3/64\) Hz to compare with the 64 -point FFT-based technique. When the curve fitting technique (Zeng & Li 2010) is applied to the serial Doppler frequency search result obtained with search step size of 1/20 Hz, the performance is improved and shows that it has the same slope to the CRLB in the moderate and high \(C/N_0\) region. Fig. 4a shows the performance of the conventional techique (Tsui2005) also; the conventional frequency discriminator uses two consecutive ACF outputs, with coherent integration length \(T_1=1\) ms, to compute

\(\delta f_{conv} = {1 \over 2(N-1)\pi T_1} \sum\limits_{m=0}^{N-2}[I_m Q_{m+1}-I_{m+1}Q_m],\)                                                                                 (28)

where \(I_m\) and \(Q_m\) represent the in-phase and quadrature-phase component of the ACF output at \(t=mT_1\), and the moving average length of \(N-1\) is assumed. Note that the performance of the conventional technique is the worst among the techniques shown in Fig. 4 for \(C=N_0\) not high enough, which is due to the short coherent integration length \(T_1\) for each ACF output, and that we can expert 3 dB more performance gain wher the coherent integration for each ACF output is doubled. Note also thet the ideal frequency estimate in Eq. (28) does not include any measurement error due to the approximation funstion. After some algebraic manipulations, the frequency error variance can be found as

\(var\{\delta f_{conv}\}={1 \over 4(N-1)\pi^2T_1^2C/N_0} \left (1+{1 \over 2T_1C/N_0} \right ).\)                                                                               (29)

Notice that when \(C/N_0\)=45 dB-Hz the frequency errors of the proposed technique for \(L_T=16, L_T=32\) and the curve fitting technique are about 1 Hz, 0.3 Hz, and 0.8 Hz, respectively, and the errors increase much slowly than in the conventional technique as the \(C/N_0\) decreases. The resulting phase estimation performnce of the techniques compared in Fig. 4a shows some similarity as shown in Fig. 4b. Note that the phase estimation is performed with a shorter data length than the data length used for frequency estimation. For example, when the frequency is estimated with \(N_F\)-point FFT \((N_F)\) ms, the maximum frequency error is \(10^3/2 N_F\) Hz with which the phase error can increase to \(\pi\) [rad]. Since a phase estimation with an arc-tangent function may result in a large arithmetic error in the presence of noise, we limit the maximum possible phase estimation error by taking a shorter data length than that for the frequency estimation. In Fig. 4b, a legend, for example, ‘64-pt FFT (64 ms)+32 ms’ denotes the phase estimation with 32 ms of data when frequency estimation is performed with \(L_T=64\), and it can be found that a shorter data length for phase estimation results in a smaller phase estimation with the proposed technique and with the curve fitting technique is much better than with other techniques for moderate and strong signals.

Fig. 4c shows the computational complexity (i.e., the number of complex multiplications) of the frequency estimation techniques, such as \(N_F\)-point FFT-based technique \((N_F=64)\), serial Doppler frequency search techique with \(T_{CO}=20\) ms and search step size \(10^3/N_F\) Hz, conventional technique Eq. (28) with \(N=20\), and the proposed technique with \(L_T=16\) ms and 32 ms, discussed in this section. When the sampling rate \(f_s=2R_c=2\times 1.023\) MHz, the computational costs for the frequency estimation techniques can be expressed as

\(M_{FFT}=N_F(f_sT_1+\log_2 N_F)\)                                                                                           (30a)

\(M_{SS}=f_s T_{co} N_{DF}\)                                                                                                           (30b)

\(M_{conv}=Nf_sT_1 \log_2(f_sT_1)\)                                                                                                 (30c)

\(M_{FRI}=f_s T_1 L_T + M_{CD,}\)                                                                                                   (30d)

where \(M_{FFT}, M_{SS}, M_{conv}\), and \(M_{FRI}\) represent the computational complexity of the \(N_F\)-point FFT based technique, the serial search technique with coherent integration length \(T_{co}\) testing \(N_{DF}=N_F/2=32\) Doppler frequency hypotheses tested with a smaller coherent integration length \(T_1\) (i.e., 500 Hz), the conventional technique Eq. (28) with \(N=20\), and the proposed technique, respectively. Notice that computational complexity for computing the arc-tangent function is not counted in the express ions in Eq. (30), and that, for high C=N0, \(M_{FFT}=64\times(2046+\log_2 64)=131328\), \(M_{SS}=20\times 2046 \times 32=1350360\), \(M_{conv}=20\times 2046\times \log_2 2046=450062\), \(M_{FRI}=2046\times 16+5\times [4\times 11^2+22\times 6^2]\times 6+[4\times 15^2+22\times 2^2]\times 2=72992\) for \(L_T=16\), where \(N_i=5\), and \(M_{FRI}=2046\times 32+2\times [4\times 21^2+22\times 12^2]\times 12+[4\times 31^2+22\times 2^2]\times 2=191704\) for \(L_T=32\), where \(N_i=2\), which match the simulation results for \(C/N_0\)=45 dB-Hz shown in Fig. 4c.

Obviously, the total complex multiplications of the curve fitting technique are larger than the serial search technique due to the additional curve fitting process. Furthermore, the three figures in Fig. 4 shows that the proposed technique can produce much more accurate Doppler frequency and phase estimates than other techniques and that the computational complexity of the proposed technique is only around \(2\times10^5\) multiplications for \(L_T\) ms of data, while other techniques have multiple times of estimation error and computational cost in the GPS pull-in state.

6. CONCLUSIONS

It has been found that a reverse directional finite rate of innovation technique can be applied to GPS pull-in search and improves the accuracy and computational cost of the pull-in search in comparison to the conventional techniques. Analytical expressions are derived to obtain the theoretical estimates of performance improvement and computational complexity with the proposed technique. The simulations have demonstrated that the proposed technique has a performance improvement when SNR is moderate or high, and that the improvement becomes significant as the resolution of frequency estimation of the proposed technique increases linearly as SNR increases, which is not possible with the conventional FFT-based frequency estimation techniques. In addition to the improvement of the frequency estimation, the proposed technique contributes to the accuracy in the phase estimation, and, therefore, improves the phase tracking performance of a receiver.