Weiss-Weinstein bound of frequency estimation error for very weak GNSS signals

Tightness remains the center quest in all modern estimation bounds. For very weak signals, this is made possible with judicial choices of prior probability distribution and bound family. While current bounds in GNSS assess performance of carrier frequency estimators under Gaussian or uniform assumptions, the circular nature of frequency is overlooked. In addition, of all bounds in Bayesian framework, Weiss-Weinstein bound (WWB) stands out since it is free from regularity conditions or requirements on the prior distribution. Therefore, WWB is extended for the current frequency estimation problem. A divide-and-conquer type of hyperparameter tuning method is developed to level off the curse of computational complexity for the WWB family while enhancing tightness. Synthetic results show that with von Mises as prior probability distribution, WWB provides a bound up to 22.5% tighter than Ziv-Zaka\"i bound (ZZB) when SNR varies between -3.5 dB and -20 dB, where GNSS signal is deemed extremely weak.


INTRODUCTION
The primary objective in a navigation system is to estimate a parameter vector .In the first category of problems, the parameters are assumed constant or static during the observation interval.There are in turn, two subcategories.The first one models the static parameter as a set of deterministic or nonrandom yet unknown variables while the second models the parameter set as random with an a priori probability distribution function (pdf).Within the collection of deterministic bounds, there are local bounds including Cramér Rao bound (CRB) (p480, Cramér, 1949) and Bhattacharyya bound (Bhattacharyya, 1946) which describes the smallest achievable estimation error variance for maximum likelihood (ML) estimators.Since these two local bounds are good performance indicators only in large-sample or high signal-to-noise ratio (SNR) settings, they are often called asymptotic or small-error bounds since if we plot them against SNR, they are tight in the high SNR region in the sense that the estimated variance approaches but never reaches the true variance.The difference between CRB and Bhattacharyya bound lies in the choice of score function (Richmond & Horowitz, 2015), which is a function of  to describe the undulation of its log-likelihood function : CRB uses first derivative of  with respect to the true parameter, , whereas Bhattacharyya bound extends to use a mixture of unlimited higher order derivatives.The latter choice of score functions help to improve the bound's ability to predict error covariance in the asymptotic region but, it is, however, this very nature of characterizing undulation of loglikelihood function around the neighborhood of the true parameter that fails these local bounds in low SNR regions.The SNR in these cases will be so low that there are no obvious peaks in the log-likelihood function and it is highly probable that wrong peaks which do not correspond to the ML estimate are picked and estimation error is thus wrongly reported as low.In this regard, a very natural remedy is to include measures of test points distinctively different from the true value, , in the score function.Test points are estimated parameters supposed to be 'outliers' that lead to the largest possible mean square error (MSE).The first of this kind is due to Barankin (1949) who in turn was inspired by the work of Riesz in functional analysis.A more engineer-friendly derivation of Barankin bound is provided by McAulay and Hofstetter (1971).Barankin bound endorses a much tighter bound in low-tomedium SNR regions and thus can be seen as a large-error or global bound in stark contrast to CRB.We adopt a delimited approach to describing distinctive bound behaviors over the entire SNR range from Weiss & Weinstein's work (1983, Figure 3) and paraphrase it in FIGURE 1 (a) where the delimiters can only be defined qualitatively.The involved bounds are materialized as CRB, and two bounds belonging to the next subcategory: Weiss Weinstein bound (WWB) and Bayesian CRB (BCRB), whose definitions are deferred to Sections 3. In FIGURE 1(b), a simulation of ML and maximum a posteriori (MAP) frequency estimators with 10,000 trials, CRB, BCRB, and WWB is shown.One can immediately observe that CRB tends to underestimate MSE while Barankin bound (region II) tends to partially overcome this shortcoming of CRB.However, there are two things left to desire.Firstly, when SNR gets even lower, i.e., in region I of FIGURE 1(a), Barankin bound is unable to characterize correctly the error variance.Secondly, even in the medium-to-low SNR region, the vanilla Barankin bound is still not tight enough (regions II & III of FIGURE 1(a)).For the first problem, the observation that in this area, noise instead of useful signals dominates the observed samples endorses modelling the estimated variable as random, i.e. with a proper a priori pdf.Since now we are actually scoring undulation of an a posteriori function and developing a bound for MAP estimators, this subcategory of static performance bounds is called Bayesian bounds, which can be divided into two classes: the Ziv-Zakaï (ZZB) family and the Weiss-Weinstein (WWB) family.For the second problem, the most efficient solution is to select test points wisely; this is the central topic of this work and will be revealed in the successive sections.
Unlike all the other bounds we have mentioned so far which are derived from covariance inequalities, the ZZB family is derived from the probability of error in a binary hypothesis testing problem.It includes the original form developed by Ziv and Zakaï (1969), and later developments by Seidman (1970), Chazan et al., (1975), Bellini & Tartara (1974, 1975), Weinstein (1988), andBrown &Liu (1993).The most important development along the line is the extended Ziv-Zakaï bound (EZZB) by Bell et al. (1997), which extended the bound to include vector parameters with arbitrary a priori distributions while its predecessors can only accept a scalar parameter with a uniform a priori pdf.
In contrast, the WWB family is derived from Cauchy-Schwarz inequality (Ibragimov & Hasminskii, 1981;Lehmann, 1983).Bobrovsky & Zakaï (1975) developed the Bayesian version of Barankin bound as a prelude.Weiss and Weinstein (Weiss & Weinstein, 1985;Weinstein & Weiss, 1988) developed the vanilla WWB.Other members of this family include Bayesian Bhattacharyya bound (Van Trees, 1968, pp.148, problem 2.4.23),Bobrovsky-Zakaï bound (BZB) (Bobrovsky & Zakaï, 1976), and Bayesian Abel lower bound (BAB) (Renaux et al., 2008).Renaux et al. (2008) also proposed a constrained optimization framework to uncover the relationship between Bayesian CRB (BCRB), BZB, BAB, and vanilla WWB in the WWB family where the proposed constraints provide an enlightening view on how to choose the test points.A more recent treatment on WWB can be found in (Chaumette et al., 2017(Chaumette et al., , 2018)).As a sidenote, the second problem of interest in navigation is the estimation of a sample function of a stochastic process which can be modelled in either a correlation approach or a state dynamics approach where the recursive version of BZB has been developed by Fritsche et al. (2018) and recursive WWB developed by Reece et al. (2005) and Rapoport & Oshman (2004, 2007a, 2007b).
Finally, the problem of bounding carrier tracking and in particular, frequency tracking using modern and theoretically tight bounds such as WWB and comparison between WWB and ZZB has been overlooked, although code tracking using ZZB under Gaussian assumptions has been investigated (Emmanuele, 2012).There are asymptotic derivations for frequency estimation, for the standard case, high dynamic case (Mcphee et al., 2023a) which accounts for up to 100g acceleration long the line-of-sight (LOS), and a generic non-line-ofsight (NLOS) case (Lubeigt et al., 2020).The misspecified CRB (Ortega, 2023;McPhee, 2023b;Lubeigt et al., 2023) have been thoroughly investigated -misspecified because the assumed signal model unintentionally and practically differs from the received signal, resulting in bias.Therefore, the misspecified CRB describes the performance of biased estimators at low SNR, where RMSE is characterized by bias (which may be large for many cases) combined with the misspecified CRB.However, because the ambiguity function of the Doppler parameter is much wider than that of the delay, the convergence of the Doppler estimate is much earlier than the operating range of the navigation system, e.g. for GPS between 2-3 dB earlier, providing a theoretical limit of frequency estimation at low SNR and hence, for new applications such as Doppler positioning using low earth orbit (LEO) satellites, where the CN0 could be as low as 10 dB-Hz, these asymptotic bounds are not fit enough to fully recover the threshold, Barankin,.Renaux (2008) investigated bounding frequency using WWB, yet with simplified signal model and overlooking the effects of pseudorandom code and assumed a Gaussian a priori pdf.All these inspire us to develop a tight bound and in particular, WWB for frequency tracking.The reason for choosing WWB is three-fold: (a) with its test points finetuned, WWB is expected to outperform ZZB in some situations, and (b) WWB is essentially free of regularity conditions (Weiss & Weinstein, 1988), which are haunting other bounds and it is in this regard, WWB is expected to be applied to a broader range of estimation problems.
In this work, we develop WWB for a specific GNSS problem, i.e. frequency tracking and compare it to the other state-of-the-art bound and in particular ZZB for very weak GNSS signals, to clarify the preference of these two bounds.A modern professional GNSS receiver can track signals at a carrier-to-noise-density ratio (C/N0) down to 10 dB-Hz (Manzano-Jurado et al., 2014;Ren & Petovello, 2017;Pany & Eissfeller, 2006) and a typical lower bound of 15 dB-Hz (Soloviev et al., 2009) is required in space applications such as Space Service Volume (IS-GPS-200K, 2019, section 3.3.1.6.1).There exists modern mass-market receiver capable of steady tracking at a received power as low as -197 dBW (U-blox, 2018) which, for a typical front-end with a noise density between -203 and -205 dBW/Hz, means the working C/N0 is as low as 6 to 8 dB-Hz.However, tracking at this level of C/N0 can be extremely challenging, almost impossible.Therefore, input C/N0 between 10 and 45 dB-Hz is assumed in our analysis to cover various GNSS applications.
The remaining sections are organized as follows.Firstly, the estimation problem is formulated using the signal and noise models introduced in section 2 where the a priori distribution, von Mises distribution is briefly discussed.The development of the Weiss Weinstein bound follows in section 3 and the analytical form of GNSS frequency WWB is presented.Detailed development is provided in Appendix A. In section 4, this WWB is optimized over its two model hyperparameters, namely the exponential index { !} and the test points {ℎ !} where  = 1, … , ;  is the number of test points.A divide and conquer optimization procedure including only two steps is introduced to ease optimization over a possibly large sets of parameters { !, ℎ !}.In particular, a new set of test points is proposed by combining side lobe positions of the estimates and a vector of points evenly spanning the frequency uncertainty region.This optimization can offer a performance advantage of up to 0.95 dB, equivalent to a 19.8 % increase in predicted root mean square error (RMSE), compared to the WWB evaluated using legacy test points.Together with the synthetic simulations in section 5, we demonstrate important properties of the WWB, namely, that a.) with a priori pdf of von Mises distribution, WWB is a tighter bound than ZZB in the SNR region of -3.5 ~ -20 dB, where GNSS signal is deemed weak, but inferior to ZZB in predicting the exact thresholding SNR, and b.) WWB and ZZB converge both asymptotically (in high SNR) and in the no-information region.We also investigate the dependence of the developed WWB on a priori pdf, input SNR and integration time and found that when von Mises distribution is concerned, the effects of the location parameter,  and the concentration parameter, , of WWB, can be decoupled within a limited range.In summary, contributions of this paper are expected to be two-fold: (1) Better prediction of RMSE of any circular frequency estimator in the mid-to-low transition region of SNR, by developing a frequency WWB using von Mises distribution as prior, and (2) A divide-and-conquer type of parameter tuning method to choose test points of WWB to achieve a tighter bound.

PROBLEM FORMULATION
Firstly, we introduce notations including the signal and noise models.Then, we formulate the problem as developing a tight bound on estimation error of the normalized digital frequency.We also confirm that very weak GNSS signals do occur within the transition and even lower SNR region that WWB is trying to characterize.Finally, we introduce von Mises as the chosen prior distribution.

Notations
We will use bold and lower-case letters to denote vectors, bold and upper-case letters to denote matrix, and italic letters to denote scalars. , {⋅} means expectation of the curly bracketed function with respect to joint pdf of  and .For an observation vector  of length  and parameter vector  of length , we will denote a priori pdf as (), joint pdf of  and  as (, ) and likelihood function as (; ) or (|) and a posteriori pdf as (|).

Signal and Noise Models
We consider receiving a complex analytical signal (Betz, 2015) () = √( − ) exp{2 %& ( − )} (1) transmitted over a carrier frequency  %& with delay  . is the power in the data/pilot component and  is the baseband signal with bandwidth  including both spreading code and navigation data.This signal model implicitly implies a narrowband assumption on  which means the effect of Doppler dilatation on the baseband signal (Lubeigt, 2020) can be neglected.Mixing in the receiver frontend, using an estimate of the received carrier frequency  F %& , performs frequency translation and adjustment of () as where the residual or the intermediate frequency (IF) after downconversion is Frequency tracking then estimates the residual frequency  '& .It is assumed that the receiver is approximately synchronized to the overlaid baseband waveform and carrier so correlation may be done by integration and dump at an interval of 1/ '() .The correlation samples, denoted as  * , can be summarized as (Betz, 2015, eqn. (18.6), (18.8), or (18.10)) where each sample comprises a deterministic part,  * , defined as with time epochs  * ≜ / '() and phase  ≜ −2 %& .The phase  is assumed to be constant within one integration interval.The complex Gaussian noise,  * , is modeled as a sampled circular complex Gaussian white noise with variance 2 , as where  *,3 is the Dirac delta function and ⟨⋅⟩ represents expected value of the bracketed item with respect to the noise pdf.The signal and noise are bandlimited to ± 4 /2 Hz.Therefore, the signal-to-noise ratio (SNR) is defined as which means its carrier-to-noise density ratio (/ 5 ) is herein.The amplitude , first appearing in (4) mainly summarizes the code correlation result, baseband power , and impacts of presence of data symbol; we will examine this signallevel term shortly and build a sketch view of the order of magnitude it takes by using SNR and  , .This completes our definition of the signal and noise models.The remaining question is how weak the target signal is that we are estimating?The answer is directly related to the necessity of developing a bound and in particular, a tighter bound in the transition and even extremely low SNR region.SNR is related to input C/N0 by input bandwidth and coherent integration time.With the assumed working C/N0 ranges from 10 to 45 dB-Hz and for a typical baseband processing chain with 4.5 MHz intermediate frequency (IF) bandwidth (Joseph, 2010, Figure 2) plotted in FIGURE 2, this means the SNR varies between -13 ~ 22 dB at the input of phase-locked loop (PLL) if the coherent integration time is 5 ms and -20 ~ 15 dB at the input of frequency locked loop (FLL) if the coherent integration time is 1 ms.Therefore, for our current problem of frequency estimation, we assume that SNR varies between -20 and 15 dB.Within this range, the transition does happen and this justifies the necessity to seek an estimation lower bound tighter than current bounds shown in FIGURE 1(b) for a wide range of GNSS applications subject to weak signals.

Problem Formulation
To simplify symbolization, we denote normalized digital frequency as which varies between [−, ] radians and finally, we frame our estimation problem as developing a Weiss Weinstein bound (WWB) on estimating  using the samples where  is determined using preset  and  , in (6); for purpose of theoretical investigation,  , is set unity.A few more words about what normalization of  '& by  '() involved in  means.Take GPS L1 C/A as an example, the primary code period is 1ms and thus  '() is 1 kHz.This means that the extreme  '& that can be estimated, usually called pull-in range, is ± 500 Hz, corresponding to  = ± .Of course, for low SNR,  > 1 integrations should be summed up coherently to achieve a reasonable accuracy of frequency estimation.For 1/ !9/ =5 ms, the estimable  '& , becomes ±100 Hz (Kaplan, 2006, pp. 170), also corresponding to  = ±.Hence, normalization facilitates different integration period used in GNSS.Note that we aim to develop an error bound for MAP estimation without having to go through the intimidating calculation of the vanilla MAP estimator; 'vanilla' means no approximations are applied and therefore, we are developing a bound for the full allowable range (Kaplan, 2006, pp.170), i.e. [− '() /2,  '() /2] ; approximated such as differential schemes of MAP or ML usually means the pull-in range halved or even quartered.

A Priori Distribution of Circular Frequency: von Mises Distribution
The a priori pdf of the circular frequency is assumed to be von Mises.This distribution was studied in directional statistics and is a close approximation to wrapped normal distribution, which is in turn the circular/periodic analog of the normal distribution (Mardia, 1972).Therefore, it is particularly appropriate to characterize the parameter of interestcircular frequency since it is wrapped around [−, ] (Nitzan, 2016).The pdf of von Mises distributed parameter, , is given as (Fisher, 1993) where  5 () is the modified Bessel function of the first kind of order zero (Mardia, 1999, pp.36);  and  are the mean and concentration parameters, respectively.In particular, large  means smaller dispersion; when  = 0 this distribution reduces to a uniform distribution bounded between [−, ].Note that von Mises distribution () becomes a close approximation to a wrapped normal distribution with parameters  and 1- @ ()/ 5 () as its defining parameters when  approaches infinity (Mardia, 1999, pp.38),where  @ () is the modified Bessel function of the first kind of order one.However, other approximations exist.For the current problem where  seldom exceeds 20, von Mises can be approximated by a normal distribution with mean  and a variance equal to −2 ln( @ ()/ 5 ()) (Berens, 2009).

WEISS-WEINSTEIN BOUND FOR GNSS FREQUENCY ESTIMATION
Signals with changing parameters are often encountered in many practical signal processing applications and endorse many ways of characterization.Such changes, known in literature in a broader sense as 'change-points' (Bacharach, 2017), are often considered random and can be reflected by shifts in the parameters of its distribution.Among all the possible estimation schemes, ML or MAP estimators often preferred for their good statistical properties, easy circuit implementation and most important of all, optimum performance in appropriate conditions.However, in the context of change-point estimation, certain regularity assumptions, usually used to prove the asymptotic normality of the ML or MAP estimators, are obviously not fulfilled, this is especially due to the discrete nature of the parameters of interest.Consequently, the study of such estimator performance requires a specific analysis.One of the accessible ways of closed-form performance evaluation is to use lower bounds with less stringent requirements on regularity and as we just mentioned in the introduction, WWB is an exact example of this family with less regularity conditions and in particular a lower bound which does not require the differentiability of its likelihood.
In this section, Weiss Weinstein bound (WWB) is developed and its link to other modern bounds is explained.
For an arbitrary estimation problem with an observation vector  and parameter vector  of length , the generic WWB is (Weiss & Weinstein, 1985;Weinstein & Weiss, 1988) where  ∈ ℝ %×% whose element will be defined shortly and  is defined as where each  !∈ ℝ D×@ and therefore  ∈ ℝ D×% . is a chosen parameter indicating the total number of test points  !whose selection method will be dealt with in Section 4.3.Intuitively, the appealing aspect of WWB comes from the fact that CRB is evaluated only around the (assumed) best estimate and thus has only one test point, while WWB is evaluated over a much wider uncertainty range of the parameter space which can be far from the best estimate, thus grasping large errors while SNR is small.An element of matrix  with  and  as its row and column index, respectively, is defined as where 0 <  !< 1, 0 <  E < 1, and function (;  @ ,  , ) is the joint likelihood ratio with as its definition.For our problem, only one parameter is to be estimated and thus the parameter vector  reduces to a scalar  and each column of  now reduces to a scalar ℎ !, and  = [ℎ @ , ℎ , , … , ℎ % ] is now a row vector of length  .With these in mind, after expanding the nominator of , we arrive at and the  matrix can be evaluated using a priori and a posteriori pdf 's.

Links to Other Bounds
Some other bounds, for example, Bayesian Cramér-Rao bound (BCRB) and Bobrovsky-Zakai bound (BZB), can be obtained as a special case of Weiss Weinstein bound.BCRB is obtained with  = 1 and  !→ 1 > and  E → 1 > when ℎ !→ 0 and ℎ E → 0 (Weinstein & Weiss, 1988).BZB is obtained when  = 1.Specifically, the Bayesian information matrix (BIM)  G for an observation vector  of length  and parameter vector  of length  is (Van Trees, 1968) and the BRCB for this problem is thus defined as readily.After some substitutions and manipulations, BCRB for this problem can be formulated as and its derivation is provided in Appendix B. We will use BCRB as one of the two benchmarks (the other is Ziv-Zakai bound) in the simulations to highlight the bounding performance of the developed bound.

TWO-STEP OPTIMIZATION OF WEISS WEINSTEIN BOUND
The advantage of WWB lies in its power to predict the threshold region of RMSE while CRB or Bayesian CRB cannot.The developed bound is a maximization of the matrix  >@  ) over a combination of two sets of parameters, {ℎ !} and { !},  = 1, … , , as in Equation ( 10).The row vector  = (ℎ @ , … , ℎ % ) is often called 'test points' in the literature (Weinstein, 1988).Test points are estimated parameters supposed to be outliers that lead to large RMSE.Therefore, evaluation of WWB is also an optimization procedure.For estimating a normalized frequency in GNSS, test points take values between [−, ].A typical value of the vector size, , can easily exceed 10 (single sided, Weinstein, 1988), which means we are facing an optimization problem including about 20 parameters (two sided) to be optimized over, a quite intimidating task.Fortunately, we have some rules-of-thumb (Weinstein, 1988;Brown & Liu, 1993;Bell et al., 1997) to follow.We propose that {ℎ !} and { !} need not be optimized simultaneously.In the sequel, we will solve this maximization problem using a divide and conquer strategy.Firstly, we fix  and optimize over { !}; then we use the optimized { !} to further optimize the bound over .But, before all these happen, it is helpful to explain the reasons why the WWB is better than other Bayesian bounds such as BCRB in predicting the threshold.The explanation helps in understanding the fundamental problem of optimization, i.e., how to select test points.

Test Point: The Reason Why WWB Outperforms BCRB
As mentioned in section 3.2, BCRB can be obtained from WWB if the length of the test point vector is 1 and the single test point, ℎ @ , approaches to zero.Therefore, we will use BCRB as an illustrative example for its simplicity.Suppose now we have a noisy signal and let us plot the logarithm of its a posteriori probability in FIGURE 4 to illustrate the behavior of MAP estimates, which the bounds we discussed so far are trying to bound.Parts of the log-likelihood that does not depend on  ‹ is deliberately omitted and here  = 0 is assumed without loss of generality.These simplifications are solely for better visual effects, since logarithm of a posteriori probability differs from the loglikelihood by a varying, i.e. not constant "noise floor" dictated by the a priori distribution.This noise floor will blur the main lobes and side lobes of the ambiguity surface once the noise gets higher, thus hiding BCRB's characteristic dependence on side lobe positions.As a side note, if the a priori distribution is uniform, the a posteriori log-likelihood reduces to log-likelihood because the noise floor is simply lifted by a constant.We will stick to the loglikelihood for the current discussion.
The reason why WWB outperforms BCRB lies in the choice of test vector: the number and locations of points.BCRB use a single test point extremely close to the main lobe while WWB use multiple test points, including the one used for BCRB computation.FIGURE 4(a) shows a noise-free ambiguity surface and in (b) we show three realizations of the corresponding noisy ambiguity surface for three values of SNR and  = 20 which means the integration time is now 20 ms for GPS L1 C/A.For SNR = 5 dB, the main lobe peaks are always correctly picked by the maximization and indexing operation, but for increasing levels of noises, the same operation tends to pick wrong peaks.This is the reason underlying the RMSE's thresholding behavior.During the transition of SNR from medium to low values, BCRB's choice of a single test point is slow to reflect this change in the log-likelihood function, since sometimes, the peak is still correctly picked, as shown in the second plot of FIGURE 4(b), when 2 out of 3 trials succeed in estimating the true value.In contrast, WWB uses multiple test points such that it has fewer chances of missing the gross errors or incorrect estimates, but still, the locations and number of test points need optimizing and the next two subsections are dedicated to optimization over WWB's two hyper parameters, namely the exponential index { !} and the test points {ℎ !} where  = 1, … , ;  is an arbitrary positive integer.Historically, selection of test points appears in prior art as assertions and several rules of thumb are obtained for specific applications.We will extend this trial-and-error way, but deal with the exponential index and test point positions separately to find an optimum set of test points in section 4.3.We will use the legacy test points in section 4.2 to optimize WWB over { !}, as the first step of the two-step optimization procedure.

Step 1: Optimization over {𝒔 𝒊 }
Previous work about WWB left the choice of { !} almost untouched and only an arbitrary value such as 0.5 is magically provided (Weinstein & Weiss, 1988), although for such a simple problem with only a single test point, i.e.  = 1, it is easy to prove that  = 0.5 is indeed the maximizer.Aiming for a more general case, we propose to use a brute-force method by numerically evaluating WWB to fix a vector of  ! that maximizes WWB.As a first step, we will use fixed, i.e. legacy test points (Xu, 2001;DeLong, 1993;Xu et al., 2004).These test points are composed of: (1) small values extremely close to the main lobe peak, e.g.ℎ != 0.01, 0.001, such that  + ℎ ! will be very close to the main lobe peak position, ; this is to ensure WWB behaves correctly at high SNR (or asymptotic region), i.e.WWB will converge just like BCRB and also, CRB in this region; (2) frequency points corresponding to the side lobe peaks of the  I (); for our problem we use all frequency points corresponding to positive side lobe peaks between [0, ].This is to ensure WWB's behavior in medium-to-low, i.e. thresholding and low SNR regions.Again, we use a noise-free MAP estimation function  I () from Equation (30) as an illustrative example.The legacy test points (groups 'C' and 'S') are shown in FIGURE 5.
With these test points, solution of WWB, formulated as Equation ( 11), is now an optimization problem only over { !}, i.e.
and as mentioned before, this could lead to a tough problem because the total number of parameters in the end could be as high as  , , where  is the total number of test points (type 'C', 'S', and 'E'), eventually rendering the optimization intractable.Therefore, we adopt a simplified assumption: i.e.,  !'s are all equal-valued.Now the only problem is to select a value within range of (0,1).
Although here the a priori pdf is parameterized with a particular choice,  = 0 and  = 2, the result is also true for other combinations of  and .We defer the discussion of how a priori pdf affects the developed WWB to section 5.In the following sections,  !=  = 0.5 will be assumed.

Step 2: Optimization over Test Point Vector, 𝐡
As the second step in the two-step optimization procedure, this subsection deals with optimization over test points {ℎ !} and presents the proposed test points.
The proposed test points are shown in FIGURE 5. We propose a trio of arguments, (C, S, E), to describe the configuration of WWB's test point: 'C' is the number of points close to the main lobe peak, 'S' the number of points corresponding to side lobe peaks, and 'E' the number of points evenly spaced between [0.1, ].The 'E' component is the additional improvement we propose to append to the legacy test point vector.Following the thinking that more points lead to better performance (Xu et al., 2004;DeLong, 1993;Xu, 2001), the proposed test points include: (1) small values extremely close to the main lobe peak; this part is identical to the legacy ones, (2) bunch of side lobe peaks selected exactly as those included in the legacy vector, and (3) an extra vector evenly spaced between [0.1, ].Points ( 1) and ( 2) are exactly the legacy points in section 4.2.Reusing (1) from the legacy points is due to the fact that these two points are enough for a good asymptotic behavior, and reusing (2) from the legacy points is because frequency points corresponding to side lobe peaks bring the biggest performance boost with least number of additional points.
As an example, FIGURE 7(a) shows an example of WWB using the legacy and proposed test points when  = 20 and increasing  values of von Mises pdf.With the introduced trio notation, in FIGURE 7(a), WWB (,,M,5) is the legacy configuration and WWB (,,M,@5) is the proposed configuration.The evenly spaced points proposed is aimed to provide additional chances of capturing 'wild' peaks during mid-to-low SNR transition, thus increasing the value of WWB in this transition region.The improvement is 0.95 dB at SNR = −4 dB, equivalent to a 19.8 % ™= 10 /.12 $/ /(2)š increase in predicted RMSE.Take GPS L1 C/A for example, this means a 3 Hz improvement in bounding accuracy:  '() ⋅ | @ −  , |/(2), where  '() = 1000 Hz.Note here  !( = 1,2) represents the RMSE's of the two compared bounds and should be converted from dB to radians first.This performance improvement can be better perceived with a close-up, i.e.FIGURE 7(b).The reason to stop continuing adding in evenly spaced points is that with increasing number of points, the performance margin gets smaller and smaller and finally diminishes.

SIMULATIONS
This section illustrates the bounding performance of WWB with the proposed model parameters through numerical simulation.The simulations consider a typical GNSS receiver architecture presented in FIGURE 2 indicating SNR varies between -20 to 15 dB for frequency estimation.

Numerical Challenges
Computation of WWB, as evident from Equation ( 16), requires ( , ) time, where  is the length of the test point vector and  represents the big O notation of computational complexity in time (in fact, also in space).This prevents a quick evaluation of WWB.For large inference network using signal bound as input (Manlaney, 2020), evaluating WWB in a quick and simple manner is still what we are after.Therefore, a reasonable number of test points are needed and it is meaningful to investigate whether all the positive side lobe peak points should be included in the test point vector.This is especially important when the data length, , which roughly two times the number of side lobe peaks, increases to a large number of for example, 1000 for extremely low SNR conditions. was assumed to be 20 in the previous examples and will still be used here to help select the appropriate number of side lobe peaks.We are not forced to use all these points, so the nine positive side lobe points, starting from the point closest to the main lobe peak, are progressively added to the test vector and the results are shown in FIGURE 8(a)-(f), with different combinations of  and .
To examine only the effects of the side lobe peaks and for better visualization, FIGURE 8 compares WWB (,,9,5) where  ∈ {1,3,5,7,9}, i.e.only half of the test points, to the legacy configuration, WWB (,,M,5) .In each case, increasing number of test points help to obtain a tighter bound, but this benefit becomes less significant with the increase of  (FIGURE 8(a)-(c)).We conclude that a properly small value of the concentration parameter of von Mises distribution, i.e., assuming frequency distributes like a flatter bell shape suffices.Changes of  seldom leads to changes of bound values (FIGURE 8(d) v.s. (e), and also, FIGURE 8(b) v.s.(f)).With all these in mind, we will proceed with a parameterization of  (,,M,@5) , i.e., 2 points close to the main lobe peak, 9 side lobe peak points, and 10 extra points evenly spaced between [0.1, ] in the ensuing simulations.

WWB Dependence on A Priori Distribution
In discussing FIGURE 8, we know by observation, that the effects of two parameters of our a priori distribution,  and  can be decoupled, but the conclusion is drawn with limited combinations of  and . ) and  ∈ (3, +∞ ).Although theoretically, this result is true for any large values of  > 3,  cannot be infinitely large since the computation of WWB will experience numerical problems for the a priori pdf now takes the shape of a Dirac delta function.The von Mises random samples are generated using open-source packages (Berens, 2009).Though the result is obtained with specific SNR and sample count , it is generally true for other SNR and  values, as seen in FIGURE 9(b).Note that when  is small (< 3) and  ∈ (− -, , -, ), the bounds in FIGURE 9(a) and (b) go in reverse directions.However, this divergence is negligible since the difference does not exceed 0.5 dB.

WWB Dependence on Input SNR and Coherent Integration Time K
Previously, FIGURE 6 shows WWB values change with input signal SNR and number of input samples,  .As defined in section 2.2,  is the number of primary PRN code integration interval and thus is also the coherent integration length in primary PRN code length.With the observations in section 5.2 that small  values (< 3) is of particular interest, we plot this dependence in FIGURE 10 with  = 1 and  = 0.The results show that larger  or larger SNR always means smaller estimation errors and for weak signal deep into −10 dB or even lower,  required to reduce estimation errors increase tremendously.This coincides with our knowledge about this type of estimators.

WWB Compared to Other Modern Bounds
Bayesian Cramér-Rao bound (BCRB) and Ziv-Zakaï bound (ZZB) are used as benchmarks to check whether WWB is indeed the tightest bound, in most cases, among modern bounds when GNSS frequency estimation is concerned.A brief development of BCRB is provided in Appendix B. FIGURE 11 depicts WWB against BCRB and MAP results with different  values while  is defaulted to zero.by assuming that the observation vector  consists of  independent and identically distributed (i.i.d.) samples from a Gaussian distribution with zero mean and variance  , as mentioned in signal model defined in Section 2 and the phase  is constant during the integration/summation.The square bracketed term is easily identified as the digital Fourier transform of the observation vector and the last term, − lnx2 5 ()y, does not depend on the parameters of interest and thus can be removed in the simulation.
In can be seen in FIGURE 11, for high SNR, both BCRB and WWB converges to MAP results and for low SNR, these two bounds flatten out at the prior variance.Within the entire SNR range, WWB is a tight lower bound for MAP error.The threshold occurs for both WWB and MAP at the same SNR while BCRB fails to do so.For low SNR, we can also observe that WWB converges to MAP error.Since we are interested SNR exceeding -20dB, we do not show comparison between WWB with BCRB or MAP in exceedingly small SNR say -20 ~ -45 dB, where WWB and BCRB eventually converges to MAP since there is no extra information provided by signal samples now; the only information the estimator can depend upon is prior knowledge of the estimated values.
As a more competent competitor, ZZB is given without proof as (Bell et al., 1997) where  RST denotes the minimum probability of error obtained from the optimum likelihood ratio test and Λ((ℎ)) is a non-increasing function of ℎ obtained by filling in any valleys in (ℎ).For the current frequency estimation problem, the right-hand side of Equation ( 35) can be computed, with some simple manipulations, as where  & is the Fisher information matrix computed as (equation 82, Appendix A) and  6 , is variance related to the a priori distribution of .Φ() is the tail probability of the standard normal distribution at point .and Γ W () is the incomplete gamma function defined by X>@  (37) herein.Comparison between WWB and ZZB are shown in FIGURE 12.The WWB configuration is again the proposed WWB (,,M,@5) .In all cases WWB and ZZB converge both asymptotically (in high SNR) and in the no-information region.In FIGURE 13, we compare both WWB and ZZB to MAP.Both of them are tight bounds in the sense that the difference between them and MAP is small within the Barankin and threshold regions.One interesting observation emerges.ZZB tend to be larger (thus tighter) than the WWB very close to the threshold, but smaller (thus looser) than the developed WWB for lower values of SNR.The larger the value of , the less this difference is (FIGURE 13).Existing work tends to discuss threshold behavior of bounds in a spanned region hence we divide the threshold region (e.g.qualitatively defined in Weiss's 4-region model mentioned in FIGURE 1) into two sub-regions: Region 1 (SNR = -10 ~ -3.5 dB): WWB > ZZB.WWB outperforms ZZB.The largest performance gain occurs at -8 dB (C/N0=22 dB-Hz) and amounts to 1.5 dB.Take GPS L1 C/A for example, this means a 67 Hz improvement in bounding accuracy: , where  '() = 1000 Hz.This is not a trivial improvement (22.5% = 10 $.2 $/ /(2) ) in bounding accuracy (ZZB: 28.8 Hz, WWB: 95.7 Hz).Note also here that  !( = , ) represents the RMSE's of the two compared bounds and should be converted from dB to radians first.
It is true that ZZB outperforms WWB in terms of predicting the exact thresholding point where outliers surge triggered by deteriorating SNR; this is a rather interesting observation.Meanwhile it is also true that within the majority of the threshold region, WWB succeeds in providing tighter bounds than ZZB.This alternating performances of the two most-discussed Bayesian bounds are repeatedly confirmed in FIGURE 12 where the same prior pdf with different parameters is used.
The phenomenon means tradeoff has to be made between these two Bayesian type bounds.For extremely low SNR (from -3.5 dB down to -20 dB), WWB is a better performance indicator, while in a relatively limited range of SNR (between -3.5~0 dB), ZZB is preferred.

CONCLUSIONS
In this work, we developed a Weiss-Weinstein bound (WWB) for circular frequency tracking error in GNSS receivers, by deriving the bound and proposing a new method to select a new set of test points that is required by WWB.We choose WWB in preference over other modern Bayesian bounds for its being free from regularity conditions or requirements on prior distribution.Later in the evaluation part we conclude that this divide-and-conquer type of selecting test points help to tighten the bound in the mid-to-low SNR region (-20 ~ -3.5 dB).
We proposed to use von Mises distribution, a cyclic counterpart of Gaussian distribution as prior distribution to specifically deal with cyclic nature of the estimated variable, i.e. carrier frequency in GNSS signals.This helps to tighten the bound in 'no signal' zone where signals are extremely weak (close to SNR=-20 dB).
We evaluated our methods through Monte Carlo simulations and compare it to Ziv-Zakaï bound and confirm that the developed bound can be tighter than Ziv-Zakai bound by up to 22.5% when SNR varies between -3.5 dB and -20 dB.
Future efforts may be invested in deriving the WWB for joint frequency and phase estimation and using the bound as a constraint condition in developing innovative code and carrier estimators (Pany, 2010) in weak GNSS conditions.
This result allows for a handy notation of the internal integral: where  5 ! is the Euler's number raised to the power of  -and With this identity, it turns out that Equation (40), which is a double integration, can be split into two separate integrations: one of them equals exactly ∫  678 !(|) 8 !(| + ℎ -)

𝐱
, and the other is . The latter will be evaluated in A.1.2.

APPENDIX B. BCRB
We can disassemble the  G term of Equation ( 27) into two parts (van Trees & Bell, 2007) as Specifically, the subscript 'D' stands for 'data' and 'P' for 'a priori', which states that  H represents the contribution from data and  I the contribution from a priori pdf of the parameter to be estimated.They are defined as from which we can readily identify the part bracketed by the outer curly braces of Equation ( 80) as the Fisher information matrix (FIM), denoted as  $ .Assuming the first sample is taken at  = 0, the FIM for our current problem can be calculated as (Wu et el., 2006) The last equation mark is due to the fact that we recognize the righthand side of this equation includes the modified Bessel function of the first kind and order 1 (Mardia, 1999, Appendix 1, eqn (A.1)).Substitute Equations ( 83) and ( 84) back into (79), we arrive at Equation (29) as BIM.

FIGURE 1
FIGURE 1 Thresholding behavior of root MSE (RMSE) for ML & MAP estimators when SNR transits from medium to low values.(a) An illustrative and delimited bound model by Weiss & Weinstein with CRB, BCRB, and WWB.(b) A simulation of GNSS L1 C/A frequency estimation (tracking) with frequency pull-in range of [−500, 500] Hz using the same bounds and estimators from (a); von Mises, a circular analog of Gaussian distribution is assumed as a priori distribution for BCRB and WWB.Here the mean and concentration parameter of von Mises is zero and 20, respectively.

FIGURE 2
FIGURE 2A typical GNSS receiver processing gain example.

FIGURE 5
FIGURE 5 Legacy and proposed test points.Legacy test points compose only 'C' and 'S' points; the proposed test points include 'E' points while retaining all 'C' and 'S' points.

FIGURE 8
FIGURE 8 WWB ((,*,+) with increasing number of test points and different combinations of  and .
FIGURE 9 (a) shows the joint effects of these two parameters on the developed WWB with SNR = 0 dB and sample counts  = 100 in milliseconds.We can see that the effects of  and  can be decoupled within a limited range of  ∈ (−