1 Introduction
Timeofarrival (TOA) estimation usually starts with the estimation of the underlying multipath communication channel. As the channel frequency response (CFR) is not bandlimited while we can only probe the channel with bandlimited signals, modeling assumptions are required. Traditionally, the channel impulse response (CIR) is modeled as an FIR filter of limited time duration, and the resulting time resolution for TOA estimation is inversely proportional to the sampling rate, i.e., to the bandwidth of the probing signal. This motivates the use of ultrawideband (UWB) systems [1, 2], but at the cost of large spectrum occupancy, high sampling, and high computational requirements at the receiver.
High resolution techniques therefore refine the channel model by considering a parametric model consisting of a small number of attenuated and delayed Diracs. Under this assumption, theoretically we only need to take an equally small number of samples in the frequency domain. The main challenge is to devise practical and robust schemes for implementing this.
In the past, many delay estimation algorithms have been proposed, and they can be classified into methods based on (i) subspace estimation
[3, 4, 5], (ii) finite rate of innovation [6, 7, 8], and (iii) compressed sampling signal reconstruction [9, 10, 11, 12, 13, 14]. Some of these methods are not quite robust to noise, while other methods require a separate receiver chain for each multipath component, which may not be practical.To improve resolution, a large frequency band (aperture) must be covered, while to limit sampling rates, the total band should not be densely sampled. This motivates the use of multiband acquisition systems, for e.g., [15] proposes estimation from a set of “dispersed” Fourier coefficients. Other methods include for e.g., bandpass sampling, multicoset sampling and modulated wideband converter (MWC) [16], where the implementation at the analog frontend is not straight forward.
In this paper, we aim at a limited complexity high resolution TOA estimation algorithm and consider coherent multiband acquisition. In a multichannel receiver, each receiver chain will coherently sample one of the subbands, which can be implemented with offtheshelf radio frequency (RF) components. By stacking the observations into Hankel matrices, the resulting data model has precisely the form of Multiple Invariance ESPRIT [17], so that the related algorithms are applicable, in particular, the Multiresolution ESPRIT algorithm [18], which was aimed at carrier frequency estimation.
Similar to [18], we propose an algorithm where the invariance structure of a single subband will provide coarse parameter estimates, while the the invariance structure of the lowest against the highest frequency subband will provide highresolution, but phase wrapped, estimates. The wrapping is resolved using the coarse estimates.
The resulting algorithm is benchmarked through simulations, by comparing its performance with the Cramér Rao Lower Bound (CRLB). The results show that the proposed approach provides high resolution estimates while reducing spectral occupancy and sampling costs compared to classical UWB approaches, paving the way for cognitive radio ranging systems.
2 Problem Formulation and Data Model
Channel model: We consider a channel model which is appropriate for modeling the multipath propagation of wideband and UWB signals. The multipath channel with propagation paths is defined by a continuoustime impulse response and its continuous time frequency transform (CTFT) as
(1) 
where we use “tilde” to represent signals at RF frequencies, and represent the gain and timedelay of the th resolvable path [19]. This model neglects the effects of frequency dependent distortions [20]. However, for the purpose of our analysis, it provides sufficient characterization of the radio signal propagation.
Continuoustime signal model: The objective in the paper is to estimate the channel parameters by probing the channel by a wideband training signal . Assume that covers (at least) separate bands , , where is the center frequency and is the bandwidth of the th subband. The CTFT of is
(2) 
The received signal is , where represents additive white Gaussian noise. The corresponding CTFT is (cf. Fig. (a)a)
(3) 
Now, consider a multibranch receiver having RF chains and complex sampling ADCs as shown in Fig. (b)b. The th RF chain bandlimits to , and performs complex downconversion to baseband, possibly followed by additional lowpass filtering before sampling at Nyquist. We model this by an equivalent lowpass filter with passband . The CTFT of the signal received in the th branch at baseband is thus given by
(4) 
where is bandlimited white Gaussian noise and , are the complex baseband equivalents of , . In particular, .
Discretetime data model: Assume for theoretical purposes, that has a finite duration , and is zero (or periodic) outside this interval.^{2}^{2}2In more general cases, a small bias will occur in the subsequent derivation. For simplicity of exposition, we will consider that the bandwidths of the signals and the sampling periods in all receiver branches are the same, that is and for all . We sample with period , and take samples on the nonzero interval such that . Let , then the point DFT of is given by^{3}^{3}3A factor is absorbed in .
(5) 
with, in particular,
Inserting the channel model (1) gives
(6) 
where .
Let us stack the samples of
into a vector
, and likewise for , , , and . The data model (5) then becomes(7) 
where denotes a pointwise multiplication. The channel model (6) can be written as
(8) 
where is the Vandermonde matrix
(9) 
and , where . Likewise,
(10) 
and , where .
Next, we apply deconvolution to the data vector (7). Assume that no entry of and is zero or close to zero.^{4}^{4}4 If there are zero entries, then we need to select a subvector of consecutive nonzero entries, and similar results will hold. As the entries of these vectors are known from training and filter design/calibration, the deconvolution step to estimate the DFT channel coefficients is written as
(11) 
which satisfies the model
(12) 
where
is a zero mean circular symmetric complex Gaussian distributed noise vector. It is common that the power spectral densities of the signal or the filters are not perfectly flat. In that case, the noise vector is not white, but the coloring is known and can be taken into account.
3 Multiresolution Delay Estimation
Our next objective is to estimate the timedelays . We begin with an algorithm for estimating these timedelays using a single frequency band, and later extend it for the multiple bands. The single band algorithm is in fact classical (cf. [3, 21, 22] and earlier references).
3.1 Single band estimation algorithm
From a single vector , we construct a Hankel matrix of size as
(13) 
Here, , and we require and . From (12), and using the shift invariance of the Vandermonde matrix (9), the Hankel matrix satisfies
(14) 
where is an submatrix of , and is a noise matrix. Furthermore,
where .
Since (14) resembles the data model of the classical ESPRIT algorithm, can be estimated by exploiting the lowrank approximation of the Hankel matrix and its shiftinvariance properties. From , the parameters immediately follow.
In particular, let be a dimensional orthonormal basis for the column span of
, obtained using the singular value decomposition, then we can write
, where is a nonsingular matrix. Next, let us define selection matrices(15) 
where
is identity matrix of size
andis a zero matrix of size
. For , and are submatrices of obtained by dropping its first and, respectively, last row. In view of the shift invariance structure of , we have(16) 
where . Finally, we form the matrix where denotes pseudoinverse.
can then be estimated directly from the eigenvalue decomposition of
, as it satisfies the model(17) 
In other words, let be an estimate of the th eigenvalue of , then the corresponding time delay estimate is . Since because , there is no phase wrapping issue here. Note that for TOA estimation, we are mostly interested in retrieving the smallest as it belongs to the lineofsight propagation, i.e., true distance.
3.2 Multiresolution estimation algorithm
The aforementioned algorithm used data from a single subband and has a limited resolution, since it is based on the shift of one row in the Hankel matrix , which results in only a small phase shift . Note that the sampling rate does not play a role in , only the total signal duration . Thus, oversampling would increase the signaltonoise ratio (SNR) but not the resolution.
The matrix is also invariant for shifts over multiple rows, and therefore, if is sufficiently large, then we can increase the resolution by considering shifts of multiple rows of . Indeed, a shift of rows using shift matrices (or by interleaving rows of [8]) leads to an estimate of . Unfortunately, phase shifts have an ambiguity of multiples of , so that approaches for increasing the resolution introduce ambiguity in the estimates for the . If is not very large, this approach is limited.
Here, we are interested in an algorithm for high resolution and unambiguous estimation of the from multiband channel estimates , where
. For simplicity of exposition, we will consider for the moment only two bands (i.e.,
), with central frequencies and . Following the procedure described in Section 3.1, we form the Hankel matrices defined in (13) and stack them in a matrixThe matrix has the model
(18) 
where , and are given in (10), and is formed by stacking on top of . Note that has a double shift invariance structure introduced by the phase shifts of the on the (i) sampling frequency within a single band, , and (ii) carrier frequency difference between two bands, , as shown in Fig. (a)a [23]. In general, the carrier frequency difference is much higher than the sampling frequency, and therefore, for . The estimation of the from will result in high resolution but ambiguous estimates, due to unknown multiples of in the phases. However, we can use the idea of multiresolution parameter estimation [18] to develop the algorithm for high resolution estimation of the without ambiguity by combining coarse and fine estimates obtained from and , respectively.
We follow a similar approach as in the previous section. Let be an orthonormal basis for the column span of , obtained using a truncated SVD. Define the selection matrices
(19)  
To estimate , we set and take submatrices consisting of the first and, respectively, the last row of each block matrix stacked in , that is and . The estimation of is based on the first and, respectively, second block matrix present in , that is and . The selected matrices have the following models:
(20)  
where and is given in (15). The Least Squares approximate solutions to the set of equations in (20) satisfy a model of the form
(21)  
Observe that and are jointly diagonalizable by the same matrix . If each submatrix in (20) has at least rows, the joint diagonalization can be computed by means of QZ iterations or Jacobi iterations [24, 25, 26]. After has been determined, the parameters and for are estimated.
Based on the phase estimates, coarse and fine timedelays of the delays are computed as
(22) 
The unknown number of cycles is determined as the best fitting integer that satisfies (22), that is,
(23) 
If the estimation errors of the and the are comparable, then the estimate based on is times more accurate and less affected by noise compared to the one based on . Therefore, the final estimate of is obtained based on , or by optimal combining of coarse and fine estimates [18].
This technique can be extended to matrices. Alternatively, we only consider pairwise estimates, gradually increasing the resolution until we are able to reliably estimate the amount of phase wraps for the largest shift ().
4 Results
4.1 Cramér Rao Lower Bound (CRLB)
We use the CRLB as a benchmark to study the performance of the algorithm derived in Section 3.2. The model (12) can be written as
(24) 
where , for are the columns of and . The is parameterized by . Under the assumption that the unknown multipath parameters and are deterministic, the CRLB for estimating of , conditioned on complex path attenuations , is given by [27]
(25) 
where , is the th column of , and . It is straightforward to extend the CRLB for the multiband case by creating the overall data model in the form (24).
4.2 Simulations
We consider a standard outdoor UWB channel model to evaluate the performance of the proposed algorithm [28]. The channel (1) has eight dominant multipath components (MPCs). The first MPC has times higher power in comparison to the second MPC. The continuous time is modeled using a GHz grid, where the channel tap delays are spaced at ps. In the simulations, we assume that the TOA is estimated using two bands with central frequencies at and bandwidths for . We use Root Mean Square Error (RMSE) as a metric for evaluation, which is obtained over independent Monte Carlo runs. These results are compared against the numerically computed CRLB (25).
In Fig. (a)a, the RMSEs of the estimated TOAs () for the first MPC are plotted against SNR for the frequency bands with bandwidths MHz. The proposed algorithm asymptotically achieves the CRLB, for increasing SNR. As expected, for larger bandwidths the proposed algorithm is more robust to noise, and offers higher resolution.
In Fig. (b)b, the RMSEs of are plotted against SNR, for various positions of the MHz wide bands. It can be seen that by increasing the distance between two bands, i.e. frequency aperture, the resolution of the increases for high SNR. However, for low SNR the proposed algorithm has better performance in scenarios where the frequency aperture is lower which is a consequence of lower error for fine timedelay estimation.
Fig. (c)c shows the RMSEs of the with respect to SNR for the following scenarios. Firstly, in scenarios S1 and S2, we consider the power of the second () and third () multipath components increased 2 or 1.5 times as compared to their value in S3, respectively. In scenarios S4 and S5 the distance between main and second () and third () MPC has been increased 2 or 3 times as compared to S3, respectively. As expected, the proposed algorithm is less robust to noise and has a lower resolution for scenarios where close MPCs have high power. It is seen, that the resolution of increases in scenarios where the main MPCs are more separated.
References
 [1] S. Gezici et al., “Localization via ultrawideband radios: a look at positioning aspects for future sensor networks,” IEEE signal processing magazine, vol. 22, no. 4, pp. 70–84, 2005.
 [2] K. Witrisal et al., “Noncoherent ultrawideband systems,” IEEE Signal Processing Magazine, vol. 26, no. 4, 2009.
 [3] A.J. Van der Veen, M. C. Vanderveen, and A. Paulraj, “Joint angle and delay estimation using shiftinvariance techniques,” IEEE Transactions on Signal Processing, vol. 46, no. 2, pp. 405–418, 1998.

[4]
X. Li and K. Pahlavan, “Superresolution TOA estimation with diversity for indoor geolocation,”
IEEE Transactions on Wireless Communications, vol. 3, no. 1, pp. 224–234, 2004.  [5] M. Pourkhaatoun and S. A. Zekavat, “Highresolution lowcomplexity cognitiveradiobased multiband range estimation: Concatenated spectrum vs. Fusionbased,” IEEE Systems Journal, vol. 8, no. 1, pp. 83–92, 2014.
 [6] I. Maravic, J. Kusuma, and M. Vetterli, “Lowsampling rate UWB channel characterization and synchronization,” Journal of Communications and Networks, vol. 5, no. 4, pp. 319–327, Dec 2003.
 [7] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE transactions on Signal Processing, vol. 50, no. 6, pp. 1417–1428, 2002.
 [8] I. Maravic and M. Vetterli, “Sampling and reconstruction of signals with finite rate of innovation in the presence of noise,” IEEE Transactions on Signal Processing, vol. 53, no. 8, pp. 2788–2805, 2005.
 [9] K. M. Cohen et al., “Channel estimation in UWB channels using compressed sensing,” in Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on. IEEE, 2014, pp. 1966–1970.
 [10] M. Mishali and Y. C. Eldar, “Blind multiband signal reconstruction: Compressed sensing for analog signals,” IEEE Transactions on signal processing, vol. 57, no. 3, pp. 993–1009, 2009.
 [11] M. Mishali, Y. C. Eldar, and A. J. Elron, “Xampling: Signal acquisition and processing in union of subspaces,” IEEE Transactions on Signal Processing, vol. 59, no. 10, pp. 4719–4734, 2011.
 [12] P. Zhang et al., “A compressed sensing based ultrawideband communication system,” in Communications, 2009. ICC’09. IEEE International Conference on. IEEE, 2009, pp. 1–5.
 [13] K. Gedalyahu and Y. C. Eldar, “Timedelay estimation from lowrate samples: A union of subspaces approach,” IEEE Transactions on Signal Processing, vol. 58, no. 6, pp. 3017–3031, 2010.
 [14] K. Gedalyahu, R. Tur, and Y. C. Eldar, “Multichannel sampling of pulse streams at the rate of innovation,” IEEE Transactions on Signal Processing, vol. 59, no. 4, pp. 1491–1504, 2011.
 [15] N. Wagner, Y. C. Eldar, and Z. Friedman, “Compressed beamforming in ultrasound imaging,” IEEE Transactions on Signal Processing, vol. 60, no. 9, pp. 4643–4657, 2012.
 [16] Y. C. Eldar, Sampling theory: Beyond bandlimited systems. Cambridge University Press, 2015.
 [17] A. L. Swindlehurst et al., “Multiple invariance ESPRIT,” IEEE Transactions on Signal Processing, vol. 40, no. 4, pp. 867–881, 1992.
 [18] A. N. Lemma, A.J. Van der Veen, and E. F. Deprettere, “Multiresolution ESPRIT algorithm,” IEEE Transactions on signal processing, vol. 47, no. 6, pp. 1722–1726, 1999.
 [19] R.M. Cramer, R. A. Scholtz, and M. Z. Win, “Evaluation of an ultrawideband propagation channel,” IEEE Transactions on Antennas and Propagation, vol. 50, no. 5, pp. 561–570, 2002.
 [20] A. F. Molisch, “Ultrawideband propagation channels,” Proceedings of the IEEE, vol. 97, no. 2, pp. 353–371, 2009.
 [21] C. Qian et al., “Enhanced PUMA for directionofarrival estimation and its performance analysis,” IEEE Transactions on Signal Processing, vol. 64, no. 16, pp. 4127–4137, 2016.
 [22] J. Steinwandt, F. Roemer, and M. Haardt, “Generalized least squares for ESPRITtype direction of arrival estimation,” IEEE Signal Processing Letters, vol. 24, no. 11, pp. 1681–1685, 2017.
 [23] T. Kazaz et al., “Joint ranging and clock synchronization for dense heterogeneous IoT networks,” in 2018 52nd Asilomar Conference on Signals, Systems, and Computers, Oct 2018.
 [24] A.J. van der Veen, P. B. Ober, and E. F. Deprettere, “Azimuth and elevation computation in high resolution DOA estimation,” IEEE Transactions on Signal Processing, vol. 40, no. 7, pp. 1828–1832, 1992.
 [25] A.J. Van Der Veen and A. Paulraj, “An analytical constant modulus algorithm,” IEEE Transactions on Signal Processing, vol. 44, no. 5, pp. 1136–1155, 1996.
 [26] G. Chabriel et al., “Joint matrices decompositions and blind source separation: A survey of methods, identification, and applications,” IEEE Signal Processing Magazine, vol. 31, no. 3, pp. 34–43, 2014.
 [27] P. Stoica and A. Nehorai, “Music, maximum likelihood, and cramerrao bound,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 5, pp. 720–741, 1989.
 [28] A. F. Molisch et al., “A comprehensive standardized model for ultrawideband propagation channels,” IEEE Transactions on Antennas and Propagation, vol. 54, no. 11, pp. 3151–3166, 2006.
Comments
There are no comments yet.