Design of Sparse Multiband Signal for Precise Positioning With Joint Low-Complexity Time Delay and Carrier Phase Estimation

This paper presents a methodology to design a sparse multiband ranging signal with a large virtual bandwidth, from which time delay and carrier phase are estimated by a low complexity multivariate maximum likelihood (ML) method. In the estimation model for a multipath channel, not all reflected paths are considered, and time delay and carrier phase are estimated in a step-wise manner to further reduce the computational load. By introducing a measure of dependence and a measure of bias for a multipath reflection, we analyse the bias, precision and accuracy of time delay and carrier phase estimation. Since these two indicators are determined by the signal spectrum pattern, they are used to formulate an optimization for signal design. By solving the optimization problem, only a few bands from the available signal spectrum are selected for ranging. Consequently, the designed signal only occupies a small amount of signal spectrum but has a large virtual bandwidth and can thereby still offer a high ranging precision with only a small bias, based on the low-complexity simplified ML method. Numerical and laboratory experiments are carried out to evaluate the ranging performance of the proposed estimation method based on sparsely selected signal bands. Relative positioning, in which we only measure a change in position, based on either the time delay estimates or the carrier phase estimates, is presented as a proof-of-concept for precise positioning. The results show that positioning based on only 7 out of 16 signal bands, sparsely placed in the available spectrum, achieves a decimeter level accuracy when using time delay estimates, and a millimeter level accuracy when using carrier phase estimates. Compared with the case of using all available bands, and without largely decreasing the positioning performance, the computational complexity when using the sparse multiband signal can be reduced by about 80%.

time delay and carrier phase have been widely used to compute position solutions, and are generally estimated for signals from a Global Navigation Satellite System (GNSS) or a ground-based positioning system. The resolution of the time delay estimate is inversely proportional to the signal bandwidth. Therefore, time delay estimation with a high accuracy requires a large signal bandwidth. A nanosecond level accuracy can be attained through an ultra-wideband (UWB) signal [1], [2], or other type of signal that has a large signal bandwidth in the order of 1 GHz.
Recently, signal design has been investigated to improve the performance of time delay estimation when spectral resources are limited. Placing the signal power more toward the frequency band edges, will improve the precision of time delay estimation [3]- [5], according to the Gabor bandwidth (GB) and the Cramer-Rao lower bound (CRLB) analysis. Hence, binary offset carrier (BOC) signals [6], [7] adopted in GNSS, in which the signal power is moved towards the edges of the available signal spectrum, offer better ranging precision than using square pulses. Similarly, a Dirac-rectangular power spectral density (PSD) is proposed in [8]. Since its power is concentrated at the edges of the spectrum, the main-lobe of the auto-correlation function gets narrower, which consequently provides a higher precision in time delay estimation. If signal design is only based on the GB analysis, the ranging signal may become very sensitive to multipath effects. Therefore, a multivariate CRLB, which not only considers the parameter of interest and noise but also the effects of reflections and interference, is considered as the performance criterion for signal design [9], [10].
Given a ranging signal, time delay estimation has been well studied in the past decades. Cross correlation between the received signal and locally generated reference signal [11], which is also referred to as matched filter (MF), is the simplest approach to determine the time delay. But the delay estimate generally contains a bias due to non-resolved reflections in a multipath channel. A multivariate maximum likelihood (ML) method [12], [13] has been developed to jointly estimate the time delay, not only for the line-of-sight (LoS) path, but also for all reflections, however, this requires a tremendous computational resource to obtain the unbiased time delay estimates. With a lot of reflections, the problem may get ill-conditioned and be troublesome to solve properly. Subspace (or eigen-decomposition) methods (e.g., MUSIC, ESPRIT [14]- [16]) have also been applied to delay estimation, which largely reduce the computational complexity compared to the multivariate ML method, and also This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ outperform and provide higher resolution than the simple MF method. Nevertheless, the delay estimate determined by a subspace method will contain a bias in a multipath channel [15]. The subspace method also requires a large amount of snapshot measurements, either from different frequencies or antennas, to compute an accurate sample variance matrix.
For the purpose of positioning, carrier phase estimation has also been investigated, based on different signals. It can provide much more precise ranging information than the time delay estimates due to the small carrier wavelength, but a phase cycle ambiguity has to be dealt with. In GNSS, the carrier phase is obtained from the Q and I outputs of the correlator between the received signal and the replica codes at the prompt branches through the arctangent [17]. Recently, carrier phase estimation on terrestrial signal-of-opportunities (SOPs) has also been studied for ranging. In [18], Yang et al. proposed to simply track the phase directly from a continuous middle sub-carrier of a DVB-T signal, but it is sensitive to multipath. In [19], Khalife et al. first estimated the Doppler frequency offset, and then by integrating it over the observation period, the carrier phase can be obtained for positioning.
Despite different approaches which have been proposed for time delay and carrier phase estimation, unbiasedness cannot always be guaranteed. The CRLB, which only represents the best performance of any unbiased estimator, has however been widely used to guide signal design. To fully explore the benefit of using such a designed signal for ranging, one needs to use the multivariate ML method, which yields unbiased estimates. Considering not only the LoS path but also all reflections in the model, results in significant computational complexity, and becomes less attractive in practice. This raises the question of how to design a ranging signal which enables time delay and carrier phase estimation, and satisfies the requirement on both precision and bias, when not all reflections are considered in the model in order to reduce the computational complexity.
In this paper, we analyse how the spectrum of the ranging signal and multipath impact time delay estimation and carrier phase estimation. The time delay and carrier phase will be jointly estimated based on the maximum likelihood method only for the LoS path and a few reflections, as a compromise between the computational complexity and the overall performance. A measure of dependence and a measure of bias are defined as indicators of the overall ranging performance, which can be used to guide the design of a ranging signal.
Given multiple available signal bands for ranging, we formulate a convex optimization problem based on the ranging performance indicators, to sparsely select only a few signal bands, so that the designed signal only occupies a small amount of frequency resource. Thereby only a correspondingly small amount of measurements in frequency domain are used for parameter estimation, which further reduces the computational load, while the user specified performance can still be achieved through a simplified model. A multiband orthogonal frequency division multiplexing (OFDM) signal, which has been adopted for UWB communication [20], is considered as an example in this paper. The remainder of the paper is organized as follows. A measurement model of the sampled channel frequency response based on a multiband signal and the estimation problem are presented in Section II. Statistical analysis on time delay estimation is provided in Section III for both the full model that considers both the LoS path and all reflections, and the simplified model with the LoS path and only a few reflections. Then, in Section IV, estimation of the complex gain and the associated carrier phase are presented both with the full and simplified model. Afterwards, an optimization-based signal band selection is proposed in Section V. Section VI characterises the performance of time delay estimation and carrier phase estimation using the full model and the simplified model, and presents an example to design a sparse multiband signal for positioning. The experimental setup and results are shown in Section VII. Finally, conclusions are drawn in Section VIII.
Notation: E{·} denotes the expectation operation. R (C) indicates a real (complex) variable. (·) denotes a random variable. (·) denotes an estimate. (·) T , (·) * and (·) H denote transposition, conjugate and Hermitian operation, respectively. | · | denotes the absolute or modulus value. Uppercase boldface letters (e.g., F ) are used for matrices, and [·] ij denotes the element in the i-th row and the j-th column of a matrix. Lowercase boldface letters (e.g., r) are used for column vectors. A normal lowercase letter or [·] i is used to denote an element in a column vector. I N is the identity matrix of size N . 1 N denotes an N -by-1 vector of one. diag(·) denotes a diagonal matrix formed by its vector argument. tr{·} denotes the trace of a matrix. {·} and {·} denote the real part and imaginary part of a complex value, respectively. x 2 = x T x denotes the Euclidean norm. x 1 = |x 1 | + . . . + |x N | denotes the sum-absolute-value (i.e., l 1 norm).

II. SIGNAL MODEL
Here we consider a multiband OFDM signal as an example, but in practice any modulation can be used in arbitrarily located signal bands. As shown in Fig. 1, there are M available bands in the allocated signal spectrum that can be used for positioning or communication. Given a medium total signal bandwidth (e.g., 100-200 MHz), all signal bands can be received simultaneously through a single RF front-end and analog-to-digital converter (ADC). Additionally, we define the bandwidth between the two signal bands at the band edges as the virtual signal bandwidth (as shown in Fig. 1), no matter how many activated bands there are in between.
After frame synchronization and fast Fourier transformation (FFT), the channel frequency response of the subcarriers across the different signal bands can be obtained based on a training symbol. The response not only contains information on the time delay but also on the carrier phase. Details on modulation and demodulation of OFDM signals in each signal band are omitted here and can be found in [21]. For communication, in order to properly demodulate the signal, a simple frequency domain equalization is used to compensate for the distortion of the entire multipath channel, based on the channel frequency response obtained from channel estimation. However, for the purpose of positioning, our aim is to estimate the time delay and the carrier phase specifically of the LoS path in a multipath channel from the sampled channel frequency response.
Given an L-path channel, the sampled complex baseband channel impulse response is given by where T s is the sample interval, τ l denotes the propagation delay of the l-th path, x l denotes the complex propagation gain, and α l denotes the modulus of x l and is a real value. In addition, f c denotes the central carrier frequency of the multiband signal. In this paper, carrier phase tracking explicitly refers to phase tracking of the central carrier frequency f c .
In principle, M signal bands could be used for ranging and positioning, however, to improve spectral efficiency and to decrease computational complexity, we can activate M a (M a ≤ M ) signal bands for this task. The selection of the activated signal bands will be introduced in Section V. In the present explanation, each OFDM signal band is assumed to have the same bandwidth and each contains N subcarriers, but this can be different in practice. In addition, the channel frequency response is assumed to satisfy the complex multivariate Gaussian probability density function (PDF). The measurement model for the sampled channel frequency response is given by where F denotes FFT, H denotes the channel estimates at different subcarriers in different signal bands (N subcarriers in each of the M a signal bands), f i = iΔf denotes the subcarrier frequency of the i-th subcarrier (−N/2 ≤ i ≤ N/2 − 1), Δf denotes the subcarrier spacing, f m denotes the centre frequency of the m-th activated signal band with m = 1, . . . , M a . The measurement noise is assumed to be white Gaussian with a variance of σ 2 which is known or can be estimated a priori. In addition, we assume that the unknown propagation delay τ and the unknown complex gain x are from two completely disjoint sets.
Since the unknown parameters are present not only in x but also in A(τ ), we jointly estimate the time delay and complex gain from the channel estimates based on the Maximum Likelihood (ML) method through the following minimization Although there are unknown parameters in A(τ ), its dimension is known a priori. In our case, the number of rows of the design matrix A(τ ) is determined by the number of the subcarriers in each signal band and the number of signal bands. The number of paths (i.e., the number of columns in A(τ )) should also be specified. The impact of selecting different number of paths for parameter estimation will be discussed in Section III and IV. It has been proven in [22] that ifτ andx are the global minimizers of (4),x must satisfŷ when the variance matrix Q H is a diagonal matrix with identical elements on the main diagonal, as defined in (2). Therefore, using the solutionx to rewrite the cost function (4), we can first estimate the propagation delay, then compute the complex gain and its corresponding carrier phase.

III. TIME DELAY ESTIMATION
Time delay estimation is the first step to obtain the solution from the ML-based cost function (4), where one should determine the number of propagation delays (i.e., the number of paths) that should be jointly estimated.
First, all reflections are considered in the design matrix assuming that the number of paths is known a priori, which is referred to as the full model. We define a measure of dependence between the LoS component and a reflection for delay estimation, and analyse how this dependence influences the accuracy (i.e., the variance) of the LoS delay estimator.
A low-complexity simplified model is also proposed in this section, in which not all reflections are considered in the design matrix. Although we may determine the number of paths in a multipath channel through model order estimation techniques, such as minimum description length (MDL) [23] and generalized Akaike information criterion (GAIC) [24], [25], they generally require a large number of data snapshots and may not provide the exact number of paths. The simplified model seems more practical than the full model to implement in practice. As less unknown parameters are estimated, the computational efficiency and possibly the precision can be improved, but the resulting delay estimator for the LoS path likely will be biased. Therefore, we define a measure of bias, and analyse how an unconsidered reflection in the simplified model impacts delay estimation, with the goal of verifying that the simplified model will eventually meet the requirements.

A. Full Model
In this subsection, we analyse how the accuracy of the delay estimator for the LoS path deteriorates when a reflection is considered in the full model for time delay estimation in an attempt to achieve unbiasedness.
Since the design matrix is partially unknown, combined with (5), and considering white Gaussian noise, the delay estimates can be equivalently derived from the minimization of the following nonlinear cost function [12], [13], [22] where the variance matrix Q H = σ 2 I NM a , as defined in (2), and the complementary projection matrix is defined by Therefore, the propagation delay can be equivalently estimated through the following cost function, with the projection matrix P A(τ ) defined in (7), In order to find the maximizer of this non-linear cost function, there are generally two different approaches [26]: direct search methods and gradient methods. The direct search methods do not require any evaluation of derivatives of the non-linear cost function, but it can be time consuming to find the solution for the cost function when it contains multiple variables. The gradient methods require derivatives to determine the direction of each iteration. For a cost function, like (6), with a projector, the reader can refer to [22]. However, the additional approximation effect, introduced by ignoring higher order derivatives of the latter, should be properly diagnosed and associated biases should be carefully dealt with (see e.g., [27], [28]). More details on solving a non-linear problem based on gradient methods can be found in [29], which extensively discussed different iterative gradient methods.
Here we consider the direct-search approach to obtain the time delay estimates. Given a multivariate cost function, a significant computation time is required to obtain the optimal solution. One may consider using the alternating projection (AP) [30] to iteratively compute the solution, which is computationally attractive for solving multivariate non-linear MLE, and also provides a better performance than the sub-optimal techniques (e.g., MUSIC [31]).
To simplify the notation and derive a closed-form expression, a two-path channel is considered here as an example to analyse how a reflection deteriorates the accuracy. The channel can contain more paths in practice, but this will not change the conclusion of this subsection. Considering white Gaussian noise, the uncertainty of unbiased ML-estimation can be derived from the CRLB [32], which is the inverse of the Fisher information matrix (FIM) F (τ ).
The propagation delay estimates and the complex gain estimates are assumed to be uncorrelated. Also, to derive the FIM for the propagation delay, the complex gain is assumed to be known a priori. In [33], the FIM for the propagation delay in a two path channel has been derived as follows, where α 2,1 = α 2 /α 1 and τ 2,1 = τ 2 − τ 1 denote the relative propagation gain and delay, respectively, f contains the frequencies for each subcarrier in each activated signal band f i + f m , denotes the element-wise dot product, and q(τ 2,1 ) = {a(τ 2,1 )} = cos(2πf τ 2,1 ).
Based on (9), the variance of the unbiased delay estimator σ 2 τ 1 for the LoS path when the full model is used, is given by In addition, for comparison, the variance of the unbiased delay estimator derived for a single path channel is given by It presents the accuracy when the propagation channel only contains the LoS path, or the precision of time delay estimation when there are more paths, but only the LoS path is considered in the model. As q(τ 2,1 ) T f f ≤ f T f , we have To simplify notations, we first define and the real part of its second derivative is given by Therefore, the variance of the unbiased delay estimator in a two path channel is rewritten as By we define ς(τ 2,1 ) as a measure of dependence for delay estimation, which indicates the dependence level between the LoS component and the reflected component. Throughout this paper, a reflection with a non-zero measure of dependence for delay estimation is referred to as a dependent reflection for delay estimation. Comparing (14) with (11), the accuracy of the estimator for a two-path channel and for a single-path channel are identical if ς(τ 2,1 ) is zero. By using the full model for time delay estimation in a practical multipath channel, the accuracy of the delay estimator for the LoS path will become poorer when more dependent reflections need to be considered in the design matrix A(τ ). It should be mentioned that considering more reflections in the model requires a considerably large computation time to obtain the optimal solutions.
On the other hand, one may not consider all reflections in the design matrix A(τ ), or in the extreme case no reflection is considered, which is referred to as the simplified model for time delay estimation. Based on the simplified model, the computational burden can be largely relaxed, and the precision of the delay estimator can be also improved as shown in (12). However, the estimator becomes biased. The resulting bias will be analysed in the following subsection.

B. Simplified Model
Using the simplified model, in which not all reflections are taken into account to determine the time delay, requires less computational time and can provide a higher precision than when using the cost function based on the full model (8), but the resulting delay estimate of the LoS path could include a bias. Again, in order to provide a closed-form expression and to keep the derivations mathematically manageable, we consider a two path channel to analyse the bias in the estimated delay with the simplified model. Hence, the measurement model (2) is changed to In addition, only for the purpose of bias analysis in delay estimation, the complex gains are assumed to be known. The propagation delay based on the ML method in the simplified model can be determined by the following minimizatioň where only one path is considered, i.e., in a two-path channel scenario. However, it is still difficult to obtain a closed-form expression for the bias based on the non-linear cost function (17). Therefore, only for the purpose of bias analysis, we linearise (17) through Taylor expansion at τ 1 which is the true time delay of the LoS path, and ignore second and higher order terms. The linearised cost function is given by where Q H = σ 2 I NM a as defined in (2). Consequently, the linear model is given by andã As higher order remainders in the Taylor expansion have been neglected in the linearisation, the estimator becomes biased even when the measurements are unbiased [27], [28]. Given the complex non-linear design matrix defined in (3), the second order remainder will contribute to the bias in the imaginary part due to non-linearity, and the third order remainder will cause a bias in the real part. Typically, the higher than second order remainders are very small and can be ignored, and the bias introduced by the second order remainder due to non-linearity is jointly determined by the signal structure (i.e., the design matrix a(τ ) itself) and the quality of the estimator [28]. Given a reasonable signal-to-noise ratio (SNR), for example larger than −10 dB, the bias will be relatively small and will not dominate the accuracy of the estimation.
A bias τ b introduced by not considering the reflection in (16) is determined by where ϕ denotes the phase from x H 1 x 2 . The estimator obtained from a complex estimation problem is complex, and the imaginary part is introduced mainly due to the non-linearity of a(τ ), which is relatively small. As the time delay estimator should be a real number, we simply take the real part of the estimator to indicate the delay bias. The justification of the approximated bias is shown in Section VI.
The actual bias of delay estimation also depends on whether the reflection is constructive or destructive to the LoS component, as well as the non-linearity of a(τ ). Without loss of generality, we analyse the envelope of the bias, or the maximum delay bias. Combined with where p(τ 2,1 ) = sin(2πf τ 2,1 ), the maximum of the delay bias τ b as a function of the relative propagation delay τ 2,1 is given by Afterwards, we define the measure of bias as The measure of bias (τ 2,1 ) = 0 means that the reflection will not cause a bias in the estimateτ 1 , even if this reflection is not considered in the estimation model. As we can see, the measure of bias depends on the signal pattern through s (τ 2,1 ), but also on the relative gain α 2,1 . Therefore, a reflected signal component, which is largely attenuated compared to the LoS component, will only cause a small bias when it is not considered in the simplified model. If the channel contains more than two paths (i.e., L ≥ 2), the bias due to the unconsidered reflections can be derived as Hence, when more reflections are neglected in the model, the resulting bias in the worst case will be the superposition of the biases derived from multiple two-path channels with the fixed LoS path, in which the reflections are treated separately. It should also be mentioned that if a(τ 1 ) and a(τ 2 ) are orthogonal (i.e., s(τ 2,1 ) = 0), though likely not realistic in practice, the multivariate cost function based on the maximum likelihood method (8) can equivalently be simplified to a single variate function (i.e., 1D MLE), which is also equivalent to the MUSIC algorithm when the noise space U n can be accurately determined from the sample covariance matrix.

C. Flop Count
A numerical operation (e.g., addition, multiplication, squareroot, etc.) can be defined as a flop, and the number of required flops can be used to evaluate the computational complexity [34]. As the direct search method is applied here to determine the time delay, the cost function (8) should be computed for each value in a search grid that contains all possible delay values. Here we only derive the required number of flops for a single grid-point, and the size of the search grid is not considered.
Considering M a signal bands, N subcarriers per band, and L paths in the estimation model, then the number of required flops is derived as follows, For more details, the reader is referred to [34], as well as Appendix A. The required flops in (26) is dominated by the term (1 + 4 L)(M a N ) 2 . As the number of measurements M a N is generally larger than the number of paths L considered in the model to avoid rank-deficiency, considering less paths can largely reduce the computational complexity. Given a fixed number of paths L, using less signal bands M a can also reduce the computational complexity considerably. Therefore, we aim to design a ranging signal which provides a good balance between computational complexity and ranging performance. This signal design will be introduced in Section V.

IV. CARRIER PHASE ESTIMATION
Once the propagation time delays are determined by MLE, we can continue to estimate the complex gains based on (5). In a similar way, as discussed in the previous section, one may use the simplified model that only contains the LoS path in A(τ ), to estimate the complex gain of the LoS path and its corresponding carrier phase, if a biased solution for the carrier phase is acceptable. Otherwise, one needs to construct a full model based on the delays from the different paths and jointly estimate their corresponding complex gains and carrier phases. Similarly, we will analyse the accuracy of complex gain estimation based on the full model and the bias introduced in the simplified model.

A. Full Model
In this subsection, we analyse the accuracy of the unbiased gain estimator. Again, a 2-path channel (16) is considered as an example here. To construct the full model, the design matrix A(τ ) in (2) contains both the LoS component a(τ 1 ) and a reflected component a(τ 2 ), and they are constructed based on the unbiased propagation delay estimates from Section III A. To simplify notation in the following derivations, a(τ 1 ) and a(τ 2 ) are replaced by a 1 and a 2 , respectively. It should be mentioned that when there are more reflections in a multipath channel, a 2 will be extended from a vector into a matrix that contains all reflected components.
Although the design matrix A(τ ) should contain all possible paths in a channel, and the complex gain can be estimated for all paths, only the complex gain of the LoS path is of interest for ranging and positioning. The carrier phase of the LoS pathΦ can be derived from the complex gainx 1 . Based on the partitioned model (16) and the MLE solution shown in (5), the complex gain of the LoS path is given bŷ and the variance of the complex gain x 1 is computed by If the channel only contains a single path, the variance is given by Now we analyse how the accuracy changes when a reflection is taken into consideration in the model for complex gain estimation. The variance of gain estimation in a 2-path channel (29) is rewritten as where sin −2 (ϑ) is a scalar in a 2-path channel. The angle ϑ measures the degree of dependence [35] between the LoS component a 1 and a reflection a 2 in A(τ ), which is written as follows Combined with (13), the variance σ 2 x 1 is inversely deteriorated by When |s(τ 2,1 )| approaches zero, the LoS component a 1 and the reflection a 2 become linearly independent in A(τ ), and the variance σ 2 x 1 will be the same as σ 2 x 1 . Hence, |s(τ 2,1 )| is defined as the measure of dependence for gain estimation in this paper, and a reflection with a non-zero measure of dependence for gain estimation is defined as a dependent reflection in gain estimation. Equivalently, if the angle ϑ as shown in Fig. 2 equals ±π/2, a 1 and a 2 are orthogonal, and if ϑ is zero, a 1 and a 2 are fully dependent.
Considering both the LoS path and a reflection in A(τ ), the variance of the complex gain of the LoS component (31) can be rewritten by As we can see, the accuracy of the estimator for the LoS path gets poor if a dependent path is added to the model, because |s(τ 2,1 )| and the corresponding sin −2 (ϑ) are large.

B. Simplified Model
Instead of jointly estimating the complex gain for all paths in a multipath channel, which could lead to a poor accuracy when dependent columns are involved in the design matrix A(τ ), one can estimate the complex gain for a few paths, i.e., A(τ ) contains only a few reflected components, with better precision and lower computational load, but at the cost of a bias. Since the design matrix A(τ ) is not fully reconstructed, it is referred to as the simplified model. We analyse the bias of the complex gain and its corresponding phase using the simplified model in this subsection.
The design matrix A(τ ) is constructed using the time delay estimates from Section III. These time delays could be biased due to unconsidered reflections in the simplified model as presented in Section III-B, and thus the computational complexity could be lower there. Without loss of generality, we consider a biased delay estimate to construct the design matrix, and also analyse how the time delay bias impacts carrier phase estimation using the simplified model.
In order to derive a closed-form expression for the resulting bias, we again use a simple 2-path channel (16) here, and the same stochastic model as shown in (2). The simplified model for gain estimation is constructed aš Then the complex gain can be derived by As a special case, if the propagation time delay is unbiasedly estimated through a full model in (8), τ b will be zero. Consequently, the precision of the gain estimator can be written as Sinceǎ 1 is a complex vector which has a similar structure as the one defined in (2), the precision ofx 1 will only be determined by the signal structure itself, instead of the biased delay estimator. With (36) and (16), the complex gain of the LoS path is given byx Once the complex gain of the LoS component is determined, the corresponding carrier phase can be obtained from its arc-tangent. However, the carrier phase estimate could become biased, because of the delay bias τ b , and also because of the unconsidered dependent reflection s(τ 2,1 + τ b )x 2 . The bias in the carrier phase varies with the number of paths and their propagation delays and gains. Then, the carrier phase derived from the simplified model is given byΦ To analyse how the bias of the time delay and an unconsidered reflection influence the carrier phase estimator, the maximum absolute phase bias is derived in a 2-path channel. The geometric interpretation is shown in Fig. 3. The maximum phase bias is given by where |Φ b1 | max denotes the maximum bias introduced by the biased delay estimate, and |Φ b2 | max denotes the maximum phase bias introduced by an unconsidered reflection. Except for extremely close-in multipath, reflected components in a multipath channel are weaker than the LoS component (i.e., α 2,1 < 1). Hence, the phase bias due to an unconsidered reflection Φ b2 can be approximated by the first order term of its Taylor expansion, As |s(τ b )| < 1 when τ b = 0, the time delay bias will enlarge the phase bias in the simplified model. Therefore, the ranging signal, which can improve the performance of time delay estimation and offer an accurate time delay estimator (i.e., keep τ b small), will consequently reduce the bias in estimating the carrier phase. The numerical results on this can be found in Section VI.

C. Flops Count
The number of required flops is computed as an indication of computational complexity for complex gain estimation, from which the carrier phase of the LoS path is obtained. When the propagation time delays for both the LoS path and the reflections have been estimated through the multivariate ML method in Section III, the number of required flops are given by where L denotes the number of path considered in the model, M a denotes the number of activated signal bands, each containing N subcarriers. The derivations can be found in Appendix A. The computational complexity is dominated by 4L(M a N ) 2 . As M a N is generally much larger than the number of the considered path L, using the simplified model that considers less paths will largely reduce the computational complexity.
It should be noted that as the design matrix for carrier phase estimation is constructed based on the time delay estimates, the actual computational complexity to compute the carrier phase should also include the one for time delay estimation (cf. Section III-C.).

V. SIGNAL DESIGN FOR PRECISE POSITIONING
In this section, we aim to design a sparse multiband ranging signal, which uses limited spectrum resources (e.g., M a signal bands out of M ), reduces the computational complexity, and in a multipath channel still achieves a user-specified precision for time delay estimation and keeps the bias small. As we determine the carrier phase based on the time delay estimate(s), in return, the bias in carrier phase estimation will also remain small. Sparsely selecting multiple signal bands for ranging and positioning is similar to the problem of sensor selection [36]- [38]. In [39], Li formulated a convex optimization problem to sparsely select frequency tones for indoor ranging, in which the CRLB of the propagation time delay for a 2-path channel is used as a criterion. However, the number of reflected components is generally larger than two in a multipath channel, and it is impractical to resolve the propagation delays for all paths. Thus, neither the accuracy of delay estimation, nor the robustness against multipath can be guaranteed in the designed sub-optimal signal.
In this paper, instead of using the CRLB derived for a 2-path channel as a criterion, we employ constraints on the precision of time delay estimation (i.e., the CRLB derived for a single path channel), on the measure of dependence for delay estimation, and on the measure of bias, to formulate an optimization problem. As discussed in Section III, the measure of dependence for delay estimation indicates how the accuracy deteriorates when a reflection is considered in the model, and the measure of bias indicates how large the bias is if we do not consider such a reflection in the simplified model.
To simplify the discussion, the signal power in each signal band is assumed to be the same, i.e., the more signal bands are used for positioning, the larger the total signal power becomes. We consider M different signal bands and introduce a binary selection vector for these bands, (43) where w m = 1(0) indicates that the signal from the m-th signal band is activated (muted), and used (not used) for positioning. The binary selection vector determines the signal pattern for ranging.
We first give a constraint on the precision of the time delay estimator. Similar to (11), but considering the frequency relation of the different signal bands and the binary selection vector, the variance of the delay estimator needs to be smaller than a user-specified threshold σ 2 τ . Hence, the constraint is given by where the frequency f m + f i , the measurement noise σ 2 and the propagation gain α 1 (or SNR) are assumed to be known to the user. Then, to consider the impact of multipath on delay estimation, we put constraints on the measure of dependence and on the measure of bias. Without obtaining the channel information a priori, we simply consider a set of reflections that cover a certain range (e.g., with a relative propagation distance from 0.6 m to 15 m). The relative delay is derived from the relative propagation distance, and the relative propagation gains are derived from the free-space path-loss (FSPL) model.
As shown in (24), a reflection close to the LoS path may have a strong relative signal power and is likely to cause a large bias when such a reflection is not considered in the model. Therefore, some reflections with strong signal power can be considered in the simplified model to mitigate the bias. One can place a constraint on the measure of dependence for delay estimation, so that the precision of the LoS estimator will not deteriorate when such a reflection is considered for delay estimation. Combined with the binary selection vector w, the measure of dependence for delay estimation (15), which is required to be smaller than a user-specified threshold c ς , is rewritten as and the user-specified set U I contains the to-be-considered reflections in the delay estimation model. Other weak reflections could be neglected in the simplified model, so that the computational complexity will not be increased significantly. Therefore, we place a constraint on the measure of bias, so that the bias of delay estimation still remains small, even though those reflections are not considered in the model. Similar to (23), the measure of bias with a binary selection vector w, which is required to be smaller than the user-specified threshold c , is given by where The user-specified set U II contains the reflections that would not need to be considered in the delay estimation model. However, due to the binary nature of the selection vector w, convexity cannot be guaranteed. Thus, we relax the binary selection vector w to g with inequalities (i.e., 0 ≤ g i ≤ 1), so that each signal band is activated with a weight.
Given M available signal bands with the frequency vector f and the signal condition (i.e., α 1 and σ 2 ), the optimization problem to sparsely select multiple signal bands for ranging, with the user-specified thresholds (i.e., σ 2 τ , c ς and c ) and the user-specified set (i.e., U I and U II ), can be formulated as follows, minimize ||g|| 1 subject to in which the l 1 norm is used in the objective function to produce a sparse selection vector. It should be mentioned that the required precision for the propagation delay, the measure of dependence for delay estimation, and the measure of bias should not be exceeded the ones when using all M available signal bands. To retrieve a binary selection vector w, we activate the signal band for positioning once the corresponding selection coefficient in g is larger than the user specified threshold c t , An example of the design of a sparse multiband signal will be provided in Section VI-C.

VI. NUMERICAL RESULTS
This section presents the numerical results on the performance analysis of time delay estimation and carrier phase estimation, based on the full model and the simplified model. This is followed by a discussion of the computational complexity in relation to the number of considered paths in the estimation model and the number of signal bands. Finally, an example of designing a sparse multiband signal for ranging is shown.

A. Full Model and Simplified Model
In this subsection, we analyse the performance of time delay estimation using both the full model and the simplified model, as well as the performance of estimating the complex gain and its corresponding carrier phase.
The estimation performance is not only determined by the number of paths considered in the model, but also by the signal pattern. Here, we consider the following signal patterns to evaluate the performance of time delay and carrier phase estimation: two edge signal bands (i.e., only the 1-st and 16-th band shown in Fig. 1), a sparse multi-band signal (consisting of the {1, 2, 5, 6, 12, 13, 16}-th signal band), all M = 16 signal bands, and 7 contiguous signal bands which provides less virtual signal bandwidth than the other three patterns. Each signal band has a bandwidth of 10 MHz and N = 64 subcarriers.
To analyse how a reflection impacts the estimation performance, a simple two path channel is considered. The relative distance between the LoS path and the reflection varies from 0 m to 15 m. The relative gain of the reflection is determined based on the FSPL model [40].
First, the delay bias (21) introduced in the simplified model, which is derived based on linearisation, is analysed. To derive a closed-form expression for the measure of bias, the Taylor expansion is applied to linearise the non-linear design matrix in (17). As stated in Section III B, apart from the unconsidered reflection and the non-linearity of the problem, an additional bias is introduced by the effect of the linear approximation on the design matrix a(τ ). Fig. 4(a) shows the bias derived from the  (17), (b) measure of dependence for delay estimation ς(τ 2,1 ) (cf. (15)) and measure of bias (τ 2,1 ) (cf. (24)).   (17), only when a reflection has a large power and is close to the LoS path. Nevertheless, the closed-form bias generally indicates how a reflection impacts the delay estimator, and can be used to define the measure of bias in (24). Fig. 4(b) shows the measure of dependence and the measure of bias. A reflection with a zero measure of dependence does not correspond to a zero measure of bias. Therefore, an independent but unconsidered reflection can still cause a delay bias in the simplified model. Compared to the measure of dependence and the measure of bias derived based on 7 contiguous bands, the indicators derived for all M = 16 bands become smaller due to a larger virtual signal bandwidth, especially for close-in reflections.
Then, Fig. 5 presents the accuracy of time delay estimation for the aforementioned signal patterns. The performance of an unbiased estimator in the full model is evaluated by the CRLB. The accuracy of a biased estimator in the simplified model is quantified by the mean-square-error (MSE), which can be decomposed in a variance-plus-bias-square. Given a fixed virtual signal bandwidth (i.e., 160 MHz), using more signal bands in the full model improves the accuracy of the delay estimator. Although the overall improvement is relatively limited (e.g., a few centimetres) as shown in Fig. 5(a), using more signal bands can increase the robustness against multipath, as the root-CRLB is less affected by the reflections. In addition, as shown in Fig. 5(b), using the simplified model, the root-MSE is dominated by the bias. Similarly, using more signal bands can generally reduce the bias resulting from not considering the reflection in the model. Also in Fig. 5, one can notice that when a reflection with a very small relative distance (e.g., less than a few decimetres), the simplified model can outperform the full model. A very close-in multipath normally complies with a large measure of dependence. Considering it in the full model leads to a poorly estimated LoS path time delay. On the other hand, close-in multipath only leads to a relatively small measure of bias as shown in Fig. 4(b). Therefore, using the simplified model, only a small bias will be introduced to the time delay estimator, and the precision will not be deteriorated by the close-in reflection.
Next, the performance of carrier phase estimation is presented. The same signal patterns as mentioned above are used in the following analysis. The complex gain and the corresponding carrier phase for the LoS path is estimated based on a reconstructed design matrix, in which the time delay estimates are initially assumed to be unbiased (i.e., τ b = 0). Using the full model (27), Fig. 6(a) shows the corresponding scaling factor sin −1 (ϑ) in the standard deviation of the LoS gain estimator, which is derived from (31). When the measure of dependence |s(τ 2,1 )| = 1, the LoS component a 1 and the reflection component a 2 will be fully dependent, and sin(ϑ) = 0 (see (33)), then the complex gain will be poorly estimated. Using a larger virtual signal bandwidth can reduce the measure of dependence for close-in reflections. However, with less signal bands (e.g., using only two edge signal bands), the accuracy of complex gain estimation becomes sensitive to the non-close-in reflection.
Alternatively, the simplified model can be applied to estimate the carrier phase, and consequently it becomes biased. Fig. 6(b) shows the maximum carrier phase bias |Φ b2 | max in (40) for the simplified model with an unbiased time delay estimate. Given a fixed virtual signal bandwidth, using more signal bands will improve the robustness against multipath, as less bias will be introduced when using the simplified model. To achieve a certain estimation performance, not all signal bands are needed. With 7 contiguous bands which cover a smaller virtual signal bandwidth, the resulting bias will be dominated by the close-in reflections due to a reduced time resolution. In addition, as shown in Fig. 4 (dashed lines), due to the large measure of dependence and the large measure of bias, the unbiasedness of the delay estimator becomes difficult to achieve.
The low complexity simplified model can also be applied for time delay estimation (Section III-B), and the design matrix for complex gain estimation will be reconstructed using a biased time delay estimate. Consequently, the biased delay estimate of the LoS path will introduce an extra bias Φ b1 , as shown in Fig. 6(c,d), but the phase bias Φ b is dominated by the unconsidered multipath (Section IV-C) Φ b2 . Although the integer phase cycle ambiguity should be correctly estimated in order to retrieve the geometric information, the ultimate bias is small in distance and likely acceptable to the user, even if only a few bands are used within the virtual signal bandwidth. In such a condition, positioning based on the carrier phase will still largely outperform the one based on the propagation time delay, which will be demonstrated in Section VII-C.
From Fig. 4(b) and Fig. 6, one can notice that the impact of the signal pattern on time delay estimation and carrier phase estimation is similar. A large virtual signal bandwidth helps to distinguish relatively close-in multipath, and using more signal bands improves the overall robustness against multipath (i.e., less bias in the simplified model, higher precision in the full model). Moreover, not all signal bands are needed to achieve a certain ranging performance.

B. Computational Complexity
In this subsection, the computational complexity for time delay and complex gain estimation is evaluated by the number of required flops. Here, the number of paths L considered in the estimation model is set to be 1, 2 and 5, and the number of signal bands M a is varied from 6 to 16. Fig. 7(a) shows the number of required flops for time delay estimation. As the direct search method is used, Fig. 7(a) only presents the number of flops required for a single element in the search grid. Fig. 7(b) shows the number of required flops for complex gain estimation.
Given a fixed amount of signal bands, considering one path less in the time delay estimation model reduces the amount of required flops by about 49%. On the other hand, considering a fixed amount of paths in the model, using one less signal band reduces the number of the required flops by about 27%. A similar behaviour is also observed for carrier phase estimation.
In practice, as a compromise between the computational efficiency and the overall performance, we can jointly estimate the time delays only for a few strong reflections (i.e., keep L limited). A reflection with a larger relative delay, generally, is associated with a smaller relative gain. Consequently, the resulting bias will be small when such a reflection is not considered in the simplified model. Not all bands are needed to achieve the user specified thresholds for the measure of dependence and the measure of bias. Therefore, one can use less spectral resources and samples (i.e., keep M a small) to even further reduce the computational complexity.

C. Example of Sparse Multiband Signal
It is assumed that there are M = 16 contiguous signal bands potentially available for positioning, and each signal band has the same signal power and the same bandwidth of 10 MHz.
In order to provide numerical values for the constraints in (47), we first analyse the measure of dependence for delay estimation ς(τ 2,1 ) and the measure of bias (τ 2,1 ) as a function of the relative propagation distance when M = 16 bands are activated for positioning. These results were presented in Fig. 4(b). There will be a strong dependence between the LoS path and the reflections with a relative propagation distance of less than 0.6 m (i.e., ς(τ 2,1 ) ≥ 0.7). Therefore, it is unfeasible to further decrease the measure of dependence based on the existing signal bandwidth to improve robustness against multipath, other than using even a larger virtual signal bandwidth. Hence, given 16 signal bands each with 10 MHz bandwidth, we only consider reflections that are 0.6 m or further away from the LoS path.
First, the standard deviation of the time based ranging error is fixed to 0.003 m, which only indicates the lower bound for delay estimation in an ideal situation where the channel only contains a LoS path. Then, we set c ς to 0.7 as the maximum threshold for the measure of dependence when the relative propagation distance ranges from 0.6 m to 2.3 m in set U I , so that the precision will not decrease significantly when a reflection is taken into consideration. As shown with (14), the variance is doubled when the measure of dependence ς(τ 2,1 ) is 0.7. Finally, we may simply neglect reflections that have a relative propagation distance from 2.3 m to 15 m in set U II , and expect that the estimation bias for the LoS will remain small. Hence, we set threshold c for the bias to 0.1 m, assuming we can accept a 0.1 m bias in time-based ranging.
The solution of the optimization problem (47) is obtained by the CVX toolbox [41]. Fig. 8(a) shows the relaxed selection vector g, and its binary selection vector w when c t = 0.5. In order to achieve a decimetre level time-based ranging accuracy Based on the activated signal bands, as shown in Fig. 8(a), Fig. 8(b) presents the value of the measure of dependence ς(τ ) and the measure of bias (τ ), cf. (15) and (24). Using a few sparse signal bands, the measure of dependence becomes larger than the one using all 16 signal bands, as shown in Fig. 4(b). However, as indicated by the measure of bias, a decimeter level bias will be introduced in the time delay estimate, if the reflection is not considered in the simplified model. Using the designed signal, the user can apply the simplified model for delay estimation, in which not all reflections are considered.
The carrier phase can also be estimated with the simplified model for precise positioning. The measure of dependence for complex gain estimation based on the designed signal is also presented in Fig. 8(b) by the blue dotted line, which links to the phase bias in the simplified model as shown in (41).
The performance of parameter estimation using the designed sparse multiband signal can also be found in Fig. 5 and Fig. 6, in which the so-called 'sparse multiband signal' is the one shown in Fig. 8(a).

VII. EXPERIMENTAL RESULTS
In this section, we provide experimental laboratory results to demonstrate the performance of propagation time delay and carrier phase estimation based on the proposed sparse ranging signal. Relative positioning is applied in this paper as a proof-of-concept for demonstration of precise positioning. In practice, the system can be expanded with more transmitters, and other positioning techniques can be used to obtain also absolute position solutions.

A. System Setup
The experiment is based on the platform we recently developed in [42], in which Ettus X310 USRPs (Universal Software Radio Peripherals) are used as radio transmitter/receiver with an effective bandwidth of 160 MHz and a sample rate of 200 MSamples/s. Here, we consider a fully synchronized system, where all transmitters and the receiver are synchronized by a 1 pulse-per-second (PPS) and a 10 MHz reference signal from a common clock distributor. Without the need to consider any clock error, three transmitters are used for a 2D positioning system. Fig. 9(a) shows the experimental system set-up in the lab. The receiver antenna is placed on a rail, so that it can be moved forth and back smoothly. The location of the transmitter (Tx) antennas and the height of the receiver (Rx) antenna are determined through a professional land-surveying Total Station by measuring angles and distances within a local coordinate system, in which coordinates typically have an error of less than 10 mm. Besides, a few locations on the rail have also been measured by the Total Station as reference points, which are used to interpolate to any other location on the rail. During the experiment, the receiver is stopped at the reference points for a few seconds.
The performance of the position solutions depends not only on the quality of the measurements but also on the positioning geometry. Given the system setup in the lab, the horizontal dilution of precision (HDoP) [17] is around 1.6, which means that the horizontal position uncertainty will be about 1.6 times larger than the range measurement uncertainty.
The central carrier frequency is fixed at 3.5 GHz and there are 16 signal bands each with a bandwidth of 10 MHz available for ranging and positioning. The signal shown in Fig. 8(a) is the sparse multiband signal used in the experiment, in which 7 out of 16 signal bands are activated for positioning. To provide a comparison, we also evaluate the performance when using all 16 bands.
To avoid inference among the transmitters, a time division multiplexing (TDM) scheme is applied, where each transmitter occupies a time slot of 100 μs per 1 ms as shown in Fig. 9(b), and each transmitter repeats its signal every 120 ms.

B. Time Delay and Carrier Phase Estimation
Based on the sparse multiband signal as designed in Section VI-C, the time delay and the carrier phase are estimated only based on the simplified model.
First, two versions of the simplified model are considered for time delay estimation. One contains only a single path in the model, which is referred to as 1D MLE. The other model contains one additional reflected path, and is referred to as 2D MLE. The alternating projection [30] is applied as a realization of the multivariate ML method to jointly estimate the delay for the LoS path and the reflection.
For the signal that occupies all signal bands, 1D MLE is applied to determine the time delay for the reason of computational complexity. Since there are less received samples for the designed sparse multiband signal, which reduces the computational burden, the time delay is estimated through both the 1D ML and 2D ML method.
As shown in Fig. 10(a), the difference between the uncalibrated propagation distance estimated from the sparse multiband signal and the one using all signal bands is small. The curve based on the delay estimates obtained from 2D MLE is less peaky than the others, since a part of the bias is mitigated by considering an extra path in the simplified model. The performance will be evaluated through the position solutions which are compared with the ground truth values determined by the Total Station. It should be mentioned that, without a calibration, the range derived from the time delay estimate contains a hardware delay bias.
Secondly, the carrier phase can be estimated once we construct the design matrix A(τ ). Here, we also use the simplified model, as presented in Section IV-B, which only contains the LoS component a(τ 1 ). As shown in Fig. 6(c), when the time delay estimation error is relatively small, the resulting phase bias also remains small. Hence, we simply construct a(τ 1 ) with the biased delay estimate derived from the 1D ML simplified model. To even reduce the computational complexity further, the time delay may not be updated for each received signal package. Therefore, an extra delay bias will be introduced by the movement of the receiver. In this experiment, since the receiver was moving slowly (a few centimetres at most within one second), the change of the delay within one second is much smaller than the inverse of the virtual signal bandwidth (i.e., about 5 ns), therefore we keep the design matrix a(τ 1 ) constant and estimate the carrier phase for every package within this one second interval, and update a(τ 1 ) only once every second.
The carrier phase estimate is ambiguous and ranges from −π to π. Therefore, each carrier phase estimate carries its own integer phase cycle ambiguity. To avoid rank deficiencies in the positioning model, which will be introduced in the next subsection, the carrier phase should be continuously tracked and the carrier phase estimates have to be compensated for changes of the phase cycle due to movement of the receiver, referred to as phase unwrapping, so that only the initial carrier phase integer cycle ambiguity remains as an unknown parameter in the measurements. Readers are referred to [43] for more details. As long as a phase jump can be properly identified and a series of carrier phase estimate can be correctly unwrapped (i.e., no cycle-slips occur), we can use the carrier phase to compute position solutions.
The properly unwrapped carrier phase estimates in units of length are presented in Fig. 10(b). The estimates from these two signal patterns are very close, and the performance will be reflected in the position solutions.

C. Single-Differenced Relative Positioning
Here a single differenced relative positioning model is used to compute the position solutions. By computing the difference between the measurements (either the propagation delay or the carrier phase) taken at two different epochs in time from the same transmitter, unknown parameters like hardware delay and integer carrier phase cycle ambiguity, will be eliminated except for the coordinates of the receiver. The location of the starting point (x r0 , y r0 , z r0 ) is measured a-priori by the Total Station. In addition, since the receiver is moved over a horizontal rail, the height does not change. A 2D positioning scenario is therefore considered, and the height of the receiver antenna is fixed to z r0 .
First, we compute the relative position solutions based on the time delay estimates. The observation model based on the propagation time delay is given by where ρ i r0 denotes the geometric distance between the Tx-i with known coordinates (x i , y i , z i ) and the receiver at the epoch t 0 , and is defined by ρ b denotes the ranging bias due to the hardware delay and is assumed to be identical for all transmitters. In order to eliminate the hardware delay bias in (48), we compute the difference between the measurements taken at two different epochs (e.g., at t 0 and t k ) and construct the positioning model as follows, Since the height of the receive antenna z rk is known and equal to z r0 , the unknown parameters in this model are x rk and y rk .
To obtain the position solution, linearization based on Taylor expansion and Gauss-Newton iteration are applied here to solve this nonlinear model [26], [44]. The observation model for the carrier phase measurement is given by  where λ denotes the wavelength of the central carrier f c , N i denotes the integer phase cycle ambiguity of the actual propagation distance between Tx-i and the receiver, N b denotes the integer phase cycle ambiguity for the hardware delay and θ i denotes the initial phase offset for the Tx-i. Similarly, by computing the difference between the measurements taken at two different epochs, the relative positioning model based on the carrier phase estimates is written by in which the unknown parameters are again only the 2D coordinates of the receiver (i.e., x rk and y rk ). The relative position solutions using the full signal bandwidth and the sparse multiband signal based on the propagation time delay estimates, are shown in Fig. 11. The root-mean square error (RMSE) of the position solutions is shown in Table. I. The RMSE in the x and y directions are only evaluated when the receiver was at one of the five reference points on the rail, since their locations are accurately measured by the Total Station and can serve as ground truth values. In addition, a linear track for the rail is determined based on the five reference points using least-squares estimation (LSE) (shown in Fig. 11 with the green dashed-line). Afterwards, the across track error, which is the orthogonal distance between the obtained solutions and the interpolated linear track, are computed for all position solutions.
The positioning performance based on the sparse multiband signal is close to the one based on the entire available bandwidth.
In this respect, it should be noted that the antenna used in this system contains a ground-plate which filters out certain close-in multipath (e.g., from a ground bounce, or a reflection from the ceiling). For the sparse multiband signal, a resolvable reflection in 2D MLE is mostly far away from the LoS path (e.g., with a relative propagation distance larger than 1.5 m). According to the measure of bias shown in Fig. 8(b), the consequent time delay bias is small even if we do not consider this reflection in the estimation model (e.g., 1D MLE). Therefore, the improvement of using the simplified model that contains two paths (i.e., 2D MLE) becomes less obvious. Overall, a decimetre level positioning accuracy can be achieved based on the sparse multiband signal shown in Fig. 8(a), in the given environment and the used system set-up.
The position solutions, using the carrier phase estimates, shown in Fig. 12(a) are well in line with the track of the rail. The receiver stopped at each reference point for a few seconds, which can provide a better evaluation of the positioning performance. Additionally, Fig. 12(b) shows the positioning solutions versus time, in which the solutions are also well aligned with the reference points when the receiver was static.
Similarly, the RMSE of the position solutions in x-direction and y-direction are computed when the receiver was at the reference points on the rail. As shown in Table I, using the sparse multiband signal, the RMSE in x and y direction are 5.7 mm and 5.1 mm, respectively, which is close to the values when using all signal bands. And the across track error using the sparse multiband signal is about 5.3 mm.

VIII. CONCLUSION
In this paper, we have proposed a methodology to sparsely select a few signal bands within a large virtual signal bandwidth for positioning. The designed ranging signal should enable time delay and carrier phase estimation satisfying the requirements on precision and bias in a multipath environment, with reduced computational complexity. Using the Maximum Likelihood (ML) principle, in a full model that considers not only the LoS component but also the reflections, precision is sacrificed for unbiasedness due to the dependent paths. The measure of dependence, which indicates the dependence between two paths, is introduced to present how the precision deteriorates by additionally considering a reflection in the model. Alternatively, a simplified model can be used to preserve the precision and reduce the computational complexity, but at the cost of introducing a bias which is quantified by the measure of bias. Therefore, the precision of the simplified model, the measure of dependence and the measure of bias in the user-specified environment are used as the criteria to formulate an optimization problem of sparsely selecting signal bands for positioning. The experimental results show that using the sparse multiband signal, for example 7 out of 16 signal bands, with the simplified model, a decimeter level positioning accuracy is achieved based on time delay estimates, and a millimeter level accuracy based on carrier phase estimates. The resulting computational complexity is reduced by about 80%, compared with the case of using all available signal bands. For the purpose of positioning, time delay and carrier phase estimation based on the simplified model through the sparse multiband signal largely reduces the computational complexity, and yields a high precision and with only a small bias.

APPENDIX A FLOP COUNT
The computational complexity is evaluated by the number of required flops. First, we compute the required flops for time delay estimation, which is based on (8). Given M a signal bands with N subcarriers in each of the bands, and considering L paths in the estimation model, the design matrix A(τ ) becomes an M a N -by-L matrix.
It should be noted that instead of computing the number of the flops separately for the imaginary part and the real part, the complex value component is treated as a single term when we compute the number of flops. In addition, for notation simplicity, the variable τ is removed in the following derivations.
One can first analyse the computational complexity of the following term which is part of the cost function (8), and can be rewritten by Given an M a N -by-M a N variance matrix Q H , the inverse of such a matrix requires 2(M a N ) 3 /3 flops [45]. However, it is only computed once, and will be applied for both time delay and carrier phase estimation. Therefore, the number of required flops for Q −1 H is not taken into consideration in the following derivations.
Let z = C H B, one has Cz = r. As shown in [34], to obtain an L-by-1 vector z requires L 2 flops. Afterwards, another L 2 flops are needed to obtain B from z.
In Therefore, to compute the cost function (8) for each grid point, the number of required flops is derived by flops tde = (1 + 4 L)(M a N ) 2 + 2(L 2 + L)(M a N ) If there is only a single path considered in the simplified model for time delay estimation, the number of flops is given by Then, the number of flops required for complex gain estimation (5), from which one can derive the LoS carrier phase, is identical to the one in B as shown in (52), and is given by Similarly, if there is only path considered for carrier phase estimation, the number of flops is derived as flops cge = 4(M a N ) 2 + 2(M a N ), L = 1.