Joint Ranging and Clock Synchronization for Dense Heterogeneous IoT Networks

TL;DR

This paper introduces a unified framework for joint range and clock-skew estimation in dense IoT networks, addressing the limitations of existing methods by incorporating clock-skew effects into the measurement model and proposing an efficient estimation algorithm.

Contribution

It develops a novel data model that includes clock-skew effects and formulates joint estimation as a 2-D frequency problem, along with a new protocol and WLS algorithm for accurate synchronization and ranging.

Findings

The proposed estimator is asymptotically efficient, approaching the Cramer Rao lower bound.

The joint estimation method outperforms existing estimators that ignore clock-skew effects.

Numerical experiments validate the effectiveness of the proposed approach.

Abstract

Synchronization and ranging in internet of things (IoT) networks are challenging due to the narrowband nature of signals used for communication between IoT nodes. Recently, several estimators for range estimation using phase difference of arrival (PDoA) measurements of narrowband signals have been proposed. However, these estimators are based on data models which do not consider the impact of clock-skew on the range estimation. In this paper, clock-skew and range estimation are studied under a unified framework. We derive a novel and precise data model for PDoA measurements which incorporates the unknown clock-skew effects. We then formulate joint estimation of the clock-skew and range as a two-dimensional (2-D) frequency estimation problem of a single complex sinusoid. Furthermore, we propose: (i) a two-way communication protocol for collecting PDoA measurements and (ii) a weighted least squares (WLS) algorithm for joint estimation of clock-skew and range leveraging the shift invariance property of the measurement data. Finally, through numerical experiments, the performance of the proposed protocol and estimator is compared against the Cramer Rao lower bound demonstrating that the proposed estimator is asymptotically efficient.

Equations

$${\nu }_0={\nu }_1(1+{\eta }_o),$$

$${\mathcal{F}}_i={\{f_i^{(k)}=f_i^{(1)}+(k-1)\Delta f_i:\{f_i^{(1)},\Delta f_i\}\in \mathbb{R}\}}_{k=1}^K,$$

$$s_0^{(k)}(t)=\Re \left\{{\underline{s}}_0{e}^{j\left(2\pi f_0^{(k)}t+{\phi }_0^{(k)}\right)}\right\},$$

$$r_{01}^{(k)}(t)={\beta }_{01}^{(k)}\underline{s}{e}^{j\left(2\pi {\mu }_{01}^{(k)}t+{\delta }_{01}^{(k)}\right)}+n_1^{(k)}(t),$$

$${\psi }_0^{(k)}$$

$${\psi }_1^{(k)}=-2\pi f_0^{(k)}{\tau }_{01}+{\delta }_{01}^{(k)},$$

$$\text{arg}\{a_k\}=2\pi {\mu }_{01}^{(k)}\Delta t+2\pi (f_1^{(k)}+f_0^{(k)}){\tau }_{01}.$$

$${\mu }_{01}^{(k)}=G^{(1)}{\nu }_1{\eta }_0+(k-1){\eta }_0\Delta f_1,$$

$$f_1^{(k)}+$$

$$a({\tau }_{\eta })=a({\tau }_{\eta }){\left[1,e^{j2\pi \Delta f_1{\tau }_{\eta }},\dots ,e^{j2\pi (N-1)\Delta f_1{\tau }_{\eta }}\right]}^T,$$

$$a_k({\eta }_o,{\tau }_{01})=a({\tau }_{\eta }){\gamma }^{k-1}[1,\varphi ,\dots ,{\varphi }^{P-1}{]}^T,$$

$$\mathbb{M}=[m_1,\dots ,m_N].$$

$$\mathbb{M}=\mathbb{A}+\mathbb{N},$$

$$\mathbb{A}=q({\eta }_0,{\tau }_{01})h^T({\eta }_0,{\tau }_{01}),$$

$$q=$$

$$h$$

$$u_1=1/{\rho }_{q}q,v_1=1/{\rho }_h{\mathbf{h}}^{\ast },$$

$$\text{centering}{\mathbb{J}}_{\varphi 1}=[{\mathbb{I}}_{P-1}{0}_{P-1}],$$

$${\mathbb{J}}_{\varphi 2}=[0_{P-1}{\mathbb{I}}_{P-1}],$$

$$u_{\varphi 2}\approx u_{\varphi 1}\varphi \text{and}{v}_{\gamma 2}\approx v_{\gamma 1}{\gamma }^{\ast }.$$

$$\hat{\varphi }=\underset{\varphi }{\mathrm{argmin}}\Vert {\mathbb{C}}_{\varphi }^{-1/2}(u_{\varphi 1}\varphi -u_{\varphi 2}){\Vert }_2^2$$

$$\hat{\gamma }=\underset{\gamma }{\mathrm{argmin}}\Vert {\mathbb{C}}_{\gamma }^{-1/2}(v_{\gamma 1}{\gamma }^{\ast }-v_{\gamma 2}){\Vert }_2^2,$$

$${\mathbb{W}}_{\varphi }[p,n]=(P\text{min}(p,n)-pn){\varphi }^{(p-n)}/P$$

$${\mathbb{W}}_{\gamma }[p,n]=(N\text{min}(p,n)-pn){\gamma }^{(n-p)}/N,$$

$${\hat{\eta }}_o=(2\pi \Delta f_1\Delta t{)}^{-1}\text{arg}(\hat{\varphi })$$

$${\hat{d}}_{01}$$

$$\text{var}(\hat{\mathbf{\theta }})\ge {\mathbb{F}}^{-1}$$

$${\mathbb{F}}_{p,k}=2{\mathbf{\sigma }}^{-2}\Re \left[\frac{\partial a^H}{\partial {\theta }_p}\frac{\partial a}{\partial {\theta }_k}\right]\in {\mathbb{R}}^{2\times 2}$$

$$\text{var}({\hat{\eta }}_0)\ge \frac{6}{SNR(2\pi \Delta f_1\Delta t{)}^{2}PN(P^2-1)},$$

$$\text{var}({\hat{d}}_{01})\ge \frac{6c^2}{SNR(4\pi \Delta f_1{)}^{2}PN(N^2-1)}.$$

Full text

I Introduction

Synchronization and localization are key requirements of future internet of things (IoT) applications. IoT networks enable distributed information processing tasks such as sensing, aggregation, and other tasks which benefit from node location information and network-wide synchronization [1, 2]. Typically, IoT nodes are low-power devices equipped with low-cost reference clock sources, i.e. local oscillators, and narrowband radio chips. The individual clocks of the nodes drift from each other due to local oscillator imperfections, environmental and voltage variations. Therefore, it is essential to periodically synchronize and calibrate the clocks in order to keep the time synchronization among the nodes in the network. Clock drifts have a direct impact on various IoT applications, for instance, on the range estimation between the nodes, which is a crucial input for most network localization techniques.

Clock synchronization and node localization in IoT networks have received considerable attention in the past. Many research efforts have approached these problems as either separate or joint estimation problems [3, 4, 5, 6, 7, 8, 9]. Existing methods can be classified into (i) time-stamping methods based on ultra-wideband (UWB) signals [3, 4, 5, 6], and (ii) phase-based methods which utilize carrier phase measurements of narrowband signals [7, 8, 9]. Methods falling in the first class offer high timing resolution and a plethora of protocols and algorithms for joint ranging and clock synchronization have been proposed [4, 5, 6]. However, in general, these methods are not applicable to IoT networks due to the narrowband radio constraints of the nodes. On the other side, the methods based on phase-based ranging, i.e. phase difference of arrival (PDoA), do not consider the impact of unknown clock-skew on range estimation [7, 8, 9]. These methods are based on a simplified and inaccurate data model which results in a biased range estimation due to the influence of clock-skew.

In this paper, we aim for joint clock-skew and range estimation from PDoA measurements. We derive a novel and precise data model which considers hardware imperfections of the IoT nodes, i.e. clock-skew of the local oscillators, and wireless channel effects. Therefore, we propose a two-way communication protocol for collecting PDoA measurements over a two-dimensional (2-D) set of equispaced time epochs and carrier frequencies. With this data, a matrix whose rows collect measurements acquired on the same carrier frequencies but different time epochs is constructed. This data matrix exhibits a structure that allows 2-D frequency estimation techniques. Furthermore, we show that the data matrix is rank one and that its principal singular vectors have a shift invariance property which enables joint estimation of the clock-skew and range.

We propose an algorithm for joint clock-skew and range estimation based on weighted least squares (WLS) using the ideas of 2-D frequency estimation [10, 11, 12]. In this approach, the shift invariance of the left singular vector provides the range estimate, while the shift invariance of the right singular vector provides the clock-skew estimate. Finally, the performance of the proposed protocol and estimator is compared against the Cramér Rao lower bound (CRLB) using numerical simulations demonstrating that the proposed estimator is asymptotically efficient and approaches the CRLB for sufficiently high signal-to-noise ratio (SNR).

Notation: Upper (lower) bold face letters are used to denote matrices (column vectors), while $(.)^{T}$, $(.)^{H}$, $(.)^{*}$, $\odot$, $\mathbb{I}_{N}$ and $\mathbb{0}_{N}$ respectively represent transpose, Hermitian transpose, complex conjugate, element-wise Hadamard product, $N\times N$ identity matrix and $N\times 1$ vector of zeros. Furthermore, $\widehat{(.)}$ denotes estimate of a parameter, $\mathbb{E}(.)$ is the expectation operator and vec(.) forms a vector from a matrix by stacking the columns of the matrix.

II Problem Formulation and System Model

Without loss of generality, consider a single sensor (node 0) and anchor node (node 1) in a fully asynchronous wireless IoT network, as shown in Fig. 1. Let us, assume that the anchor node has a relatively stable clock oscillator and known position, while the sensor node has an unknown position and a non-ideal oscillator with frequency drift. The clock behavior of the sensor node is considered to be characterized by the first-order affine clock model [13]

[TABLE]

where $\eta_{o}$ is the clock-skew of the sensor node measured in parts per million (ppm), while $\nu_{0}$ and $\nu_{1}$ are the frequencies of the oscillator signals at the sensor and anchor node, respectively.

Here, we assume that the nodes are equipped with narrowband radio transceivers allowing two-way communication. In addition, the radio transceivers support estimation of the phase difference between the carrier frequency of the received signal and its own local oscillator frequency.

For simplicity, consider that the nodes are distributed over a two-dimensional space. Let the vectors $\mathbb{x}_{i}\in\mathbb{R}^{2\times 1}$, $i=0,1$ collect the coordinates of the nodes, where the coordinates of the sensor node $\mathbb{x}_{0}$ are unknown. The range between the anchor and sensor node is defined as ${d_{01}=d_{10}=\|\mathbb{x}_{1}-\mathbb{x}_{0}\|_{2}}$, where $\|\cdot\|_{2}$ denotes the Euclidean norm.

Frequency synthesizer model. In order to transmit a signal in a desired frequency band, each radio transceiver generates a carrier signal. The carrier signal is generated by the frequency synthesizer of the transceiver which is driven by the clock signal of the local oscillator. Modern radio transceivers support communication on a number of carrier frequencies which can be selected by changing the gain of the divider in the frequency synthesizer.

Therefore, we can assume that all the nodes have the same frequency synthesizer with a set of $K$ equispaced gains defined as $G^{(k)}=G^{(1)}+(k-1)\Delta G$, $k=1,\ldots,K$, where $G^{(k)}\in\mathbb{Q}$ is the $k$th gain and $\Delta G$ is the step of the frequency divider. The carrier frequency generated at the output of the frequency synthesizer for the $k$th gain is given by $f_{i}^{(k)}=G^{(k)}\nu_{i}$, $i=0,1$ [14]. The set of all equispaced carrier frequencies supported by the frequency synthesizer can be written as

[TABLE]

where $\Delta f_{i}=\Delta G\nu_{i}$, $i=0,1$ is the step of the frequency synthesizer, and it depends on the clock oscillator signal frequency (1).

Signal model. Consider that the sensor node transmits a single tone unmodulated carrier signal at the $k$th carrier frequency

[TABLE]

where $\underline{s}_{0}\in\mathbb{R}$ is the amplitude of the complex envelope of $s_{0}^{(k)}(t)$ and $\varphi_{0}^{(k)}$ is the unknown phase offset introduced by the process of switching the carrier frequency [15].

The transmitted signal $s_{0}^{(k)}(t)$ is narrowband, and therefore, it is reasonable to consider flat-fading effects in the channel model. The signal received at the anchor node after propagation through the channel and down-conversion by $f_{1}^{(k)}$ is given by

[TABLE]

where $\beta_{01}^{(k)}\in\mathbb{C}$ is the complex path attenuation of the channel at $f_{0}^{(k)}$, $\mu_{01}^{(k)}=f_{0}^{(k)}-f_{1}^{(k)}$ and $\delta_{01}^{(k)}=\varphi_{0}^{(k)}-\varphi_{1}^{(k)}$ are the unknown $k$th carrier frequency and phase offsets, respectively, while $n_{1}^{(k)}(t)\backsim\mathcal{CN}(0,\sigma_{1}^{2})$ denotes the zero-mean complex Gaussian noise present at the anchor node. The complex path attenuation is defined as $\beta_{01}^{(k)}=\alpha_{01}^{(k)}e^{-j2\pi f_{0}^{(k)}\tau_{01}}$ where $\alpha_{01}^{(k)}\in\mathbb{R}_{+}$ is the channel attenuation, ${\tau}_{01}=d_{01}/c=d_{10}/c$ is the unknown propagation delay between two nodes and $c$ is the known propagation speed of the radio signal. Using the frequency synthesizer and clock models, the carrier frequency offset is given by ${\mu_{01}^{(k)}=\mu_{01}^{(1)}+(k-1)\Delta\mu_{01}}$ where ${\Delta\mu_{01}=\Delta f_{1}\eta_{o}}$.

The objective in this paper is then to estimate the unknown parameters $\eta_{o}$ and $d_{01}$ given the two-way communication between nodes and PDoA functionalities of the radio transceivers.

III Communication Protocol and Data Model

In the following, we first derive a detailed data model for PDoA measurements considering a classical two-way protocol for ranging. Then, based on the derived model we propose a novel 2-D PDoA protocol for joint ranging and synchronization.

III-A *Classical PDoA ranging protocol *

In the classical two-way PDoA protocol (cf. Fig. 2a) the sensor node initiates the communication and sends a message using the signal $s_{0}^{(k)}(t)$, i.e. using carrier frequency $f_{0}^{(k)}$, to the anchor node. Then, the anchor node receives the message as the signal $r_{01}^{(k)}(t)$ and replies back to the sensor node by sending a message using signal $s_{1}^{(k)}(t)$, i.e. using carrier frequency $f_{1}^{(k)}$. After the exchange, both nodes change their carrier frequencies to $f_{i}^{(k+1)}=f_{i}^{(k)}+\Delta f_{i}$, $i=0,1$, and the same two-way message exchange pattern is repeated. The phase difference of the carrier signals, $\psi_{0}^{(k)}$ and $\psi_{1}^{(k)}$, using the $k$th carrier frequency are measured at both sensor and anchor nodes, respectively.

Now, considering the noiseless case and assuming that channel reciprocity111The received signals at anchor and sensor nodes differs only in the signs of phase and frequency offset. conditions hold, using (4) $\psi_{0}^{(k)}$ and $\psi_{1}^{(k)}$ are given by

[TABLE]

where $\Delta t$ is the deterministic time epoch between measurements collected at anchor and sensor nodes, while all other nondeterministic timing differences between nodes are absorbed in $\delta_{01}^{(k)}$. In general, the time epoch $\Delta t$ is controllable by the anchor node and it has values in the order of tens of microseconds. In the classical PDoA protocol, it is assumed that $\Delta t$ is fixed during the recollection of the measurements.

In this paper, we focus on indoor localization scenarios where the channel coherence time is typically of the order of several hundreds of milliseconds [16]. Hence, we can assume that $N\leq K$ two-way messages have been exchanged according to the PDoA protocol within the channel coherence time. For the sake of simplicity, the $N$ phase difference measurements recorded at sensor and anchor nodes are transformed in their negative complex exponential form and collected in the vectors ${\mathbb{b}_{0}=\left[e^{-j\psi_{0}^{(1)}},\ldots,e^{-j\psi_{0}^{(N)}}\right]^{T}\in\mathbb{C}^{N\times 1}}$, ${\mathbb{b}_{1}=\left[e^{-j\psi_{1}^{(1)}},\ldots,e^{-j\psi_{1}^{(N)}}\right]^{T}\in\mathbb{C}^{N\times 1}}$.

For ranging purposes, the phase offset represent nuisance parameter which can be eliminated from the acquired measurements by considering ${\mathbb{a}=\mathbb{b}_{0}\odot\mathbb{b}_{1}}$ instead. The argument of the $k$th element in $\mathbb{a}$ is given by

[TABLE]

Using the frequency synthesizer model and (5), we can write

[TABLE]

Therefore, the vector $\mathbb{a}$ has the model

[TABLE]

where ${a(\tau_{\eta})=e^{j2\pi G^{(1)}\nu_{1}\tau_{\eta}}}$ is the the first element in $\mathbb{a}(\tau_{\eta})$ and ${\tau_{\eta}=\eta_{o}\Delta t+(2+\eta_{o})\tau_{01}}$. Note that $\mathbb{a}(\tau_{\eta})$ has a shift invariance structure. This structure is precisely the one that has a uniform linear array (ULA) response vector in array processing [17]. However, in this case the phase shift of the elements in $\mathbb{a}(\tau_{\eta})$ is caused by equispaced carrier frequency switching, i.e. frequency hopping.

III-B 2-D PDoA ranging and synchronization protocol

The shift invariance of $\mathbb{a}(\tau_{\eta})$ only allows for the estimation of a single parameter $\tau_{\eta}$. However, $\eta_{o}$ and $\tau_{01}$, i.e. $d_{01}$, cannot be uniquely determined from $\mathbb{a}(\tau_{\eta})$. For example, estimation of the $\tau_{01}$ from $\mathbb{a}(\tau_{\eta})$ results in an estimate biased by the clock-skew. To alleviate this, here, we are interested in a protocol for collecting measurements that allows joint clock-skew and range estimation.

In the classical PDoA protocol, measurements are collected over the set of equispaced carrier frequencies while the time epoch $\Delta t$ is fixed during message exchange. In the 2-D PDoA protocol, we propose to collect the measurements over a 2-D set of equispaced time epochs and carrier frequencies (cf. Figs. 2b and 2c). In this case, the sensor node transmits a single message per two-way exchange, while the anchor node transmits $P$ messages based on the equispaced time epochs. The set of the equispaced time epochs for the $k$th carrier frequency is given by $\Delta t^{(k,p)}=p\Delta t/k$, where $k=1,\ldots,N$ and $p=1,\ldots,P$. Note that these time epochs depend on the index of the carrier frequency, i.e. $k$.

The $P$ phase difference measurements recorded at the sensor node for the $k$th carrier frequency are transformed in their negative complex exponential form and collected in the vector $\mathbb{b}_{k}\in\mathbb{C}^{P\times 1}$. As before, we follow a similar approach for nuisance parameters elimination. The vector that collects noiseless PDoA measurements recorded at the $k$th carrier frequency is written as ${\mathbb{a}_{k}=e^{-j\psi_{1}^{(k)}}\mathbb{b}_{k}\in\mathbb{C}^{P\times 1}}$ which satisfies the model

[TABLE]

where $a(\tau_{\eta})$ is defined in (8), $\gamma=e^{j2\pi\Delta f_{1}(2+\eta_{o})\tau_{01}}$ and $\phi=e^{j2\pi\Delta f_{1}\eta_{o}\Delta t}$.

Remark: (Practical implementation): The 2-D PDoA protocol requires that during a single two-way message exchange no carrier frequency switching occurs. This constraint ensures that the phase offset between two nodes remains constant during time hopping. However, there is no constraint on the frequency hopping sequence. This makes the proposed protocol attractive for implementation as an adaption of existing medium access control protocols such as time-slotted channel hopping (TSCH) or WirelessHART [18].

IV Joint Clock-skew and Range Estimation

In the following, we show how to jointly estimate clock-skew, i.e. $\eta_{o}$, and range, i.e. time delay $\tau_{01}$, from collected measurements.

The noise-corrupted version of $\mathbb{a}_{k}$ is given by $\mathbb{m}_{k}=\mathbb{a}_{k}+\mathbb{n}_{k}$, where $\mathbb{n}_{k}$ is a zero-mean complex Gaussian distributed noise vector222The phase estimation errors in the PLLs are Thikonov, i.e. von Mises distributed [19]. However, for large signal to noise ratio, the Thikonov distribution can be approximated by a Gaussian distribution.. From a set of $N$ noisy 2-D PDoA measurements, we construct a measurement matrix of size $P\times N$ as

[TABLE]

The measurement matrix satisfies the model

[TABLE]

where $\mathbb{A}=[\mathbb{a}_{1},\ldots,\mathbb{a}_{N}]$ and ${\mathbb{N}\in\mathbb{C}^{P\times N}}$ is the noise matrix. Using (9), it is straightforward to show that $\mathbb{A}$ can be modeled as

[TABLE]

where

[TABLE]

Model (11), using relation (12), resembles the signal model for 2-D frequency estimation of a single complex sinusoid in white Gaussian noise. This is a classical signal processing problem for which numerous methods have been proposed [20, 10, 11, 12]. Although the maximum likelihood estimator proposed in [20] can attain optimum performance, it has high computational requirements due to the multidimensional search. Here, we are interested in suboptimal but computationally more attractive (practical) methods. To do so, inspired by [12], we develop an algorithm for joint clock-skew and range estimation.

From (12), we can observe that $\mathbb{A}$ has rank one and that the vectors $\mathbb{q}$ and $\mathbb{h}$ span its column and row space, respectively. Since $\mathbb{q}$ and $\mathbb{h}$ exhibit shift invariance, it is possible to estimate $\gamma$ and $\phi$ from the low-rank approximation of $\mathbb{M}$. Then, from $\phi$ and $\gamma$, the parameters $\eta_{o}$ and $\tau_{01}$, i.e. $d_{01}$, immediately follow.

In particular, let $\mathbb{u}_{1}$ and $\mathbb{v}_{1}$ be the principal orthonormal basis vectors for the column and row span of the rank-one approximation of $\mathbb{M}$, respectively. These vectors can be obtained using the singular value decomposition (SVD) of $\mathbb{M}$ and can be expressed as

[TABLE]

where $\rho_{q}$ and $\rho_{h}$ are unknown complex constants. Now, let us define the selection matrices:

[TABLE]

To estimate $\phi$, we take subvectors consisting of the first and, respectively, the last $P-1$ elements of the $\mathbb{u}_{1}$. That is, we consider ${\mathbb{u}_{\phi 1}=\mathbb{J}_{\phi 1}\mathbb{u}_{1}}$ and ${\mathbb{u}_{\phi 2}=\mathbb{J}_{\phi 2}\mathbb{u}_{1}}$, respectively. We follow the same process for the estimation of $\gamma$, i.e. we take subvectors ${\mathbb{v}_{\gamma 1}=\mathbb{J}_{\gamma 1}\mathbb{v}_{1}}$ and ${\mathbb{v}_{\gamma 2}=\mathbb{J}_{\gamma 2}\mathbb{v}_{1}}$, respectively. From the shift invariance property of $\mathbb{u}_{1}$ and $\mathbb{v}_{1}$ we have that

[TABLE]

In the case of white noise, the approximate solutions to the relations in (16) can be found using least squares (LS). However, here we are adopt the weighted least squares (WLS) approach [11] and formulate problem (16) as

[TABLE]

where ${\mathbb{C}_{\phi}=\mathbb{E}(\mathbb{r}_{\phi}\mathbb{r}_{\phi}^{H})}$ and ${\mathbb{C}_{\gamma}=\mathbb{E}(\mathbb{r}_{\gamma}\mathbb{r}_{\gamma}^{H})}$ are the covariance matrices of the residuals ${\mathbb{r}_{\phi}=\mathbb{u}_{\phi 1}\phi-\mathbb{u}_{\phi 2}}$ and ${\mathbb{r}_{\gamma}=\mathbb{v}_{\gamma 1}\gamma^{*}-\mathbb{v}_{\gamma 2}}$, respectively. Therefore, the weighting matrices are the inverse of the covariance of the residuals, i.e. ${\mathbb{W}_{\phi}=\mathbb{C}_{\phi}^{-1}}$ and ${\mathbb{W}_{\gamma}=\mathbb{C}_{\gamma}^{-1}}$. The optimal $\mathbb{W}_{\phi}$ and $\mathbb{W}_{\gamma}$ for the considered problem are given in closed-form by [21]

[TABLE]

where $p=1,\ldots,P$ and $n=1,\ldots,N$. Note that $\mathbb{W}_{\phi}$ and $\mathbb{W}_{\gamma}$ depend on the unknown parameters $\phi$ and $\gamma$. Therefore, first we estimate $\phi$ and $\gamma$ using LS and then these estimates are used for construction of $\mathbb{W}_{\phi}$ and $\mathbb{W}_{\gamma}$. Finally, the WLS is used to obtain $\hat{\phi}$ and $\hat{\gamma}$. Based on the WLS estimates of $\phi$ and $\gamma$ the unknown parameters are computed as

[TABLE]

Note that first the clock-skew is estimated and later this estimate is used for the estimate of the range.

V Results

V-A Cramér Rao Lower Bound

To assess the performance of the proposed estimators (19) we derive the CRLB for joint clock-skew and range estimation using the model (11). For an unbiased estimator $\widehat{\bm{\theta}}$, the CRLB is the lower bound on the error variance, that is

[TABLE]

where ${\text{var}(\widehat{\bm{\theta}})=\mathbb{E}((\widehat{\bm{\theta}}-\bm{\theta})(\widehat{\bm{\theta}}-\bm{\theta})^{T})}$ and $\mathbb{F}$ is the Fisher information matrix. We assume that the proposed estimators (19) are approximately unbiased for sufficiently large SNR and well designed measurement matrix $\mathbb{M}$ [22].

In the case of 2-D frequency estimation of the sum of the sinusoids, the Fisher information matrix is given by [23]

[TABLE]

where $\mathbb{F}_{p,k}$ is the ($p$, $k$)th element of $\mathbb{F}$, $\sigma^{2}$ is the variance of the noise, $\partial/\partial\theta_{p}$ is the partial derivative with respect to the $p$th element of $\bm{\theta}$, $\mathbb{a}=\text{vec}(\mathbb{A})\in\mathbb{C}^{PN\times 1}$ is the vector formed by stacking the columns of $\mathbb{A}$. The resulting Fisher information matrix is invertible, so closed-form expressions for the CRLBs are given by

[TABLE]

where $SNR=\sigma^{-2}$.

V-B Simulations

In the following, simulations are used to compare the performance of the proposed protocol and algorithm with state-of-the-art estimators for the same problem. We consider two nodes, i.e. anchor and sensor, which are deployed randomly within a range of $140m$. The carrier frequency step $\Delta f_{i}$, $i=0,1$ and time epoch $\Delta t$ are set to $0.5MHz$ and $80\mu s$, respectively. The clock-skew of the sensor node $\eta_{o}$ is set to $80ppm$. The phase difference of arrival measurements are corrupted with zero-mean Gaussian noise and all results presented are averaged over $10^{3}$ independent Monte Carlo runs.

Figs. 3a and 4a show the root mean square error (RMSE) of clock-skew and range against the signal-to-noise-ratio for various estimators. The number of time and frequency hops is equal and set to $10$. All the algorithms are independently applied to the same set of PDoA measurements. As shown in the figures, the proposed algorithm outperforms the approximate iterative quadratic maximum-likelihood (AIQML) [21], weighted linear predictor (WLP) [24] and ESPRIT [10] for both clock-skew and range. Furthermore, for sufficiently high SNR the proposed algorithm is asymptotically efficient and approaches the theoretical bounds, i.e. the CRLB.

Fig. 3b and 4b show the RMSE of clock-skew and range against the SNR for a different number of PDoA measurements collected over time epochs ($P$) and carrier frequencies, ($N$). It is shown that by increasing the number of PDoA measurements collected over time epochs the accuracy of the clock-skew estimates is increased, while the accuracy of the range estimates increases with $N$.

Fig. 3c and 4c show the RMSE of clock-skew and range against the number of PDoA measurements, while SNR is set to $10dB$. In all scenarios the number of time and frequency hops is equal. Similar as in the previous scenarios, the proposed algorithm outperforms AIQML, WLP and ESPRIT. In addition, it can be seen that the proposed estimator achieves the CRLB.

VI Conclusions

In this paper, we investigated the problem of joint ranging and clock-skew estimation using PDoA measurements. A novel and precise data model for PDoA measurements is derived. The derived model enables joint clock-skew and range estimation from PDoA measurements collected over a 2-D equispaced time-frequency grid. We have proposed a novel protocol for collection of PDoA measurements and an algorithm based on WLS to jointly estimate the clock-skew and range. The presented algorithm leverages shift invariance properties of principal singular vectors of the collected measurements. The proposed estimator is asymptotically efficient and reaches the CRLB for sufficiently high SNR.