TDOA , Frequency and Phase Offsets Estimation Taking Into Account Carrier Phase of Arrival

This paper deals with estimation of time difference of arrival (TDOA) and phase and frequency offsets between the channels in a receiving system with two distributed time synchronized, but phase and frequency unlocked channels. The system receives a radio signal with a known waveform. We analyze the impacts of using the carrier phase of arrival (CPOA) on parameter estimation accuracy. Depending on which of the parameters are unknown, three cases are considered. For the case when phase, frequency, and time shifts are all unknown, we derive the Fisher Information Matrix (FIM) in closed form, and Cramer-Rao bound (CRB) closed form expressions for the unknown phase and frequency offsets case, and the unknown TDOA case. Two maximum likelihood type (ML) statistically efficient estimation algorithms are proposed. The CRBs and simulation results show that in the unknown TDOA case the accuracy of TDOA estimation can be significantly increased by using CPOA.


Introduction
Consider the complex representation exp(jω c t) s(t) of a radio signal (ω c is the circular carrier frequency, s(t) is the complex envelope).When delayed by any D, its baseband representation becomes exp(−jω c D) s (t − D).This paper focuses on the analysis of effects of taking into account the carrier phase term exp(−jω c D) on the performance of time difference of arrival (TDOA) and frequency and phase offsets estimation in a distributed dual-channel receiving system whose receivers are time-synchronized, but have independent local oscillators.The estimation accuracy improvement by using the carrier phases has not been systematically studied in the state-of-the-art literature, although it could be beneficial to Distributed Beamforming (DB) and Distributed Multiuser Multiple-Input Multiple-Output systems (DM MIMO) applications, [1][2][3][4][5].
The problem of time, frequency and phase shifts estimation in dual-channel systems has been extensively studied in the literature, see [6][7][8][9][10][11][12].An estimator for TDOA, phase and discrete Doppler shifts using Gaussian random signals was proposed in [7].A signal model with TDOA, arbitrary Doppler shifts, and random phase using unknown deterministic signals was analyzed in [8].In these papers and in a number of TDOA/FDOA localization related papers such as [10], the carrier phase of arrival (CPOA) term was modeled as a part of an unknown or random phase term or complex attenuation.The Root-Mean-Square (RMS) error of TDOA estimation is then inversely proportional to the complex envelope bandwidth.
Unlike this, we consider a fully coherent propagation scenario, i.e.CPOA of the signal due to propagation from the transmitter, Tx, to receiver Rx i is modeled by exp(−jω c D pi ), where D pi is the propagation time.In this scenario, if the term exp(−jω c D pi ) were modeled as a part of an unknown/random phase term/complex attenuation, a part of distance (delay) related information would be lost.Therefore, this term is a separate term in our model.Unwrapped phase ω c D pi carries information about distances between the transmitter and receivers.Since it is measured modulo 2π, this model is ambiguous.Weiss and Weinstein in [12] have analyzed TDOA estimation of analog bandpass signals, which implicitly takes into account CPOA.They have recognized the ambiguity problem and derived a composite theoretical bound on TDOA estimation accuracy.For high SNRs it coincides with the Cramer-Rao bound (CRB).However, we estimate TDOA based on digital baseband versions of those signals.
Possible solutions to the ambiguity problem were analyzed in [13] and [14].
The assumptions of our model are: the frequency shifts are due to local oscillator mismatch (not Doppler), with an arbitrary (usually smaller than the DFT resolution, unlike [7]) and constant value in the observation interval, time shifts are not integer multiples of the sampling interval in general, the distance between receiving antennas can be significantly larger than the carrier wavelength, and the signal is arbitrarily wideband and known (unlike [7], [8], [10], so phase, frequency, and time shifts are all modeled even in the referent channel).
We analyze three cases.In case 1 all three types of parameters (phase, frequency, and time shifts) are unknown, which corresponds to the over-the-air calibration of a 2channel receiving system using a transmitter at an unknown location, or to the TDOA estimation with phase and frequency non-calibrated channels.In case 2 TDOA is known (phase and frequency offsets are unknown), which corresponds to a transmitter at a known location used to calibrate the receiving system (e.g.over-the-air calibration of sensors on a master-slave principle in DB systems).In case 3 phase and frequency offsets are known, which corresponds to a TDOA estimation by a fully synchronized receiving system (a distributed receiving system calibrated by another transmitter as in case 2 or a system with collocated receivers connected via calibrated coaxial cables or Radio Frequency over Fiber (RFoF) links to distributed antennas).For example, this could be used to estimate the coefficients for a distributed beamformer.
The contributions of the paper are as follows.We derive the Fisher Information Matrix (FIM) in closed form for case 1, and we derive CRB expressions in closed form for cases 2 and 3. Two maximum likelihood-type (ML) estimation algorithms are proposed, one for the unknown sequence scenario and the other for the known sequence scenario.The algorithms are statistically efficient and partially numerically optimized.The CRBs and simulations show how using the information embedded in the carrier phase influences the accuracy of estimation for each case.The most important result is that the accuracy of TDOA estimation in case 3 can be dramatically improved, given that the correct lobe of the criterion function is disambiguated from the other lobes.In this case the RMS error is inversely proportional to the carrier frequency.
The remainder of the paper is organized as follows: In Sec. 2 our signal model is presented.Section 3 includes a derivation of CRB expressions for the three cases.Two ML-type estimation algorithms and the ambiguity problem are described in Sec. 4. The performance of the algorithms is evaluated in Sec. 5.The conclusions are given in Sec. 6.
In the rest of the paper, operators () T , () H , () * , () , Re, and E will denote transpose, Hermitian, complex conjugate, the first derivative, the real part and expectation, respectively.

Signal Model
A periodic pilot signal s(t) with period T, bandlimited to [−B/2, B/2], B is in [Hz], is upconverted and transmitted by a stationary transmitter at nominal carrier frequency ν c [Hz].It is received at two time-synchronized distributed stationary receivers, with independent local  The signal s(t) is formed in the Tx from a complex known sequence q = q 0 , q 1 , . . ., q N −1 T in the following way.Let the DFT of q be then (θ is time normalized with 1/ν s ), then s(t) = s(tν s ).Note that s(θ) is continuous-time and bandlimited to [−0.5, 0.5].
Figure 1 shows the system model used in the paper.After the Tx generates the RF signal, it propagates to Rx i with attenuation a i (which is considered known) and propagation time D pi , where i ∈ {1, 2} is the channel index.The receivers are connected to a fusion center via digital links.The Tx local oscillator phase misalignment at the sample s(0) is ϕ LOTX .In Fig. 1 the signals at the receivers are written relative to the Rx time axis.The local oscillator frequency offset in Rx i relative to ω c = 2πν c is ω i in [rad/s].Its phase misalignment at the beginning of the (Rx) observation interval is ϕ LOi .Although in theory s(t) is periodic, in practice the Tx sends only a finite number of periods, each corresponding to the sequence q, in Fig. 2 denoted by −2, −1, 0, 1.The figure also shows the propagation delays and the misalignment of the Tx time axis relative to the Rx time axis, denoted D 0 in seconds.The Rx time axis origin is determined by the beginning of the observation interval.The Rx system chooses the observation interval arbitrarily as long as the signal exists in both channels for the entire duration of the interval.Then, the origin of the Rx time axis is at the beginning of that observation interval, see Fig. 2. Therefore, the time shift of the signal s(t) in Rx i relative to the Rx time axis is T. Throughout the rest of the paper, we will use time that is normalized with the sampling interval, 1/ν s , and the quantities defined as and τ i = D i ν s in samples, for i ∈ {1, 2}.Then, the samples that the receivers acquire in T, after scaling by 1/a i , are: where is the normalized circular carrier frequency; n 1 and n 2 are independent circularcomplex AWGN processes with variance σ 2 ; for simplicity, it is assumed that the channel SNRs are known and equal; signal attenuations, a i , are real-valued and included in σ are the relative phase offset, frequency offset and time delay (TDOA), respectively.Note that the argument of s(•) need not be an integer.Given the samples in ( 3)-( 4), the goal is to estimate ϕ, Ω, and τ, or a subset of them, depending on which of them are unknown.
The discrete time matrix form of the signal model ( 3)-( 4), similar to [8], is: where where exp is applied to a matrix element by element, and Diag is the diagonal matrix with the elements of the vector argument.Note that D τ contains two terms, exp(−jΩ c τ), which models CPOA, and exp(−j2π/N • mτ), which models the complex envelope time shift.Even though CPOAs, i.e.Ω c τ 1 and Ω c τ 2 , are not elements of the unknown parameter vector, they are analyzed implicitly through propagation times, i.e. τ 1 and τ 2 .

Cramer-Rao Bounds
Since the given signal model differs from the existing ones, CRBs for estimation of unknown parameters are derived in this section for the three previously defined cases.
In case 1, all six parameters are unknown and the vector containing them is α 1 = (ϕ 1 , ϕ, Ω 1 , Ω, τ 1 , τ).The Fisher Information Matrix (FIM), [15], which is symmetric, is where the loglikelihood function, G, is where The expressions for the first and the second partial derivatives of G are given in the appendix.Since all the terms in the derivatives except x 1 (k) and x 2 (k) are deterministic, it follows To calculate the CRBs for appropriate parameters, the FIM, I 1 , needs to be inverted.Since the rows (and columns) of I 1 and I −1 1 correspond to the unknown parameters in the order in which they appear in α 1 , the CRBs for ϕ, Ω and τ are respectively The expressions ( 26)-( 28) are evaluated numerically.
In case 2, time shifts τ 1 and τ are known.The unknown parameter vector is α 2 = (ϕ 1 , ϕ, Ω 1 , Ω), so the last two rows and columns of I 1 are removed to obtain I 2 .If we assume the energy of s(k) is roughly equally distributed in time domain, the approximate CRBs are then where is the signal-to-noise ratio.These CRBs do not depend on Ω c , and are twice the CRBs in [15] (eq.15.72), for a cisoid in a single complex-AWGN channel.

Estimation Algorithms
This section introduces two algorithms, A and B, where A does not rely on knowing the sequence q, and B does.
Let us define block matrices Then the ML criterion reduces to finding the minimum of X − Qq 2 , where • is the Frobenius norm.If the sequence q is considered unknown, the estimate of q that minimizes the norm is q = Q † X, where and the fact that multiplying a signal with Q 1 or Q 2 does not change its energy, the ML reduces to φ, Ω1 , Ω, τ = arg max where Note that J 1 resembles cross-correlation, and that the criterion (33) resembles maximizing the energy at the output of a two-channel beamformer [16].Also note that J 1 does not depend on ϕ 1 or τ 1 , but depends on Ω 1 .Recall that we are only interested in estimating relative parameters, ϕ, Ω, τ.By further numerical optimization, by eliminating the search over ϕ, but only if ϕ is unknown, we get Ω1 , Ω, τ = arg max If some of the parameters are known, their true values are input into either (34) (if ϕ is known) or ( 35)-(36) (if ϕ is unknown) and they are not searched for.This concludes the derivation of Algorithm A. The advantages of this algorithm are: that it can be applied to a non-cooperative transmitter (receivers do not need to know the transmitted sequence), and the local oscillator and D/A converter in it need not be aligned to each other (which is the case with the majority of commercial transmitters).However, this algorithm has to be executed in the fusion center.An algorithm does not have to use all of the available information in the model for which the CRB is derived.This is why we compare the performance of Algorithm A to the previous CRB (for known signals).
The algorithm is implemented as follows.It performs a 3D search over a grid defined as the Cartesian product of the elements of column vectors Ω 1 , Ω and τ.The vectors contain N Ω 1 , N Ω and N τ values of interest of Ω 1 , Ω and τ, respectively.The inputs are x 1 , x 2 , N, Ω c , and the search grid.The steps of the algorithm are: Φ( j, k, l) = arg(s p ) 10.

end for
11.
end for 12. end for 13. ( j 0 , k 0 , l 0 ) = arg max f 14.Ω1 , Ω, τ, φ = Ω 1 ( j 0 ) , Ω (k 0 ) , τ(l 0 ), Φ( j 0 , k 0 , l 0 ) 15. return Ω1 , Ω, τ, φ Starting with the ML criterion, but for a known sequence q, we derive Algorithm B. Algorithm B is similar to A, except it is calculated for each pair (x i , q), i ∈ {1, 2}, and, instead of compensating for offsets in x i , the offsets are induced into q (mind the order of the operations).Finally, the difference of the estimates is calculated.More formally: φi , Ωi , τi = arg max It can be shown that this is equivalent to the joint ML search for (ϕ 1 , Ω 1 , τ 1 , ϕ, Ω, τ).A numerical optimization similar to the one applied to (34) can be applied to (37) as well.One advantage of Algorithm B is that signal processing can be distributed among the receiving channels, whereas Algorithm A is centralized.However, for Algorithm B to work correctly, the carrier has to have an integer number of cycles during one sequence period, i.e. ν c T has to be an integer.
The algorithm performs a 2D search over a grid defined as the Cartesian product of the elements of column vectors Ω i and τ i , which have N Ω and N τ elements, respectively.In addition to the input arguments of Algorithm A, q is an input for B as well.The steps of the algorithm are: for k = 1 to N Ω do 5.

end for
10.

end for
11.
Ωi , τi , φi = Ω i ( j 0 ) , τ i (k 0 ), Φi ( j 0 , k 0 ) 13. end for 14.Ω, τ, φ = Ω2 − Ω1 , τ2 − τ1 , φ2 − φ1 15. return Ω, τ, φ Let us discuss the ambiguity problem which arises from using CPOAs for estimation.Figure 3 a) shows a 3D slice representation of the criterion function J 1 for the true value of Ω 1 .Here the ambiguity manifests itself as multiple lobes along the τ-axis.Figures 3 b), c) and d) show the criterion function for estimating τ in case 3, for two different carriers and the sum of their criterion functions as a way to suppress ambiguity.The cross-correlation of a bandpass signal in [12] has a similar shape to b) and c).Our simulation results, not shown here, feature an SNR threshold for meansquare error (MSE) when the ambiguity becomes dominant.All the results shown in the next section are obtained under the assumption that the ambiguity is resolved.Development of an ambiguity-robust algorithm will be a part of future work.

Numerical Results
In order to examine the effects of taking into account CPOA on the estimation with actual algorithms, we performed extensive Monte Carlo simulations.We used an adaptive grid to speed up the execution of the algorithms.In all search dimensions, the grid was narrowed so that it spanned only five points in each dimension in any given iteration.If the argument of the maximum of the criterion function was at an edge of the grid, calculation was repeated after shifting the grid by three points in the direction of that edge.This was repeated until the argument of the maximum was not at any edge of the grid, which ended the current estimation.The number of simulation runs at each SNR value was 8192.The spacings in the grid were determined by observing the convergence of results with their change.The grid spanned a narrow interval around expected parameter values, in order to achieve a satisfactory precision at a moderate computation complexity.As a consequence, we included only one criterion function lobe into the grid, so the ambiguity problem due to using CPOAs did not influence the results.All of the results shown in figures in this section were generated using the first of the modulatable orthogonal sequences proposed in [17] of different lengths N. It is chirp-like so the parameter estimation errors are correlated, see [8].

Figures 4 and 5 present MSE curves for algorithms
A and B versus SNR for frequency offset Ω and relative time delay τ estimation, respectively.All 3 parameters ϕ, Ω and τ are considered unknown (case 1).If not otherwise stated, we consider that ϕ and ϕ 1 are both either known or unknown, and similarly for pairs (Ω, Ω 1 ) and (τ, τ 1 ).The normalized   circular carrier frequency, Ω c , was set to 2π100, whereas N ∈ {256, 1024, 4096}.Appropriate CRB curves are displayed as well.Algorithm B closely follows the CRB over the entire range of observed SNRs and it outperforms Algorithm A at lower SNRs, which deviates in this region due to not relying on knowledge of the transmitted signal.The results do not depend on Ω c (based on numerical results).The information embedded in CPOA is lost because the initial phases are unknown, so there are no estimation accuracy gains in this case.As expected, the MSEs decrease with the increase of N and SNR.
Figure 6 a) shows the dependance of ϕ MSE on Ω c , in case 1 (in which τ is unknown) for N = 4096, and Ω c /(2π) ∈ {10, 100, 1000}.These values for Ω c roughly correspond to UWB, LTE and GSM cases, respectively.Note that, since CRBs do not take into account that ϕ ∈ [−π, π), CRB curves surpass the corresponding MSE curves.A small error in TDOA estimation can produce a large error in the initial phase estimate due to the exp(−jΩ c τ) phase term and this error increases linearly with Ω c .For comparison, the MSE curves for case 2 (time delays are known) are displayed as well.We show results only for Ω c = 2π100 as these curves are not dependent on Ω c , see (29) (in this case the term exp(−jΩ c τ) induces no errors).This implies that knowing the time delays as accurately as possible is crucial for phase offset estimation.For comparison, we also present the MSE curves for case 1, for Ω c = 2π100, as these curves do not depend on carrier frequency.The results for case 3 show that τ RMS error is by 1 to 3 orders of magnitude lower than 1/Ω c , which presents a significant reduction of TDOA RMS error compared to the estimation techniques not using the carrier phase, which have an RMS error 1 to 3 orders of magnitude lower than 1/B, similar to our methods in case 1 (see Fig. 5).The main reason for this accuracy gain is the correlation between τ and ϕ due to the exp(−jΩ c τ) term, so this benefit from taking into account carrier phase information is only available when ϕ is known.
Let the carrier frequency be ν c = 2 GHz, the signal bandwidth B = 5 MHz, the signal-to-noise ratio 15 dB, and N = 1024.This implies that Ω c = 2π400.In case 3 Algorithm A then gives τ RMS error of 2.18 × 10 −6 samples (therefore it implicitly shows the dependence on B), which is 4 × 10 −13 s, which in propagation distance equals to 120 µm.This is 0.0008λ c , where λ c is the carrier wavelength.
Instead of a chirp-like sequence, realizations of Gaussian random processes have also been tested.On average, the results were only slightly better.Based on the estimated mean-to-standard-deviation ratios from simulations (not shown here), the estimators are unbiased.

Conclusion
In this paper we have analyzed the impacts of taking into account CPOA on the time, phase, and frequency shifts estimation accuracy in a distributed frequency-unlocked timesynchronized 2-channel receiving system.CRBs and numerical results were used to determine under which conditions gains were achievable and to quantify them.If all three types of parameters are unknown, CPOA creates no gains in the accuracy.On the contrary, the initial phase error is increased.In the case of unknown phase and frequency offsets, taking into account CPOA does not influence the estimation accuracy.In the case of unknown TDOA only, CPOA drastically improves the accuracy, so that an RMS error 1 to 3 orders of magnitude lower than the inverse of the carrier frequency can be achieved, if the ambiguity problem is solved.If CPOA were not used in this case, the achievable error would be 1 to 3 orders of magnitude lower than the inverse of the bandwidth at best.The proposed algorithms have MSEs that approach the CRBs.We believe that this analysis shows how far the estimation accuracies can be pushed by exploiting CPOA, and which conditions have to be met to have TDOA estimation accuracy good enough to be used, for example, for beamforming by a distributed antenna array.

Fig. 2 .
Fig. 2. Temporal relationships between signals and Tx and Rx axes.oscillators, IQ demodulated and sampled at the Nyquist rate, ν s = B.Each Rx acquires N samples in a common observation interval of length T = N/ν s whose beginning is chosen independently from the start of the transmission.

Fig. 3 .
Fig. 3. Criterion functions showing the ambiguity problem and a possible solution.

Fig. 4 .
Fig. 4. CRB and MSE for frequency offset estimation when all the parameters are unknown.

Fig. 5 .
Fig. 5. CRB and MSE for TDOA estimation when all the parameters are unknown.

Fig. 6 .
Fig. 6. a) Dependence of phase offset CRB and MSE on carrier frequency.b) Dependence of TDOA CRB and MSE on carrier frequency.

Figure 6 b
Figure 6 b) shows τ MSE for case 3, for N = 4096, and Ω c /(2π) ∈ {10, 100, 1000}.It can be seen that the MSE decreases with Ω c .For comparison, we also present the MSE curves for case 1, for Ω c = 2π100, as these curves do not depend on carrier frequency.The results for case 3 show that τ RMS error is by 1 to 3 orders of magnitude lower than 1/Ω c , which presents a significant reduction of TDOA RMS error compared to the estimation techniques not using the carrier phase, which have an RMS error 1 to 3 orders of magnitude lower than 1/B, similar to our methods in case 1 (see Fig.5).The main reason for this accuracy gain is the correlation between τ and ϕ due to the exp(−jΩ c τ) term, so this benefit from taking into account carrier phase information is only available when ϕ is known.