Multiresolution time-of-arrival estimation from multiband radio channel measurements

TL;DR

This paper introduces a multiband sampling method combined with a multiresolution ESPRIT-based algorithm for high-resolution time-of-arrival estimation in multipath radio channels, reducing spectral and sampling costs.

Contribution

It proposes a practical multibranch receiver and a multiresolution TOA estimation algorithm leveraging multiple shift invariance structures, advancing high-resolution multipath parameter estimation.

Findings

Achieves high-resolution TOA estimates with reduced spectral occupancy.

Performance approaches the Cramér-Rao Lower Bound in simulations.

Outperforms traditional UWB methods in resolution and efficiency.

Abstract

Achieving high resolution time-of-arrival (TOA) estimation in multipath propagation scenarios from bandlimited observations of communication signals is challenging because the multipath channel impulse response (CIR) is not bandlimited. Modeling the CIR as a sparse sequence of Diracs, TOA estimation becomes a problem of parametric spectral inference from observed bandlimited signals. To increase resolution without arriving at unrealistic sampling rates, we consider multiband sampling approach, and propose a practical multibranch receiver for the acquisition. The resulting data model exhibits multiple shift invariance structures, and we propose a corresponding multiresolution TOA estimation algorithm based on the ESPRIT algorithm. The performance of the algorithm is compared against the derived Cram\'er Rao Lower Bound, using simulations with standardized ultra-wideband (UWB) channel models. We show that the proposed approach provides high-resolution estimates while reducing spectral occupancy and sampling costs compared to traditional UWB approaches.

Equations

$$\tilde{h}(t)=\sum _{k=1}^K{\tilde{\alpha }}_k\delta (t-{\tau }_k)\text{and}\tilde{H}(\Omega )=\sum _{k=1}^K{\tilde{\alpha }}_k{e}^{-j\Omega {\tau }_k},$$

$$\tilde{S}(\Omega )=\{\begin{array}{ll}{\tilde{S}}_i(\Omega ),& \Omega \in {\mathcal{W}}_i,i=1,\cdots ,L\\ \text{mbox}arbitrary,& \text{otherwise}.\end{array}$$

$$\tilde{X}(\Omega )=\tilde{H}(\Omega )\tilde{S}(\Omega )+\tilde{N}(\Omega ).$$

$$X_i(\Omega )=\{\begin{array}{ll}{G}_i(\Omega )H_i(\Omega )S_i(\Omega )+N_i(\Omega ),& \Omega \in {\mathcal{B}}_i\\ 0,& \text{otherwise}\end{array}$$

$$X_i[n]=G_i[n]H_i[n]S_i[n]+N_i[n],n=0,\cdots ,N-1$$

$$H_i[n]=H_i\left(\frac{n}{N}{\Omega }_s\right)=\tilde{H}(\frac{n}{N}{\Omega }_s+{\Omega }_i).$$

$$H_i[n]=\sum _{k=1}^K{\tilde{\alpha }}_k{e}^{-j{\Omega }_i{\tau }_k}{e}^{-jn{\Omega }_t{\tau }_k}$$

$${\mathbf{x}}_i={\mathbf{h}}_i\odot {\mathbf{g}}_i\odot {\mathbf{s}}_i+{\mathbf{n}}_i.$$

$${\mathbf{h}}_i=\mathbf{M}{\mathbf{\Theta }}_i\mathbf{\alpha },$$

$$\mathbf{M}=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ {\Phi }_1& {\Phi }_2& \cdots & {\Phi }_K\\ ⋮& ⋮& \ddots & ⋮\\ {\Phi }_1^{N-1}& {\Phi }_2^{N-1}& \cdots & {\Phi }_K^{N-1}\end{array}\right],$$

$${\mathbf{\Theta }}_i=\left[\begin{array}{ccc}{\Theta }_{i,1}& & 0\\ & \ddots & \\ 0& & {\Theta }_{i,K}\end{array}\right],\mathbf{\alpha }=\left[\begin{array}{c}{\tilde{\alpha }}_1\\ ⋮\\ {\tilde{\alpha }}_K\end{array}\right]$$

$${\mathbf{h}}_i=\{\text{diag}({\mathbf{g}}_i\odot {\mathbf{s}}_i){\}}^{-1}{\mathbf{x}}_i,$$

$${\mathbf{h}}_i=\mathbf{M}{\mathbf{\Theta }}_i\alpha +{\mathbf{n}}_i^{\prime },$$

$${\mathcal{H}}_i=\left[\begin{array}{cccc}{H}_i[0]& H_i[1]& \cdots & H_i[Q]\\ H_i[1]& H_i[2]& \cdots & H_i[Q+1]\\ ⋮& ⋮& \ddots & ⋮\\ H_i[P-1]& H_i[P]& \cdots & H_i[N-1]\end{array}\right].$$

$${\mathcal{H}}_i={\mathbf{M}}^{\prime }{\mathbf{\Theta }}_i\mathbf{A}+{\mathbf{N}}_i,$$

$$\mathbf{A}=[\mathbf{\alpha },\mathbf{\Phi }\mathbf{\alpha },{\mathbf{\Phi }}^2\mathbf{\alpha },\cdots ,{\mathbf{\Phi }}^{Q-1}\mathbf{\alpha }]$$

$${\mathbf{J}}_1^{(1)}=[{\mathbf{I}}_{P-r}{0}_{P-r,r}],{\mathbf{J}}_2^{(1)}=[0_{P-r,r}{\mathbf{I}}_{P-r}],$$

$${\mathbf{U}}_1={\mathbf{M}}_1^{\prime }{\mathbf{T}}^{-1},{\mathbf{U}}_2={\mathbf{M}}_1^{\prime }\mathbf{\Phi }{\mathbf{T}}^{-1}$$

$$\mathbf{\Psi }=\mathbf{T}\mathbf{\Phi }{\mathbf{T}}^{-1}.$$

$$\mathcal{H}=\left[\begin{array}{c}{\mathcal{H}}_1\\ {\mathcal{H}}_2\end{array}\right].$$

$$\mathcal{H}=\left[\begin{array}{c}{\mathbf{M}}^{\prime }\\ {\mathbf{M}}^{\prime }\mathbf{\Theta }\end{array}\right]{\mathbf{\Theta }}_1\mathbf{A}+\mathbf{N},$$

$$\text{centering}{\mathbf{J}}_{\Phi 1}^{(r)}={\mathbf{I}}_2\otimes [{\mathbf{I}}_{P-r}{0}_{P-r,r}],$$

$${\mathbf{J}}_{\Phi 2}^{(r)}={\mathbf{I}}_2\otimes [0_{P-r,r}{\mathbf{I}}_{P-r}],$$

$${\mathbf{U}}_{\Phi 1}$$

$${\mathbf{U}}_{\Phi 2}$$

$$\mathbf{\Psi }$$

$$\mathbf{{\rm Y}}$$

$${\tau }_k={\Omega }_t^{-1}{\varphi }_k=({\Omega }_2-{\Omega }_1{)}^{-1}({\theta }_k+2\pi n_k).$$

$$n_k=\text{round}\left\{\frac{1}{2\pi }\left({\Omega }_t^{-1}({\Omega }_2-{\Omega }_1){\varphi }_k-{\theta }_k\right)\right\}.$$

$${\mathbf{h}}_i={\mathbf{B}}_i\mathbf{\alpha }+{\mathbf{n}}_i,$$

$$\text{CRLB}(\mathbf{\tau })=\frac{{\sigma }^2}{2}\left\{\Re {\left[{\mathbf{D}}^H{\mathbf{P}}_{{\mathbf{B}}_i}^{\perp }\mathbf{D}\odot {\mathbf{R}}_{\mathbf{\alpha }}\right]}^{-1}\right\},$$

Full text

**Index Terms— ** time-of-arrival, multiresolution estimation, cognitive radio, multiband sampling, multipath channel estimation

1 Introduction

Time-of-arrival (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 ultra-wideband (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 front-end 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 sub-bands, which can be implemented with off-the-shelf 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 sub-band will provide coarse parameter estimates, while the the invariance structure of the lowest against the highest frequency sub-band will provide high-resolution, 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 $K$ propagation paths is defined by a continuous-time impulse response $\widetilde{h}(t)$ and its continuous time frequency transform (CTFT) $\widetilde{H}(\Omega)$ as

[TABLE]

where we use “tilde” to represent signals at RF frequencies, $\widetilde{\alpha}_{k}\in\mathbb{R}$ and $\tau_{k}\in\mathbb{R}_{+}$ represent the gain and time-delay of the $k$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.

Continuous-time signal model: The objective in the paper is to estimate the $2K$ channel parameters by probing the channel by a wideband training signal $\widetilde{s}(t)$. Assume that $\widetilde{s}(t)$ covers (at least) $L$ separate bands $\mathcal{W}_{i}=[\Omega_{i}-\frac{1}{2}B_{i},\Omega_{i}+\frac{1}{2}B_{i}]$, $i=1,\cdots,L$, where $\Omega_{i}$ is the center frequency and $B_{i}$ is the bandwidth of the $i$th sub-band. The CTFT of $\widetilde{s}(t)$ is

[TABLE]

The received signal is $\widetilde{x}(t)=\widetilde{s}(t)\ast\widetilde{h}(t)+\widetilde{n}(t)$, where $\widetilde{n}(t)$ represents additive white Gaussian noise. The corresponding CTFT is (cf. Fig. 1a)

[TABLE]

Now, consider a multibranch receiver having $L$ RF chains and complex sampling ADCs as shown in Fig. 1b. The $i$th RF chain bandlimits $\widetilde{X}(\Omega)$ to $\Omega\in\mathcal{W}_{i}$, and performs complex downconversion to baseband, possibly followed by additional lowpass filtering before sampling at Nyquist. We model this by an equivalent lowpass filter $G_{i}(\Omega)$ with passband $\mathcal{B}_{i}=[-\frac{1}{2}B_{i},\frac{1}{2}B_{i}]$. The CTFT of the signal $x_{i}(t)\in\mathbb{C}$ received in the $i$th branch at baseband is thus given by

[TABLE]

where $N_{i}(\Omega)$ is bandlimited white Gaussian noise and $\{H_{i}(\Omega)$, $S_{i}(\Omega),G_{i}(\Omega)\}$ are the complex baseband equivalents of $\{\widetilde{H}(\Omega)$, $\widetilde{S}(\Omega),\widetilde{G}(\Omega)\}$. In particular, $H_{i}(\Omega)=\widetilde{H}(\Omega+\Omega_{i})$.

Discrete-time data model: Assume for theoretical purposes, that $x_{i}(t)$ has a finite duration $T$, and is zero (or periodic) outside this interval.222In 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 $B_{i}=B$ and $T_{s,i}=T_{s}$ for all $i\in[1,L]$. We sample $x_{i}(t)$ with period $T_{s}$, and take $N$ samples on the nonzero interval such that $T=NT_{s}$. Let $\Omega_{s}=2\pi/T_{s}$, then the $N$-point DFT of $x_{i}(t)$ is given by333A factor $1/T_{s}$ is absorbed in $S_{i}[n]$.

[TABLE]

with, in particular,

[TABLE]

Inserting the channel model (1) gives

[TABLE]

where $\Omega_{t}=\frac{1}{N}\Omega_{s}=\frac{2\pi}{T}$.

Let us stack the $N$ samples of $X_{i}[n]$ into a vector $\bm{x}_{i}$, and likewise for $\bm{g}_{i}$, $\bm{h}_{i}$, $\bm{s}_{i}$, and $\bm{n}_{i}$. The data model (5) then becomes

[TABLE]

where $\odot$ denotes a pointwise multiplication. The channel model (6) can be written as

[TABLE]

where $\bm{M}$ is the $N\times K$ Vandermonde matrix

[TABLE]

and $\Phi_{k}=e^{-j\phi_{k}}$, where $\phi_{k}=\Omega_{t}\tau_{k}$. Likewise,

[TABLE]

and $\Theta_{i,k}=e^{-j\theta_{i,k}}$, where $\theta_{i,k}=\Omega_{i}\tau_{k}$.

Next, we apply deconvolution to the data vector (7). Assume that no entry of $\bm{g}_{i}$ and $\bm{s}_{i}$ is zero or close to zero.444 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

[TABLE]

which satisfies the model

[TABLE]

where $\bm{n}^{\prime}_{i}$ 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 $K$ time-delays $\{{\tau_{k}}\}_{k=1}^{K}$. We begin with an algorithm for estimating these time-delays 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 $\bm{h}_{i}$, we construct a Hankel matrix of size $P\times Q$ as

[TABLE]

Here, $P=N-Q-1$, and we require $P>K$ and $Q\geq K$. From (12), and using the shift invariance of the Vandermonde matrix (9), the Hankel matrix satisfies

[TABLE]

where $\bm{M}^{\prime}$ is an $P\times K$ submatrix of $\bm{M}$, and $\bm{N}_{i}$ is a noise matrix. Furthermore,

[TABLE]

where $\mathbf{\Phi}=\text{diag}([\phi_{1}\cdots\phi_{K}])$.

Since (14) resembles the data model of the classical ESPRIT algorithm, $\bm{\Phi}$ can be estimated by exploiting the low-rank approximation of the Hankel matrix and its shift-invariance properties. From $\bm{\Phi}$, the parameters $\tau_{k}$ immediately follow.

In particular, let $\bm{U}$ be a $K$-dimensional orthonormal basis for the column span of $\bm{\mathcal{H}}_{i}$, obtained using the singular value decomposition, then we can write $\bm{M}^{\prime}=\bm{U}\bm{T}$, where $\bm{T}$ is a $K\times K$ nonsingular matrix. Next, let us define selection matrices

[TABLE]

where $\bm{I}_{P-r}$ is identity matrix of size $(P-r)\times(P-r)$ and $\mathbf{0}_{P-r,r}$ is a zero matrix of size $(P-r)\times r$. For $r=1$, $\bm{U}_{1}=\bm{J}_{1}^{(1)}\bm{U}$ and $\bm{U}_{2}=\bm{J}_{2}^{(1)}\bm{U}$ are submatrices of $\bm{U}$ obtained by dropping its first and, respectively, last row. In view of the shift invariance structure of $\bm{M}^{\prime}$, we have

[TABLE]

where $\bm{M}_{1}^{\prime}=\bm{J}_{1}^{(1)}\bm{M}^{\prime}$. Finally, we form the matrix $\bm{\Psi}=\bm{U}_{1}^{\dagger}\bm{U}_{2}$ where $\dagger$ denotes pseudo-inverse. $\bm{\Phi}$ can then be estimated directly from the eigenvalue decomposition of $\bm{\Psi}$, as it satisfies the model

[TABLE]

In other words, let $\hat{\lambda}_{k}$ be an estimate of the $k$th eigenvalue of $\bm{\Psi}$, then the corresponding time delay estimate is ${\tau}_{k}=\text{arg}\{\lambda_{k}\}/\Omega_{t}$. Since $\Omega_{t}\tau_{k}<2\pi$ because $\tau_{k}<T$, there is no phase wrapping issue here. Note that for TOA estimation, we are mostly interested in retrieving the smallest $\tau_{k}$ as it belongs to the line-of-sight propagation, i.e., true distance.

3.2 Multiresolution estimation algorithm

The aforementioned algorithm used data from a single sub-band and has a limited resolution, since it is based on the shift of one row in the Hankel matrix $\bm{\mathcal{H}}_{i}$, which results in only a small phase shift $\Omega_{t}\tau_{k}$. Note that the sampling rate does not play a role in $\Omega_{t}$, only the total signal duration $T$. Thus, oversampling would increase the signal-to-noise ratio (SNR) but not the resolution.

The matrix $\bm{M}$ is also invariant for shifts over multiple rows, and therefore, if $N$ is sufficiently large, then we can increase the resolution by considering shifts of multiple rows of $\bm{\mathcal{H}}_{i}$. Indeed, a shift of $r$ rows using shift matrices $\bm{J}_{1}^{(r)},\bm{J}_{2}^{(r)}$ (or by interleaving rows of $\bm{\mathcal{H}}_{i}$ [8]) leads to an estimate of $\bm{\Phi}^{r}$. Unfortunately, phase shifts have an ambiguity of multiples of $2\pi$, so that approaches for increasing the resolution introduce ambiguity in the estimates for the $\tau_{k}$. If $T$ is not very large, this approach is limited.

Here, we are interested in an algorithm for high resolution and unambiguous estimation of the $\tau_{k}$ from multiband channel estimates $\bm{h}_{i}$, where $i=1,\cdots,L$. For simplicity of exposition, we will consider for the moment only two bands (i.e., $i=1,2$), with central frequencies $\Omega_{1}$ and $\Omega_{2}$. Following the procedure described in Section 3.1, we form the Hankel matrices $\bm{\mathcal{H}}_{i}$ defined in (13) and stack them in a matrix

[TABLE]

The matrix $\bm{\mathcal{H}}$ has the model

[TABLE]

where $\bm{\Theta}=\bm{\Theta}_{2}\bm{\Theta}_{1}^{-1}$, $\bm{\Theta}_{1}$ and $\bm{\Theta}_{2}$ are given in (10), and $\bm{N}$ is formed by stacking $\bm{N}_{1}$ on top of $\bm{N}_{2}$. Note that $\bm{\mathcal{H}}$ has a double shift invariance structure introduced by the phase shifts of the $\tau_{k}$ on the (i) sampling frequency within a single band, $\phi_{k}$, and (ii) carrier frequency difference between two bands, $\theta_{k}=\theta_{2,k}-\theta_{1,k}$, as shown in Fig.  1a [23]. In general, the carrier frequency difference is much higher than the sampling frequency, and therefore, $\theta_{k}\gg\phi_{k}$ for $k=1,\cdots,K$. The estimation of the $\tau_{k}$ from $\bm{\Theta}$ will result in high resolution but ambiguous estimates, due to unknown multiples of $2\pi$ in the phases. However, we can use the idea of multiresolution parameter estimation [18] to develop the algorithm for high resolution estimation of the $\tau_{k}$ without ambiguity by combining coarse and fine estimates obtained from $\bm{\Phi}$ and $\bm{\Theta}$, respectively.

We follow a similar approach as in the previous section. Let $\bm{U}$ be an orthonormal basis for the column span of $\bm{\mathcal{H}}$, obtained using a truncated SVD. Define the selection matrices

[TABLE]

To estimate $\bm{\Phi}$, we set $r=1$ and take submatrices consisting of the first and, respectively, the last row of each block matrix stacked in $\bm{U}$, that is $\bm{U}_{\Phi 1}=\bm{J}_{\Phi 1}^{(1)}\bm{U}$ and $\bm{U}_{\Phi 2}=\bm{J}_{\Phi 2}^{(1)}\bm{U}$. The estimation of $\bm{\Theta}$ is based on the first and, respectively, second block matrix present in $\bm{U}$, that is $\bm{U}_{\Theta 1}=\bm{J}_{\Theta 1}\bm{U}$ and $\bm{U}_{\Theta 2}=\bm{J}_{\Theta 2}\bm{U}$. The selected matrices have the following models:

[TABLE]

where $\bm{M}^{\prime\prime}=\bm{J}_{1}^{(1)}\bm{M}^{\prime}$ and $\bm{J}_{1}^{(1)}$ is given in (15). The Least Squares approximate solutions to the set of equations in (20) satisfy a model of the form

[TABLE]

Observe that $\bm{\Psi}$ and $\bm{\Upsilon}$ are jointly diagonalizable by the same matrix $\bm{T}$. If each submatrix in (20) has at least $K$ rows, the joint diagonalization can be computed by means of QZ iterations or Jacobi iterations [24, 25, 26]. After $\bm{T}$ has been determined, the parameters $\phi_{k}$ and $\theta_{k}$ for $k=1,\cdots,K$ are estimated.

Based on the phase estimates, coarse and fine time-delays of the delays are computed as

[TABLE]

The unknown number of cycles $n_{k}$ is determined as the best fitting integer that satisfies (22), that is,

[TABLE]

If the estimation errors of the $\phi_{k}$ and the $\theta_{k}$ are comparable, then the $\tau_{k}$ estimate based on $\theta_{k}$ is $\Omega_{t}^{-1}(\Omega_{2}-\Omega_{1})$ times more accurate and less affected by noise compared to the one based on $\phi_{k}$. Therefore, the final estimate of $\phi_{k}$ is obtained based on $\theta_{k}$, or by optimal combining of coarse and fine estimates [18].

This technique can be extended to $L$ matrices. Alternatively, we only consider pairwise estimates, gradually increasing the resolution until we are able to reliably estimate the amount of $2\pi$ phase wraps for the largest shift ($\Omega_{L}-\Omega_{1}$).

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

[TABLE]

where $\mathbf{B}_{i}=\mathbf{M}\bm{\Theta}_{i}=[\mathbf{m}_{1}\Theta_{i,1},\;\mathbf{m}_{2}\Theta_{i,2},\;\dots\;\mathbf{m}_{K}\Theta_{i,K}]$, $\mathbf{m}_{k}$ for $k=1,\cdots,K$ are the columns of $\mathbf{M}$ and $\mathbf{n}_{i}\sim\mathcal{C}\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{N})$. The $\mathbf{B}_{i}$ is parameterized by $\bm{\tau}=[\tau_{1},\cdots,\tau_{K}]^{T}$. Under the assumption that the unknown multipath parameters $\bm{\alpha}$ and $\bm{\tau}$ are deterministic, the CRLB for estimating of $\bm{\tau}$, conditioned on complex path attenuations $\bm{\alpha}$, is given by [27]

[TABLE]

where $\mathbf{D}=\left[\dfrac{\partial\mathbf{b}_{i}(\tau_{1})}{\partial\tau_{1}},\dots,\dfrac{\partial\mathbf{b}_{i}(\tau_{K})}{\partial\tau_{K}}\right]$, $\mathbf{b}_{i}(\tau_{k})=\mathbf{m}_{k}\theta_{i,k}$ is the $k$th column of $\mathbf{B}_{i}$, $\mathbf{P}_{\mathbf{B}_{i}}^{\perp}=\mathbf{I}_{N}-\mathbf{B}_{i}(\mathbf{B}_{i}^{H}\mathbf{B}_{i})^{-1}\mathbf{B}_{i}$ and $\mathbf{R}_{\bm{\alpha}}=E\{\bm{\alpha}\bm{\alpha}^{H}\}$. 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 $8$ times higher power in comparison to the second MPC. The continuous time is modeled using a $3$ GHz grid, where the channel tap delays are spaced at $333.33$ ps. In the simulations, we assume that the TOA is estimated using two bands with central frequencies at $\widetilde{f}_{i}=\Omega_{i}/(2\pi)$ and bandwidths $\Delta f_{i}=B_{i}/(2\pi)$ for $i\in[1,2]$. We use Root Mean Square Error (RMSE) as a metric for evaluation, which is obtained over $10^{4}$ independent Monte Carlo runs. These results are compared against the numerically computed CRLB (25).

In Fig. 2a, the RMSEs of the estimated TOAs (${\tau}_{1}$) for the first MPC are plotted against SNR for the frequency bands with bandwidths $[160,200,300]$ 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. 2b, the RMSEs of ${\tau}_{1}$ are plotted against SNR, for various positions of the $200$ MHz wide bands. It can be seen that by increasing the distance between two bands, i.e. frequency aperture, the resolution of the ${\tau}_{1}$ 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 time-delay estimation.

Fig. 2c shows the RMSEs of the ${\tau}_{1}$ with respect to SNR for the following scenarios. Firstly, in scenarios S1 and S2, we consider the power of the second ($P_{MPC2}$) and third ($P_{MPC3}$) 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 ($\Delta_{1,2}$) and third ($\Delta_{1,3}$) 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 ${\tau}_{1}$ increases in scenarios where the main MPCs are more separated.