Abstract
Tightness remains the primary goal in all modern estimation bounds. For very weak signals, tightness is enabled by appropriately selecting the prior probability distribution and bound family. While current bounds in global navigation satellite systems (GNSSs) assess the performance of carrier frequency estimators under Gaussian or uniform assumptions, the circular nature of frequency is overlooked. Of all bounds in the Bayesian framework, the Weiss–Weinstein bound (WWB) stands out because it is free from regularity conditions or restrictions on the prior distribution. Therefore, the WWB is extended for the current frequency estimation problem. A divide-and-conquer type of hyperparameter tuning method is developed to mitigate issues of computational complexity for the WWB family while enhancing tightness. Synthetic results show that for a von Mises prior probability distribution, the WWB provides a bound up to 22.5% tighter than the Ziv–Zakaï bound when the signal-to-noise ratio varies between –3.5 dB and –20 dB, where the GNSS signal is deemed extremely weak.
1 INTRODUCTION
Lower bounds are a fundamental tool for predicting, before actual implementation, the closeness of an estimator to optimality criteria such as the maximum likelihood (ML) or maximum a posteriori (MAP). The procedure for tightening bounds reveals, in a principled way, the principal components leading to smaller errors. This unique analytical capability offered by bounding theory is unmatched by repeated experimental or simulation approaches. Frequency tracking is pivotal for any global navigation satellite system (GNSS) receiver. This work aims to develop a bound on GNSS frequency tracking that is tighter than previous bounds to facilitate the discovery of new tracking loop designs for GNSS-challenged settings.
The primary objective in a navigation system including GNSS receivers is to estimate a parameter vector θ. In the first category of problems, the parameters are assumed constant or static during the observation interval. Within this category, there are two subcategories. In the first subcategory, static parameters are modeled as a set of deterministic or nonrandom yet unknown variables. In the second subcategory, the parameter set is modeled as random with an a priori probability distribution function (pdf). Within the collection of deterministic bounds, there are local bounds, including the Cramér Rao bound (CRB) (Cramér, 1946) and Bhattacharyya bound (Bhattacharyya, 1946), which describe the smallest achievable estimation error variance for ML estimators. Because these two local bounds show good performance only in settings with large samples or high signal-to-noise ratios (SNRs), they are often called asymptotic or small-error bounds. If we plot these bounds against the 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 the CRB and Bhattacharyya bound lies in the choice of score function (Richmond & Horowitz, 2015), which is a function of θ that describes the undulation of its log-likelihood function ϕ. The CRB uses the first derivative of ϕ with respect to the true parameter, θ, whereas the Bhattacharyya bound uses a mixture of unlimited higher-order derivatives. The latter choice of score function helps to improve the bound’s ability to predict error covariance in the asymptotic region; however, this feature of characterizing the undulation of the log-likelihood function around the neighborhood of the true parameter causes these local bounds to fail in low-SNR regions. In such cases, the SNR is so low that there are no obvious peaks in the log-likelihood function; thus, it is highly probable that false peaks that do not correspond to the ML estimate will be chosen and the estimation error will thus be falsely reported as low.
In this regard, a natural remedy is to include measures of test points distinctively different from the true value, θ, in the score function. Test points are estimated parameters assumed to be outliers that lead to the largest possible mean square error (MSE). The notion of test point was first reported by Barankin (1949), who was inspired by the work of Riesz in functional analysis. A more engineer-friendly derivation of the Barankin bound was provided by McAulay and Hofstetter (1971). The Barankin bound endorses a much tighter bound in low-to medium-SNR regions and thus can be seen as a large-error or global bound, in stark contrast to the CRB.
We adopt a delimited approach for describing distinctive bound behaviors over the entire SNR range from the work of Weiss and Weinstein (1983), as summarized in Figure 1(a), where the delimiters can only be defined qualitatively. The involved bounds are materialized as deterministic bounds, including the CRB and Barankin bound, and three bounds belonging to the second subcategory: the Weiss– Weinstein bound (WWB), Ziv–Zakaï bound (ZZB), and Bayesian CRB (BCRB), whose definitions are presented in Sections 3.1 (WWB), 3.2 (BCRB), and 5.4 (ZZB). In Figure 1(b), simulations of ML and MAP frequency estimators with 10,000 trials as well as the CRB, BCRB, Barankin bound, ZZB, and WWB are shown. One can immediately observe that the CRB tends to underestimate the MSE whereas the Barankin bound (region II in Figure 1(a)) tends to partially overcome this limitation of the CRB. However, shortcomings remain. First, as the SNR continues to decrease, i.e., in region I of Figure 1(a), the Barankin bound is unable to correctly characterize the error variance. Second, even in the low-to medium-SNR region, the conventional Barankin bound is still not tight enough (regions II and III in Figure 1(a)). For the first shortcoming, modeling the estimated variable as random, i.e., with an appropriate a priori pdf, is motivated by the observation that, in this area, noise instead of useful signals dominates the observed samples. In this case, we are actually scoring the undulation of an a posteriori pdf and developing a bound for MAP instead of ML estimators; thus, this subcategory of static performance bounds is denoted as Bayesian bounds, which can be divided into two classes: the Ziv–Zakaï family and the Weiss–Weinstein family. For the second shortcoming, the most efficient solution is to select appropriate test points; this aim is the central topic of this work and will be revealed in the following sections.
Thresholding behavior of RMSE for ML and MAP estimators when the SNR transits from medium to low values: (a) an illustrative and delimited bound categorization inspired by Weiss and Weinstein (1983); (b) a simulation of GNSS L1 C/A frequency estimation (tracking) with a frequency pull-in range of [−500, 500] Hz instantiating the same bounds and estimators from (a), where a von Mises distribution, i.e., a circular analog of the Gaussian distribution, is assumed as an a priori distribution for the BCRB and WWB, with the mean and concentration parameter equal to 0 and 5, respectively
Unlike the other bounds we have mentioned thus far, which are derived from covariance inequalities, the ZZB family is derived from the probability of error in a binary hypothesis testing problem. This bound family includes the original form developed by Ziv and Zakaï (1969) and later developments by Seidman (1970), Chazan et al. (1975), Bellini and Tartara (1974, 1975), Weinstein (1988), and Brown and Liu (1993). The most important development along this line is the extended ZZB (EZZB) reported by Bell et al. (1997), which extends 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 the Cauchy–Schwarz inequality (Ibragimov & Hasminskii, 1981; Lehmann, 1983). Bobrovsky and Zakaï (1975) developed a Bayesian version of the Barankin bound as a prelude. Weiss and Weinstein (Weiss & Weinstein, 1985; Weinstein & Weiss, 1988) developed the original WWB. Other members of this family include the Bayesian Bhattacharyya bound (Van Trees, 1968), 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 the BCRB, BZB, BAB, and original WWB in the WWB family, where the proposed constraints provide an enlightening view on how to select test points. A more recent treatment on the WWB has been provided by Chaumette et al. (2017, 2018). As a sidenote, in contrast to static parameter estimation, the second category of problems of interest in navigation is the estimation of a sample function of a stochastic process, which can be modeled via either a correlation approach or a state dynamics approach; in this context, the recursive version of the BZB was developed by Fritsche et al. (2018) and a recursive WWB was developed by Reece et al. (2005) and Rapoport and Oshman (2004, 2007a, 2007b).
Both the WWB and ZZB have long been applied to digital communications systems (Villares, 2005) and positioning (Graff et al., 2021) both in theory and in practice (Manlaney, 2020; Watson et al., 2020). In terms of GNSSs, Denis (2009) and Emmanuele (2012) investigated bounding GNSS code tracking with the ZZB. Closas (2009) compared the BCRB with the conventional CRB when bounding and differentiating errors of two-step positioning from errors of direct positioning. Later, Gusi-Amigó (2018) adapted the EZZB to bound direct positioning errors. Gifford et al. (2022) investigated time-of-arrival error bounds for wideband and ultra-wideband ranging systems provided with different a priori information on the multipath phenomena.
Finally, the problem of bounding carrier tracking and, in particular, frequency tracking using modern and theoretically tight bounds such as the WWB has been overlooked. There exist asymptotic, i.e. CRB, derivations for frequency estimation for the standard case and high dynamic case (Mcphee et al., 2023a) that accounts for up to 100g acceleration along the line of sight (LOS) and a generic non-LOS (NLOS) case (Lubeigt et al., 2020). The mis-specified CRB (Ortega, 2023; McPhee, 2023b; Lubeigt et al., 2023) has been thoroughly investigated; this bound is mis-specified because the assumed signal model unintentionally differs from the received signal in practice, resulting in bias. Therefore, the mis-specified CRB describes the performance of biased estimators at low SNR, where the root MSE (RMSE) is characterized by bias (which may be large for many cases) combined with the mis-specified CRB. However, because the ambiguity function of the Doppler parameter is much wider than that of the delay, the Doppler estimate converges much sooner than the operating range of the navigation system. For example, for the Global Positioning System (GPS), the Doppler estimate converges 2–3 dB earlier, providing a theoretical limit of frequency estimation at low SNR. Hence, for new applications such as Doppler positioning using low earth orbit satellites, where the carrier-to-noise-density ratio (C/N0) can be as low as 10 dB-Hz, these asymptotic bounds are not suitable to fully recover MAP estimation errors in the threshold, Barankin, and no-information regions shown in Figure 1(a). Renaux (2008) investigated bounding frequency using the WWB; however, this approach utilized a simplified signal model, overlooked the effects of pseudorandom code, and assumed a Gaussian a priori pdf, not accounting for the circular nature of the frequency variable.
All of these oversights for tightly bounding GNSS frequency tracking errors using appropriate forms of the WWB or ZZB open opportunities to investigate which bound is tighter. In the current work, we investigate these bounds by first proposing a method to obtain the tightest WWB for GNSS frequency tracking and then comparing this optimal WWB with the ZZB for very weak GNSS signals. Simulation results confirm two advantages of the WWB: (a) with its test points finetuned, the WWB outperforms the ZZB in some SNR regions, and (b) the WWB is essentially free of regularity conditions (Weiss & Weinstein, 1988), which hinder the use of other bounds, allowing the WWB to be readily applied to a broader range of estimation problems.
The major task of the current work is to develop the WWB and compare it with other state-of-the-art bounds, particularly the ZZB, for very weak GNSS signals to clarify the performance of these two bounds. A modern professional GNSS receiver can track signals at C/N0 values as low as 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). There exist modern mass-market receivers 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, corresponds to a working C/N0 as low as 6–8 dB-Hz. However, tracking at this level of C/N0 can be extremely challenging, almost impossible. Therefore, an input C/N0 of 10–45 dB-Hz is assumed in our analysis to cover various GNSS applications, with a focus on those working with very weak signals.
The remainder of this paper is organized as follows. First, the estimation problem is formulated using the signal and noise models introduced in Section 2, where a von Mises a priori distribution is briefly discussed. The development of the WWB follows in Section 3, and the analytical form of the GNSS frequency WWB is presented. A detailed development is provided in Appendix A. In Section 4, this WWB is further optimized over its two model hyperparameters, namely, the exponential index {si} and the test points {hi}, where i = 1,…, R and R is the number of test points. A two-step divide-and-conquer optimization procedure is introduced to ease optimization over a potentially large set of parameters {si, hi}. In particular, a new set of test points is proposed by combining side lobe peaks of the ambiguity surface 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 RMSE, compared with the WWB evaluated via legacy test points. Together with synthetic simulations in Section 5, we demonstrate important properties of the WWB, namely, that (a) with a von Mises a priori pdf, the WWB is tighter (22.5%) than the ZZB in the SNR range of −20 to −3.5 dB ( C/N0 of 10–26.5 dB-Hz), where the GNSS signal is deemed weak, but inferior to the ZZB in predicting the exact thresholding SNR, and (b) both the WWB and ZZB converge asymptotically in the high-SNR region and the no-information region characterized by a very low SNR. We also investigate the dependence of the developed WWB on the a priori pdf, input SNR, and integration time and find that when a von Mises distribution is used, the effects of the location parameter, μ, and the concentration parameter, κ, on the WWB can be decoupled within a limited range. In summary, the contributions of this paper include (a) an enhanced prediction of the RMSE of any circular frequency estimator in the mid-to-low transition region of SNR, achieved by developing a frequency WWB using a von Mises a priori distribution. The lower end of this region is deemed very weak for GNSS signal processing. The choice of the a priori pdf is further augmented by (b) a divide-and-conquer type of hyperparameter tuning method for selecting WWB test points to achieve a tighter bound while limiting computational complexity.
This open-source work is available at https://github.com/TMBOC/WWB, where the MATLAB codes used to generate the figures are provided.
2 PROBLEM FORMULATION
First, we introduce our notation in addition to the signal and noise models. Then, we formulate the problem of developing a tight bound on the estimation error of a normalized digital frequency. We also confirm that very weak GNSS signals do occur within the threshold and lower-SNR regions that we seek to characterize via the WWB. Finally, we introduce the von Mises distribution, for use as an a priori distribution, and reveal its connection with the more familiar normal or uniform distribution.
2.1 Notation
We use bold and lower-case letters to denote vectors, bold and upper-case letters to denote matrices, and italic letters to denote scalars. Ex,θ{·} indicates the expectation of the curly bracketed function with respect to the joint pdf of x and θ. For an observation vector x of length K and parameter vector θ of length D, we denote the a priori pdf as p(θ), the joint pdf of x and θ as p(x, θ), the likelihood function as p(x; θ) or p(x|θ), and the a posteriori pdf as p(θ|x).
2.2 Signal and Noise Models
We consider receiving a complex analytical signal (Betz, 2016):
1
transmitted over a carrier frequency fRF with delay τ. C is the power in the data/pilot component, and γ is the baseband signal with bandwidth B including both spreading code and navigation data. This signal model implicitly implies a narrowband assumption on γ, which means that the effect of Doppler dilatation on the baseband signal (Lubeigt, 2020) can be neglected. Mixing in the receiver front-end, based on an estimate of the received carrier frequency 
, leads to a frequency translation and adjustment of s(t):
2
where the residual or intermediate frequency (IF) after downconversion is 
. Frequency tracking then estimates the residual frequency fIF. It is assumed that the receiver is approximately synchronized to the overlaid baseband waveform and carrier such that correlation may be performed by an integrate-and-dump method at an interval of 1/fINT. Each of K arbitrary correlation samples, denoted as xk, can be summarized as follows (Betz, 2016, eqn. (18.6), (18.8), or (18.10)):
3
where each sample comprises a deterministic component, rk, defined as follows:
4
with time epochs 
and phase 
. The phase φ is assumed constant within a single integration interval. The complex Gaussian noise, nk, is modeled as sampled circular complex Gaussian white noise with variance 2σ2:
5
where δk,m is the Dirac delta function, overbars indicate the complex conjugate, and 〈·〉 represents the expected value of the bracketed item with respect to the noise pdf. Therefore, the K noise samples as a whole are centered and characterized by a covariance matrix 2σ2IK, where IK is a diagonal matrix of dimension K. The signal and noise are band-limited to ±fs/2 Hz. Therefore, the SNR is defined as follows:
6
and thus, the linear ratio (C/N0)linear is as follows:
7
Note that we append the subscript “linear” to the symbol C/N0 to differentiate the carrier-to-noise-density ratio on a linear scale, i.e. (C/N0)linear, from the carrier-to-noise-density ratio on a logarithmic scale, i.e. C/N0. The amplitude A, first appearing in Equation (4), primarily summarizes the code correlation result, baseband power C, and impacts of the presence of data symbols. Note that there also exists a ratio of energy per symbol to noise power spectral density, Es/N0 = SNR · fs/fINT, which incorporates the coherent integration gain. Below, we will examine the order of magnitude of C/N0 and its corresponding SNR within a typical GNSS receiver. This completes our definition of the signal and noise models. It now remains for us to define the weakness of the target signal that we are estimating. This feature is directly related to the necessity of developing a bound, particularly a tighter bound in the threshold and extremely-low-SNR regions.
The SNR is related to the input C/N0 by the input bandwidth and coherent integration time. For an assumed working C/N0 range of 10–45 dB-Hz and a typical baseband processing chain with 4.5-MHz IF bandwidth (Joseph, 2010) (plotted in Figure 2), the SNR varies between −13 and 22 dB at the input of the phase-locked loop (PLL) if the coherent integration time is 5 ms or between −20 and 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 the SNR varies between −20 and 15 dB. Within this range, a transition occurs, justifying the need to seek a lower estimation bound that is tighter than the current bounds shown in Figure 1(b) for a wide range of GNSS applications subject to weak signals.
A typical example of GNSS receiver processing gain
2.3 Problem Formulation
To simplify our notation, we denote the normalized digital frequency as follows:
8
which varies within [−π, π] radians. Finally, we frame our problem as developing a WWB on estimating θ using the following samples:
9
where A is determined from preset values for SNR and σ2 in Equation (6). For the purpose of theoretical investigation, σ2 is set to unity. To clarify the normalization of fIF by fINT in θ, we take the GPS L1 coarse acquisition (C/A) signal as an example. Here, the primary code period is 1 ms; thus, fINT is 1 kHz. Consequently, the maximum fIF that can be estimated, typically called the pull-in range, is ±500 Hz, corresponding to θ = ±π. For low SNRs, K > 1 integrations should be coherently summed to achieve a reasonable accuracy of frequency estimation. For 1/fINT = 5 ms, the estimable fIF is ±100 Hz (Kaplan, 2006), also corresponding to θ = ±π. Hence, normalization facilitates the use of different integration periods in GNSS. Note that we aim to develop an error bound for MAP estimation without having to apply the complex calculation of the conventional MAP estimator, where “conventional” indicates that no approximations are applied. Therefore, we are developing a bound for the full allowable range (Kaplan, 2006), i.e., [−fINT/2, fINT/2]; in approximation methods such as differential MAP or ML schemes, the pull-in range is usually halved or even quartered.
2.4 A Priori Distribution of Circular Frequency: von Mises Distribution
The a priori pdf of the circular frequency is assumed to be a von Mises distribution. This distribution has been evaluated in directional statistics and is a close approximation to the wrapped normal distribution, which, in turn, is the circular/periodic analog of the normal distribution (Mardia, 1972). Therefore, the von Mises distribution is particularly appropriate for characterizing the parameter of interest, i.e., the circular frequency, as it is wrapped around [−π, π] (Nitzan, 2016). The pdf of the von Mises distributed parameter, θ, is given as follows (Fisher, 1993):
10
where I0(κ) is the modified Bessel function of the first kind of order zero (Mardia, 1999) and μ and κ are the location/mean and concentration parameters, respectively. In particular, large κ values indicate less dispersion; when κ = 0, this distribution reduces to a uniform distribution bounded between [−π, π]. Figure 3(a) illustrates a von Mises pdf subject to μ = 0 for multiple κ values (i.e., 0, 1, 2, 5, 20) meaningful for a typical GPS L1 C/A use case, i.e., fINT = 1 ms and fIF ∈ [−500, 500] Hz. Figure 3(b) shows the mixed effects of the mean μ and concentration κ.
Von Mises pdf subject to the two defining parameters, i.e., mean μ and concentration parameter κ, for a typical GPS L1 C/A use case, i.e., fINT = 1 ms and fIF ∈ [−500, 500] Hz: (a) μ = 0 and κ ∈ {0, 1, 2, 5, 20}; (b) (μ, κ) ∈ {(0, 0), (π, 1), (π/2, 2), (0, 5), (−π/2, 20)}
Note that the von Mises distribution p(θ) closely approximates a wrapped normal distribution with μ and 1-I1(κ)/I0(κ) as its defining parameters when κ approaches infinity (Mardia, 1999), where I1(κ) is the modified Bessel function of the first kind of order one. However, other approximations exist. For the current problem in which κ seldom exceeds 20, the von Mises distribution can be approximated by a normal distribution with mean μ and a variance equal to −2ln(I1(κ)/I0(κ)) (Berens, 2009).
3 WWB FOR GNSS FREQUENCY ESTIMATION
Static random parameter estimation can be seen as a special case of estimating signals with changing parameters, which are often encountered in many practical signal processing applications and require various characterization methods. Such changes, known in the literature in a broader sense as “change-points” (Bacharach, 2017), are often considered random and can be reflected by shifts in distribution parameters. Among all possible estimation schemes, ML and MAP estimators are often preferred for their good statistical properties, easy circuit implementation, and, most importantly, optimum performance in appropriate conditions. However, in the context of change-point estimation, certain regularity assumptions, usually used to prove the asymptotic normality of ML or MAP estimators, are obviously not fulfilled; these violations are usually due to the discrete nature of the parameters of interest. Consequently, evaluating the performance of such estimators requires a specific analysis. One accessible method for conducting a closed-form performance evaluation is to use lower bounds with less stringent requirements on regularity; as mentioned in Section 1, the WWB is a prime example of this family free of regularity conditions and, in particular, a lower bound that does not require the differentiability of its likelihood. In this section, the WWB is developed, and its relation to other modern bounds is explained.
For an arbitrary estimation problem with an observation vector x and parameter vector θ of length D, the generic WWB is as follows (Weiss & Weinstein, 1985; Weinstein & Weiss, 1988):
11
where 
, whose elements will be defined shortly, and H is defined as follows:
12
where each 
and therefore 
. R is a chosen hyperparameter indicating the total number of test points hi, whose selection method will be discussed in Section 4.3. Intuitively, the appealing aspect of the WWB comes from the fact that the CRB is evaluated only around the (assumed) best estimate and thus has only one test point, whereas the WWB is evaluated over a much wider uncertainty range of the parameter space, which can be far from the best estimate, thus including large errors while the SNR is small. An element of matrix Q with i and j as its row and column indices, respectively, is defined as follows:
13
where 0 < si < 1, 0 < sj < 1, and the function L(x; θ1, θ2) is the joint likelihood ratio with the following definition:
14
For our problem, only one parameter is to be estimated, and thus, the parameter vector θ reduces to a scalar θ, each column of h reduces to a scalar hi, and H = [h1, h2,…, hR] becomes a row vector of length R. With these reductions in mind, after expanding the numerator of Q, we obtain the following:
15
and the Q matrix can be evaluated using the a priori and a posteriori pdf.
3.1 WWB for GNSS Frequency Estimation
In Appendix A, the WWB is developed in detail for frequency estimation in a generic GNSS setting. In short, [Q]ij can be computed as follows:
16
where the exponents μij,n, γij,n, n ∈ {1, 2, 3, 4}, μi, and γi are respectively defined by the following:
17
18
19
20
and:
21
22
23
24
and:
25
and:
26
μj and γj are obtained by simply replacing the subscript i with j in μi and γi accordingly.
3.2 Relations to Other Bounds
Other members of the same family, for example, the BCRB and BZB, can be obtained as a special case of the WWB. The BZB is obtained when s = 1 and is a Bayesian analog to the Barankin bound. The BCRB is obtained for R = 1, si → 1−, and sj → 1− when hi → 0 and hj → 0 (Weinstein & Weiss, 1988). Specifically, the Bayesian information matrix (BIM) JB for an observation vector x of length K and parameter vector θ of length D is as follows (Van Trees, 1968):
27
The BRCB for this problem is thus readily defined as follows:
28
After some substitutions and manipulations, the BCRB for this problem can be formulated as the inverse of JB, which is defined as follows:
29
A corresponding derivation is provided in Appendix B. We will use the BCRB as one of two benchmarks (the other is the ZZB) in the simulations to highlight the bounding performance of the developed bound.
4 TWO-STEP OPTIMIZATION OF THE WWB
The advantage of the WWB lies in its power to predict the thresholding behavior of the RMSE, which cannot be achieved by the CRB or BCRB. The developed bound is a maximization of the matrix HQ−1HT over a combination of two sets of hyperparameters: {hi} and {si}, with i = 1,…, R, as specified in Equation (11). The elements in the row vector H = (h1, …, hR) are often called “test points” in the literature (Weinstein, 1988). Test points are estimated parameters assumed to be outliers that lead to a large RMSE. Therefore, an evaluation of the WWB is also an optimization procedure. To estimate a normalized frequency in a GNSS, test points take values within [−π, π]. A typical value of vector size R can easily exceed 10 (single-sided) (Weinstein, 1988), corresponding to an optimization problem including approximately 20 parameters (two-sided) to be optimized, an intimidating task considering the sheer number of combinations of i and j involved in computing [Q]ij in Equation (13). Fortunately, we have some rules of thumb (Weinstein, 1988; Brown & Liu, 1993; Bell et al., 1997) to follow. We propose that {hi} and {si} need not be optimized simultaneously. In Sections 4.2 and 4.3, we will solve this maximization problem using a divide-and-conquer strategy. First, we fix H and optimize over {si}; then, we use the optimized {si} to further optimize the bound over H. However, it is helpful to first explain why the WWB is superior to other Bayesian bounds such as the BCRB in predicting the threshold. Such an explanation will aid in elucidating the fundamental problem of optimization, i.e., how to select test points.
4.1 Test Points: Why the WWB Outperforms the BCRB
As mentioned in Section 3.2, the BCRB can be obtained from the WWB if the length of the test point vector is 1 and this single test point, h1, approaches zero. Therefore, because of its simplicity, we will use the BCRB as an illustrative example. For a noisy signal, let us plot the logarithm of its a posteriori probability:
30
This probability is illustrated in Figure 4 to illustrate the behavior of MAP estimates, whose RMSE we are trying to approximate. Components of the log a posteriori probability that do not depend on θ are deliberately omitted, and here, κ = 0 is assumed without a loss of generality. These simplifications are solely for providing enhanced visual effects, as the logarithm of an a posteriori probability differs from the log-likelihood 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 as the noise increases, thus hiding the characteristic dependence of the BRCB on side lobe positions. As a side note, if the a priori distribution is uniform, the a posteriori log-likelihood reduces to the log-likelihood because the noise floor is simply lifted by a constant.
Ambiguity surface Ax(θ) with K = 20, for a true frequency set to 0.4π and a known phase of φ = π/6: (a) Ax(θ ) with xk = rk; (b) noisy Ax(θ) with xk = rk + nk
The reason why the WWB outperforms the BCRB lies in the choice of the test vector: the number and locations of points. The BCRB uses a single test point assumed to be extremely close to the main lobe whereas the WWB uses multiple test points, including the one used for the BCRB computation. Figure 4(a) shows a noise-free ambiguity surface, and Figure 4(b) shows three realizations of the corresponding noisy ambiguity surface for three values of SNR and K = 20, corresponding to an integration time of 20 ms for GPS L1 C/A. For SNR = 5 dB, the main lobe peaks are always correctly chosen by the maximization and indexing operations; however, for increasing levels of noise, the same operation tends to select false peaks, which leads to RMSE thresholding behavior. During the transition of SNR from medium to low values, the single test point for the BCRB is slow to reflect this change in the log a posteriori pdf; in some cases, the correct peak is chosen, as shown in the second plot of Figure 4(b), where two out of three trials succeed in estimating the true value.
In contrast, the WWB uses multiple test points such that it has a lower probability of missing gross errors or incorrect estimates; yet, the locations and number of test points still require optimization, and Sections 4.2 and 4.3 are dedicated to optimization over the two WWB hyperparameters, namely, the exponential index {si} and the test points {hi}, where i = 1,…, R and R is an arbitrary positive integer. Historically, the selection of test points has been reported in prior works as assertions, with several rules of thumb for specific applications. We will extend upon this trial-and-error method, 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 the WWB over {si}, as the first step of the two-step optimization procedure.
4.2 Step 1: Optimization Over {si}
Previous works on the WWB have primarily neglected the choice of {si}, providing only an arbitrary value such as 0.5 (Weinstein & Weiss, 1988); however, for a simple problem with only a single test point, i.e., R = 1, it is easy to prove that s = 0.5 is indeed the optimal value. Aiming for a more general case, we propose to use a brute-force method by numerically evaluating the WWB to fix a vector of si that maximizes the WWB. As a first step, we will use fixed, i.e., legacy, test points (Xu, 2001; DeLong, 1993; Xu et al., 2004), including the following:
Small values extremely close to the main lobe peak, e.g., hi = 0.01π, 0.001π, such that θ + hi is very close to the main lobe peak position, θ. These test points are included to ensure that the WWB behaves correctly at high SNR (or in the asymptotic region), i.e., the WWB will show a convergence similar to that of the BCRB and CRB in this region. This group of points is denoted “C.”
Frequency points corresponding to the side lobe peaks of Ax(θ); for our problem, we use all frequency points corresponding to positive side lobe peaks within [0.1π, π], half of the pull-in range but excluding the points close to the main lobe peak. These test points are included to assess the behavior of the WWB in medium-to-low, i.e., thresholding, and low-SNR regions. This group of points is denoted “S.”
Again, we use a noise-free MAP estimation function Ax(θ) from Equation (30) as an illustrative example. The legacy test points (groups “C” and “S”) are shown in Figure 5. The proposed test points belonging to group “E” will be introduced in Section 4.3.
Legacy and proposed test points
Legacy test points include only “C” and “S” points; the proposed test points include “E” points while retaining all “C” and “S” points.
With these test points, the maximization problem, formulated as Equation (11), is now an optimization problem over only {si}:
31
As mentioned before, this scenario could lead to a difficult problem because the final number of parameters could be as high as R2, where R is the total number of test points, eventually rendering the optimization computationally intractable. Therefore, we adopt a simplified assumption:
32
where all si values are equal. Now, all that remains is to select a value within the range of (0, 1). Because determining the derivatives of the WWB is almost analytically impossible, we use numerical methods by comparing WWBs with different s values. Specifically, we examine s ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. Results show that (1) s = 0.5 is the optimal value, and (2) s and (1-s) produce identical WWBs. An example of this result is shown in Figure 6. For better visualization, only WWBs with s = 0.1 and s = 0.5 are shown.
WWB with s = 0.1 and s = 0.5
For the sample numbers, K ∈ {20, 40, 60} and SNR ∈ [−20, 10] (dB); for the von Mises distribution, μ = 0 and κ = 2.
Although the a priori pdf is parameterized in this case with a particular choice, μ = 0 and κ = 2, this result holds for other combinations of μ and κ. We defer the discussion of how the a priori pdf affects the developed WWB to Section 5.2. In the following subsections, si = s = 0.5 will be assumed.
4.3 Step 2: Optimization Over Test Point Vector, h
As the second step in the two-step optimization procedure, this section deals with optimization over test points {hi} 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 test points: “C” is the number of points close to the main lobe peak, “S” is the number of points corresponding to side lobe peaks, and “E” is the number of points evenly spaced within [0.1π, π], which is half of the pull-in range but excluding the points close to the main lobe peak. The “E” component is the additional improvement we propose to append to the legacy test point vector. Therefore, following the viewpoint that more points lead to better performance (Xu et al., 2004; DeLong, 1993; Xu, 2001), the proposed complete set of test points includes the following:
small values extremely close to the main lobe peak, identical to the legacy points
side lobe peaks selected exactly as those included in the legacy vector
an additional vector evenly spaced within [0.1π, π]
Points in groups (1) and (2) are exactly the legacy points in Section 4.2. The points in (1) are reused from the legacy points because these points are sufficient for good asymptotic behavior. We reuse the points in group (2) from the legacy points because frequency points corresponding to side lobe peaks bring the biggest performance boost for the lowest number of additional points.
Figure 7(a) shows an example of the WWB using the legacy and proposed test points when K = 20 for increasing κ values of the von Mises pdf. With the introduced trio notation, in Figure 7(a), WWB(2,9,0) is the legacy configuration, and WWB(2,9,10) is the proposed configuration. The proposed evenly spaced points aim to provide additional chances of capturing “wild” peaks during the mid-to-low SNR transition, thus increasing the value of the WWB in this transition region.
(a) WWB: Legacy vs. proposed configurations for K = 20 and μ = 0, showing a maximal improvement of 0.95 dB at SNR = –4 dB; (b) zoomed-in view of (a) around SNR = −4 dB
The improvement is 0.95 dB for SNR = −4 dB, equivalent to a 19.8% 
increase in the predicted RMSE. Taking the GPS L1 C/A as an example, this result corresponds to a 3-Hz improvement in bounding accuracy: fINT · |RMSE1 – RMSE2|/(2π), where fINT = 1000 Hz. Note here that RMSEi(i = 1, 2) represents the RMSE of the two compared bounds and should first be converted from dB to radians. This performance improvement can be better perceived with a close-up view, as shown in Figure 7(b). The reason that we stop adding evenly spaced points is because as the number of points increases, the performance margin decreases and finally vanishes.
To summarize, a test point vector incorporating the legacy “C” and “S” points and the proposed “E” points is the optimal choice balancing bounding accuracy and efficiency. In Section 5, we evaluate the performance of this proposal, denoted as WWB(2,9,10) in Figure 7.
5 SIMULATIONS
This section illustrates the bounding performance of the derived frequency WWB (Section 3) with the proposed hyperparameters (Section 4) through numerical simulation. The simulations consider the typical GNSS receiver architecture presented in Figure 2, indicating that the SNR varies from –20 to 15 dB for frequency estimation.
5.1 Numerical Challenges
Computation of the WWB, as evident from Equation (16), requires O(R2) time, where R is the length of the test point vector and O represents the big O notation of computational complexity in time (and in space). This time requirement impedes a rapid evaluation of the WWB. For a large inference network using a signal bound as input (Manlaney, 2020), we still seek to evaluate the WWB in a quick and simple manner. Therefore, a reasonable number of test points is needed, and it is meaningful to investigate whether all positive side lobe peak points should be included in the test point vector. This is especially important when the data length, K, roughly twice the number of side lobe peaks, increases to a large number, such as 1000 for extremely low-SNR conditions. K was assumed to be 20 in the previous examples; we will retain this assumption here to help select the appropriate number of side lobe peaks. We are not required to use all of these points; thus, the nine positive side lobe points, starting from the point closest to the main lobe peak, are progressively added to the test point vector, and the results, i.e., bound values against SNR, are shown in Figures 8(a)–(f), subject to six different combinations of μ and κ.
WWB(2,n,0) for an increasing number of test points and different combinations of μ and κ (a) μ = 0, κ = 0 (b) μ = 0, κ = 2 (c) μ = 0, κ = 5 (d) μ = π/4, κ = 0 (e) μ = π/2, κ = 0 (f) μ = π/2, κ = 2
To examine only the effects of the side lobe peaks and for better visualization, Figure 8 compares WWB(2,n,0), where n ∈ {1, 3, 5, 7}, i.e., only half of the test point configurations, against the legacy configuration, WWB(2,9,0). In each case, an increased number of test points helps to obtain a tighter bound. The difference between WWB(2,7,0) and WWB(2,9,0) is 14 Hz (= fINT · |10^(−3.5/10) −10^(−4.5/10)|/(2π)), as shown in Figures 8(a) and (c), where SNR = −7 dB. Consequently, when the number of side lobe peaks, n, is between 7 and 9, similar bounding accuracy is obtained in some cases. In all other settings, the bounding difference is not negligible. Meanwhile, the benefit brought by an increased number of test points becomes less significant with increasing κ (Figures 8(a)–(c)). We conclude that a suitably small value of the concentration parameter κ for a von Mises distribution is sufficient, assuming that the frequency is distributed like an appropriate flatter bell shape. Changes in μ seldom lead to changes in bound values (Figure 8(d) vs. (e) and Figure 8(b) vs. (f)). With these findings in mind, we will proceed with a parameterization of WWB(2,9,10), i.e., with 2 points close to the main lobe peak, 9 side lobe peak points, and 10 additional points evenly spaced within [0.1π, π], in the ensuing simulations.
5.2 WWB Dependence on a Priori Distribution
In Figure 8, we observe that the effects of two parameters of our a priori distribution, μ and κ, appear to be decoupled in some cases: changes in these two parameters do not cause a marked difference in bound values. We will justify this observation by further considering the WWB dependence on extended 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 K = 100 in milliseconds. We can see that the effects of μ and κ can be decoupled within a limited range of μ ∈ (−π/2, π/2) and κ ∈ (3, +∞). Theoretically, this result is true for any large values of κ > 3; however, κ cannot be infinitely large because the computation of WWB will experience numerical problems as the shape of the a priori pdf now resembles a Dirac delta function. For this and later figures, random von Mises samples were generated via open-source packages (Berens, 2009).
WWB dependence on the von Mises distribution with parameters μ and κ: (a) K = 100, SNR = 0 dB; (b) K = 200, SNR = −10 dB
Although the above result was obtained for a specific SNR and sample count K, this result generally holds for other SNR and K values, as shown in Figure 9(b). Note that when κ is small (< 3) and μ ∈ (−π/2, π/2), the bounds in Figures 9(a) and (b) move in reverse directions. However, this divergence is negligible, as the difference does not exceed 0.5 dB or approximately 7 Hz for fINT = 1 kHz.
5.3 WWB Dependence on Input SNR and Coherent Integration Time K
Figure 6 showed that WWB values change with the input signal SNR and number of input samples, K. As defined in Section 2.2, K is the number of primary pseudorandom noise (PRN) code integration intervals and thus is also the number of coherent integration periods in primary PRN code length. With the observations in Section 5.2 that small κ values (< 3) are of particular interest, we plot this dependence in Figure 10 with κ = 1 and μ = 0.
WWB dependence on input SNR and primary PRN code period for (κ, μ) = (1, 0)
The results show that increases in K or SNR always result in smaller estimation errors; moreover, for weak signals of −10 dB or lower, the K values required to reduce estimation errors increase tremendously. This finding coincides with our knowledge about this type of estimator.
5.4 Comparison of the WWB with Modern Bounds
The BCRB and ZZB are used as benchmarks to assess whether the WWB is indeed the tightest bound, in most cases, among modern bounds for GNSS frequency estimation. A brief development of the BCRB is provided in Appendix B. Figure 11 depicts the WWB against BCRB and MAP results for different κ values with μ set to zero.
Comparison of the WWB vs. MAP and BCRB, with κ ∈{0, 1, 2, 5, 20} and μ = 0
The MAP result here is an average of 10,000 Monte Carlo trials. An individual run of MAP estimation can be performed as follows:
33
where ln p(x, θ) reduces to the following form using the likelihood p(x | θ) and prior distribution p(θ), i.e., von Mises distribution, as follows:
34
Here, it is assumed that the observation vector x consists of K independent and identically distributed samples from a Gaussian distribution with zero mean and variance σ2, as mentioned for the signal model defined in Section 2. The phase φ is also assumed to be 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, −ln(2πI0(κ)), does not depend on the parameters of interest and thus can be removed in the simulation.
Figure 11 shows that, for high SNR, both the BCRB and WWB converge to the MAP results; for low SNR, these two bounds flatten out at the prior variance. Throughout the entire SNR range, the WWB is a tight lower bound for MAP error. Thresholding behavior occurs for both the WWB and MAP for a given SNR, but the BCRB does not show thresholding behavior. For low SNRs, we also observe that the WWB converges to the MAP error. Because we are interested in SNRs exceeding –20 dB, we do not show a comparison between the WWB and BCRB or MAP error for exceedingly small SNRs such as –20 to –45 dB, where the WWB and BCRB eventually converge to the MAP error because no additional information is provided by signal samples at this point. Here, the only information the estimator can depend upon is prior knowledge of the estimated values.
As a more competent competitor, the ZZB is given without proof as follows (Bell et al., 1997):
35
where Pmin denotes the minimum probability of error obtained from the optimum likelihood ratio test and Λ(f(h)) is a non-increasing function of h obtained by filling any valleys in f(h). For the current frequency estimation problem, Equation (35) can be computed, with some simple manipulations, as follows:
36
where JF is the Fisher information matrix (FIM) computed as Equation (82) (Appendix A) and 
is the variance related to the a priori distribution of θ. Φ(z) is the tail probability of the standard normal distribution at point z, and Γa(z) is the incomplete gamma function defined herein as follows:
37
A comparison between the WWB and ZZB is shown in Figure 12. The WWB configuration is again the proposed WWB(2,9,10). In all cases, the WWB and ZZB both converge asymptotically (for high SNR) and in the no-information region.
WWB and ZZB with a von Mises distribution as the a priori pdf for a concentration parameter κ ∈ {0, 1, 2, 5, 20} and mean μ = 0
In all cases, the WWB is tighter than the ZZB for most SNR values, particularly in low-SNR regions. However, the ZZB seems to outperform the WWB within a limited range of SNR values in the THRESHOLD region. We examine this interesting behavior in Figure 13 by condensing the a priori pdf to a specific case to better visualize the difference between the WWB and ZZB.
WWB vs. ZZB for κ = 1, μ = 0
RMSEs for MAP estimates are also plotted as performance benchmarks.
In Figure 13, the ZZB tends to be larger (thus tighter) than the section of the WWB very close to the threshold, but smaller (thus looser) than the developed WWB for lower SNRs. As κ increases, this difference decreases. Existing work tends to explore the thresholding behavior of bounds in a spanned region; hence,we divided the thresholding regions (e.g., qualitatively defined in Weiss’s four region model mentioned in Figure 1, including NO-INFORMATION, BARANKIN, and THRESHOLD regions) into two subregions.
In subregion 1 (where SNR = –20 to –3.5 dB), the WWB outperforms the ZZB. The largest performance gain occurs at –8 dB (C/N0 = 22 dB-Hz), reaching 1.5 dB. Taking GPS L1 C/A as an example, this gain corresponds to a 28-Hz improvement in bounding accuracy: fINT · |RMSEWWB – RMSEZZB|/(2π), where fINT = 1000 Hz, RMSEZZB = −3.713 dB, and RMSEWWB = −2.209 dB. This is not a trivial improvement 
in bounding accuracy (ZZB: 67.7Hz, WWB: 95.7 Hz and a larger value means being more capable to predict the MAP error). Note here that RMSEi(i = WWB, ZZB) represents the RMSEs of the two compared bounds and should first be converted from dB to radians.
In subregion 2 (where SNR = –3.5 to 0 dB), the ZZB outperforms the WWB, with smaller performance gains compared with the previous case.
The finding that the ZZB outperforms the WWB in terms of predicting the exact thresholding point at which outliers surge, triggered by a deteriorating SNR, is interesting. Meanwhile, within the majority of the threshold region, the WWB succeeds in providing tighter bounds than the ZZB. These alternating performances of the two most-discussed Bayesian bounds are repeatedly confirmed in Figure 12, where the same prior pdf is used with different parameters.
This phenomenon indicates that a tradeoff must be made between these two Bayesian-type bounds. For extremely low SNR (–3.5 to –20 dB), the WWB is a better performance indicator, while in a relatively limited SNR range (–3.5 to 0 dB), the ZZB is superior.
6 CONCLUSIONS
In this work, we developed a WWB for circular frequency tracking error in GNSS receivers, by tightening the bound in three medium-to low-SNR regions, namely, no-information, Barankin, and threshold regions, with a focus on very low SNRs. We chose the WWB over other modern Bayesian bounds because it is free from regularity conditions or requirements on the prior distribution. From our evaluation, we concluded that the proposed divide-and-conquer approach for selecting test points can help to tighten the bound in the mid-to low-SNR region (–20 to –3.5 dB), equivalent to a C/N0 of 10–26.5 dB-Hz with a front-end bandwidth of 4.5 MHz. Specifically, we applied two changes to obtain a tighter WWB.
First, we used a von Mises distribution, i.e., a cyclic counterpart of the Gaussian distribution, as the prior distribution to specifically handle the cyclic nature of the estimated variable, i.e., carrier frequency in GNSS signals. This choice helps to tighten the bound in the no-information region, where signals are extremely weak (close to SNR = –20 dB).
Second, we augmented the legacy test points with “E” points, i.e., points evenly distributed within the pull-in range, to achieve a tighter WWB in the Barankin and threshold regions (SNR = –20 to –3.5 dB).
During the investigations, two useful observations emerged. (1) The inclusion of nearly all side lobe peaks of the log a posteriori pdf is necessary to provide a basic bounding accuracy. This accuracy may be further enhanced by including an appropriate number of points evenly distributed throughout the pull-in range in WWB evaluation. (2) The WWB is insensitive to changes in the parameters μ and κ of the von Mises pdf over a certain parameter range. The WWB is more sensitive to the coherent integration time K, input SNR, and choice of test points.
We evaluated our methods through Monte Carlo simulations, compared the WWB with the ZZB, and confirmed that the developed bound can be tighter than the ZZB by up to 22.5% when the SNR varies between –3.5 dB and –20 dB; SNRs at the lower end of this range are deemed very weak and present a challenge for GNSS frequency tracking.
Future efforts may focus on deriving the WWB for joint frequency and phase estimation and using the bound as a constraint condition for developing innovative code and carrier estimators (Pany, 2010) in very weak GNSS conditions.
HOW TO CITE THIS ARTICLE
Zhang, X., Zhan, X., Huang, J., Liu, J., & Xiao, Y. (2024). Weiss–Weinstein bound of frequency estimation error for very weak GNSS signals. NAVIGATION, 71(3). https://doi.org/10.33012/navi.654
CONFLICT OF INTEREST
The authors declare no potential conflicts of interest.
ACKNOWLEDGEMENTS
This work was supported by the National Key R&D Program of China (2022YFB3904401).
APPENDIX A WWB
The Q matrix of the WWB was provided in Equation (15) and is re-written here for ease of reference:
38
We will evaluate each expectation in the denominator and numerator separately, using the von Mises distribution as an a priori distribution.
A.1 Denominator of Equation (38)
The denominator is as follows:
39
By inspection, the denominator is a product of two expectations with identical likelihood ratio structures; only the test points (hi and hj) and thus the subscripts differ. We shall derive the first expectation and obtain the second expectation with a change of subscript. The first part of the product can be evaluated using a double integration as follows:
40
where Ωθ and Ωx are all possible values of θ and x, respectively. By the end of Appendix A.1.1, we will prove that this double integral can be evaluated via two independent integrals. The remaining derivations in this appendix will benefit greatly from this fact.
A.1.1 Inner Integral
The inner integral of Equation (40) can be evaluated as follows:
41
Let us set the following:
42
where ak is set as 
, which is a deterministic component (as we are deriving formulas of probability conditioned on the frequency to be estimated), and yk is thus a complex (circular) Gaussian variable because of the circular symmetric and Gaussian assumption for xk. Hence, we have the following:
43
and:
44
Substituting this result into Equation (41), we arrive at the following:
45
where:
46
Here, 
and 
represent the real and imaginary parts of a complex number. By substituting this expression into Equation (45), taking the natural logarithm, and recognizing 
, we obtain the following:
47
The summation in square brackets can be evaluated as follows:
48
This result allows for a convenient notation of the internal integral:
49
where 
is Euler’s number raised to the power of μi:
50
With this identity, we find that Equation (40), which is a double integration, can be split into two separate integrations: one integration equals exactly 
, and the other equals 
. The latter will be evaluated in Appendix A.1.2.
A.1.2 Outer Integral
The outer integral of the double integration can be evaluated by using a priori information about f as follows:
51
To comply with the notation adopted for the result of the first integral, we denote 
as γi:
52
where a change of variables is applied in the second equality: 
and 
. The second integral is then 
, and the left-hand side of Equation (40) can finally be evaluated as 
. With a simple change of subscript from i to j, μj and γj can be readily obtained. The denominator of [Q]ij can thus be evaluated as follows:
53
The integration and calculation of I0(κ) necessary to compute γi and γj are easily performed with available built-in MATLAB functions.
A.2 Numerator of Equation (38)
The numerator of Equation (38) takes the following form:
54
which is an algebraic sum of four expectations with respect to p(x, θ). The processes of transforming these four expectations into computation-friendly forms are similar. We derive the final form for the fourth expectation in detail in Appendix A.2.1 and then provide a brief treatment of the other three expectations.
A.2.1 Fourth Expectation from Equation (54)
We have the following:
55
Again, this is a double integration and can eventually be evaluated by a product of two separate integrals. We will prove this assertion by evaluating the inner and outer integrals, similar to Appendix A.1:
56
As in Appendix A.1.1, we set the following:
57
In this case, ak is defined as follows:
58
to facilitate derivation. Thus, we have the following:
59
which leads to the following:
60
where Dk is dependent on k:
61
Substituting this result into Equation (56), we arrive at the following:
62
The last equality is based on Equation (46), with 
. By taking the natural logarithm of both sides of Equation (62), we obtain the following:
63
where μij, 4 is defined as follows:
64
Because μij, 4 does not depend on θ, Equation (55) can be computed as a product of two integrals, one of which takes the value of 
:
65
Now, the integral over Ωθ can be computed by using the a priori pdf as follows:
66
Here, without a loss of generality, we assume hi > hj > 0, and a change of variables is applied in the second equality: 
and 
. By inserting Equation (66) into Equation (65), we can compute the fourth expectation of Equation (54) as follows:
67
The first, second, and third expectations can be evaluated in a similar manner.
A.2.2 First Expectation from Equation (54)
With similar manipulations, the first expectation in Equation (54) can be evaluated as follows:
68
where γij,1 and μij,1 are respectively defined as follows:
69
and:
70
During the development, hi > hj > 0, and a change of variables is applied in the second equality of Equation (69): 
and 
. The adoption of intermediate equalities such as 
and 
was helpful in arriving at Equation (70).
A.2.3 Second Expectation from Equation (54)
The second expectation from Equation (54) can be evaluated as follows:
71
where 
is defined as follows:
72
With a change of variables, 
and 
, the following can be shown:
73
μij,2 in Equation (71) is defined as follows:
74
During the development, the intermediate variable yk is set exactly as Equation (57), where ak is set as follows:
75
to arrive at Equation (74).
A.2.4 Third Expectation from Equation (54)
The third expectation from Equation (54) can be evaluated as follows:
76
where 
and μij,3 can be easily evaluated by simply swapping the subscripts in Equations (73) and (74), respectively, as follows:
77
and:
78
This completes our development.
APPENDIX B BCRB
We can disassemble the JB term of Equation (27) into two parts (van Trees & Bell, 2007):
79
Specifically, the subscript “D” denotes “data” and “P” denotes “a priori’,” indicating that JD represents the contribution from data and JP represents the contribution from the a priori pdf of the parameter to be estimated. These components are defined as follows:
80
and:
81
from which we can readily identify the term within the outer curly brackets of Equation (80) as the FIM, denoted as JF. Assuming that the first sample is taken at t = 0, the FIM for our current problem can be calculated as follows (Wu et el., 2006):
82
leading to the determination of JD:
83
The JP term can be computed as follows:
84
The final expression arises from the fact that the right-hand side of this equation includes the modified Bessel function of the first kind of order one (Mardia, 1999). By substituting Equations (83) and (84) into Equation (79), we arrive at Equation (29) as the BIM.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.