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 , , , , and respectively represent transpose, Hermitian transpose, complex conjugate, element-wise Hadamard product, identity matrix and vector of zeros. Furthermore, denotes estimate of a parameter, 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 
where is the clock-skew of the sensor node measured in parts per million (ppm), while and 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 , collect the coordinates of the nodes, where the coordinates of the sensor node are unknown. The range between the anchor and sensor node is defined as , where 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 equispaced gains defined as , , where is the th gain and is the step of the frequency divider. The carrier frequency generated at the output of the frequency synthesizer for the th gain is given by , . The set of all equispaced carrier frequencies supported by the frequency synthesizer can be written as
where , 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 th carrier frequency
where is the amplitude of the complex envelope of and is the unknown phase offset introduced by the process of switching the carrier frequency .
The transmitted signal 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 is given by
where is the complex path attenuation of the channel at , and are the unknown th carrier frequency and phase offsets, respectively, while denotes the zero-mean complex Gaussian noise present at the anchor node. The complex path attenuation is defined as where is the channel attenuation, is the unknown propagation delay between two nodes and is the known propagation speed of the radio signal. Using the frequency synthesizer and clock models, the carrier frequency offset is given by where .
The objective in this paper is then to estimate the unknown parameters and 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. (a)a) the sensor node initiates the communication and sends a message using the signal , i.e. using carrier frequency , to the anchor node. Then, the anchor node receives the message as the signal and replies back to the sensor node by sending a message using signal , i.e. using carrier frequency . After the exchange, both nodes change their carrier frequencies to , , and the same two-way message exchange pattern is repeated. The phase difference of the carrier signals, and , using the 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) and are given by
where is the deterministic time epoch between measurements collected at anchor and sensor nodes, while all other nondeterministic timing differences between nodes are absorbed in . In general, the time epoch 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 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 . Hence, we can assume that two-way messages have been exchanged according to the PDoA protocol within the channel coherence time. For the sake of simplicity, the phase difference measurements recorded at sensor and anchor nodes are transformed in their negative complex exponential form and collected in the vectors , .
For ranging purposes, the phase offset represent nuisance parameter which can be eliminated from the acquired measurements by considering instead. The argument of the th element in is given by
Using the frequency synthesizer model and (5), we can write
Therefore, the vector has the model
where is the the first element in and . Note that has a shift invariance structure. This structure is precisely the one that has a uniform linear array (ULA) response vector in array processing . However, in this case the phase shift of the elements in is caused by equispaced carrier frequency switching, i.e. frequency hopping.
Iii-B 2-D PDoA ranging and synchronization protocol
The shift invariance of only allows for the estimation of a single parameter . However, and , i.e. , cannot be uniquely determined from . For example, estimation of the from 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 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. (b)b and (c)c). In this case, the sensor node transmits a single message per two-way exchange, while the anchor node transmits messages based on the equispaced time epochs. The set of the equispaced time epochs for the th carrier frequency is given by , where and . Note that these time epochs depend on the index of the carrier frequency, i.e. .
The phase difference measurements recorded at the sensor node for the th carrier frequency are transformed in their negative complex exponential form and collected in the vector . As before, we follow a similar approach for nuisance parameters elimination. The vector that collects noiseless PDoA measurements recorded at the th carrier frequency is written as which satisfies the model
where is defined in (8), and .
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 .
Iv Joint Clock-skew and Range Estimation
In the following, we show how to jointly estimate clock-skew, i.e. , and range, i.e. time delay , from collected measurements.
The noise-corrupted version of is given by , where
is a zero-mean complex Gaussian distributed noise vector222The phase estimation errors in the PLLs are Thikonov, i.e. von Mises distributed . However, for large signal to noise ratio, the Thikonov distribution can be approximated by a Gaussian distribution.. From a set of noisy 2-D PDoA measurements, we construct a measurement matrix of size as
The measurement matrix satisfies the model
where and is the noise matrix. Using (9), it is straightforward to show that can be modeled as
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  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 , we develop an algorithm for joint clock-skew and range estimation.
From (12), we can observe that has rank one and that the vectors and span its column and row space, respectively. Since and exhibit shift invariance, it is possible to estimate and from the low-rank approximation of . Then, from and , the parameters and , i.e. , immediately follow.
In particular, let and be the principal orthonormal basis vectors for the column and row span of the rank-one approximation of
, respectively. These vectors can be obtained using the singular value decomposition (SVD) ofand can be expressed as
where and are unknown complex constants. Now, let us define the selection matrices:
To estimate , we take subvectors consisting of the first and, respectively, the last elements of the . That is, we consider and , respectively. We follow the same process for the estimation of , i.e. we take subvectors and , respectively. From the shift invariance property of and we have that
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  and formulate problem (16) as
where and are the covariance matrices of the residuals and , respectively. Therefore, the weighting matrices are the inverse of the covariance of the residuals, i.e. and . The optimal and for the considered problem are given in closed-form by 
where and . Note that and depend on the unknown parameters and . Therefore, first we estimate and using LS and then these estimates are used for construction of and . Finally, the WLS is used to obtain and . Based on the WLS estimates of and the unknown parameters are computed as
Note that first the clock-skew is estimated and later this estimate is used for the estimate of the range.
V-a Cramér Rao Lower Bound
). For an unbiased estimator
, the CRLB is the lower bound on the error variance, that is
In the case of 2-D frequency estimation of the sum of the sinusoids, the Fisher information matrix is given by 
where is the (, )th element of , is the variance of the noise, is the partial derivative with respect to the th element of , is the vector formed by stacking the columns of . The resulting Fisher information matrix is invertible, so closed-form expressions for the CRLBs are given by
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 . The carrier frequency step , and time epoch are set to and , respectively. The clock-skew of the sensor node is set to . The phase difference of arrival measurements are corrupted with zero-mean Gaussian noise and all results presented are averaged over independent Monte Carlo runs.
Figs. (a)a and (a)a 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 . 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) , weighted linear predictor (WLP)  and ESPRIT  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. (b)b and (b)b show the RMSE of clock-skew and range against the SNR for a different number of PDoA measurements collected over time epochs () and carrier frequencies, (). 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 .
Fig. (c)c and (c)c show the RMSE of clock-skew and range against the number of PDoA measurements, while SNR is set to . 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.
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.
-  Y.-C. Wu, Q. Chaudhari, and E. Serpedin, “Clock synchronization of wireless sensor networks,” IEEE Signal Processing Magazine, vol. 28, no. 1, pp. 124–138, 2011.
-  M. Z. Win et al., “Efficient Multisensor Localization for the Internet of Things: Exploring a New Class of Scalable Localization Algorithms,” IEEE Signal Processing Magazine, vol. 35, no. 5, pp. 153–167, 2018.
-  Y. Wang, X. Ma, and G. Leus, “Robust time-based localization for asynchronous networks,” IEEE Transactions on Signal Processing, vol. 59, no. 9, pp. 4397–4410, 2011.
-  R. T. Rajan and A.-J. van der Veen, “Joint ranging and clock synchronization for a wireless network,” in Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 2011 4th IEEE International Workshop on. IEEE, 2011, pp. 297–300.
-  S. P. Chepuri et al., “Joint clock synchronization and ranging: Asymmetrical time-stamping and passive listening,” IEEE Signal Processing Letters, vol. 20, no. 1, pp. 51–54, 2013.
-  R. T. Rajan and A.-J. van der Veen, “Joint ranging and synchronization for an anchorless network of mobile nodes,” IEEE Transactions on Signal Processing, vol. 63, no. 8, pp. 1925–1940, 2015.
-  M. Pelka, C. Bollmeyer, and H. Hellbrück, “Accurate radio distance estimation by phase measurements with multiple frequencies,” in Indoor Positioning and Indoor Navigation (IPIN), 2014 International Conference on. IEEE, 2014, pp. 142–151.
-  G. von Zengen et al., “No-cost distance estimation using standard WSN radios,” in Computer Communications, IEEE INFOCOM 2016-The 35th Annual IEEE International Conference on. IEEE, 2016, pp. 1–9.
-  O. Oshiga, S. Severi, and G. T. de Abreu, “Superresolution multipoint ranging with optimized sampling via orthogonally designed golomb rulers,” IEEE Transactions on Wireless Communications, vol. 15, no. 1, pp. 267–282, 2016.
-  S. Rouquette and M. Najim, “Estimation of frequencies and damping factors by two-dimensional ESPRIT type methods,” IEEE Transactions on signal processing, vol. 49, no. 1, pp. 237–245, 2001.
-  H.-C. So and F. K. Chan, “A generalized weighted linear predictor frequency estimation approach for a complex sinusoid,” IEEE Transactions on Signal Processing, vol. 54, no. 4, pp. 1304–1315, 2006.
-  H.-C. So et al., “An efficient approach for two-dimensional parameter estimation of a single-tone,” IEEE Transactions on Signal Processing, vol. 58, no. 4, pp. 1999–2009, 2010.
-  E. Tzoreff, B.-Z. Bobrovsky, and A. J. Weiss, “Single Receiver Emitter Geolocation Based on Signal Periodicity With Oscillator Instability.” IEEE Trans. Signal Processing, vol. 62, no. 6, pp. 1377–1385, 2014.
-  B. Razavi, RF Microelectronics (2Nd Edition) (Prentice Hall Communications Engineering and Emerging Technologies Series), 2nd ed. Upper Saddle River, NJ, USA: Prentice Hall Press, 2011.
-  D. Tse and P. Viswanath, Fundamentals of wireless communication. Cambridge university press, 2005.
-  H. S. Rahul, S. Kumar, and D. Katabi, “JMB: scaling wireless capacity with user demands,” in Proceedings of the ACM SIGCOMM 2012 conference on Applications, technologies, architectures, and protocols for computer communication. ACM, 2012, pp. 235–246.
-  H. Krim and M. Viberg, “Two decades of array signal processing research: the parametric approach,” IEEE signal processing magazine, vol. 13, no. 4, pp. 67–94, 1996.
-  T. Watteyne, M. Palattella, and L. Grieco, “Using IEEE 802.15. 4e time-slotted channel hopping (TSCH) in the internet of things (IoT): Problem statement,” Tech. Rep., 2015.
-  Y. S. Shmaliy, “Von mises/tikhonov-based distributions for systems with differential phase measurement,” Signal Processing, vol. 85, no. 4, pp. 693–703, 2005.
-  M. P. Clark and L. L. Scharf, “Two-dimensional modal analysis based on maximum likelihood,” IEEE Transactions on Signal Processing, vol. 42, no. 6, pp. 1443–1452, 1994.
-  H.-C. So and K. Chan, “Approximate maximum-likelihood algorithms for two-dimensional frequency estimation of a complex sinusoid,” IEEE Transactions on Signal Processing, vol. 54, no. 8, pp. 3231–3237, 2006.
-  J. Liu, X. Liu, and X. Ma, “First-order perturbation analysis of singular vectors in singular value decomposition,” IEEE Transactions on Signal Processing, vol. 56, no. 7, pp. 3044–3049, 2008.
-  Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” IEEE Transactions on Signal Processing, vol. 40, no. 9, pp. 2267–2280, 1992.
-  S. Kay and R. Nekovei, “An efficient two-dimensional frequency estimator,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 38, no. 10, pp. 1807–1809, Oct 1990.