# CSP-Based Time-Difference-of-Arrival Estimation

Let’s discuss an application of cyclostationary signal processing (CSP): time-delay estimation. The idea is that sampled data is available from two antennas (sensors), and there is a common signal component in each data set. The signal component in one data set is the time-delayed or time-advanced version of the component in the other set. This can happen when a plane-wave radio frequency (RF) signal propagates and impinges on the two antennas. In such a case, the RF signal arrives at the sensors with a time difference proportional to the distance between the sensors along the direction of propagation, and so the time-delay estimation is also commonly referred to as time-difference-of-arrival (TDOA) estimation.

Consider the diagram shown to the right. A distant transmitter emits a signal that is well-modeled as a plane wave once it reaches our two receivers. An individual wavefront of the signal arrives at the two sensors at different times.

The line segment AB is perpendicular to the direction of propagation for the RF signal. The angle $\theta$ is called the angle of arrival (AOA). If we could estimate the AOA, we can tell the direction from which the signal arrives, which could be useful in a variety of settings. Since the triangle ABC is a right triangle, we have

$\displaystyle \cos (\theta) = \frac{x}{d}. \hfill (1)$

When $\theta = 0$, the wavefronts first strike receiver 2, then must propagate over $x$ meters before striking receiver 1. On the other hand, when $\theta = 90^\circ$, each wavefront strikes the two receivers simultaneously. In the former case, the TDOA is maximum, and in the latter it is zero. The TDOA can be negative too, so that $180^\circ$ azimuthal degrees can be determined by estimating the TDOA.

In general, the wavefront must traverse $x$ meters between striking receiver 2 and striking receiver 1,

$\displaystyle x = d \cos(\theta). \hfill (2)$

Assuming the speed of propagation is $c$ meters/sec, the TDOA is given by

$\displaystyle D = \frac{x}{c} = \frac{d\cos{\theta}}{c} \mbox{\rm \ \ seconds}. \hfill (3)$

In this post I’ll review several methods of TDOA estimation, some of which employ CSP and some of which do not. We’ll see some of the advantages and disadvantages of the various classes of methods through inspection, simulation, and application to collected data.

### Signal Models

Throughout this post, the data models of interest correspond to the following physical situation:

$\displaystyle x(t) = s(t) + n_x(t), \hfill (4)$

$\displaystyle y(t) = s(t-D) + n_y(t), \hfill (5)$

where $x(t)$ and $y(t)$ are the observed data from two physically separated sensors (antennas and their receivers), $D$ is the time difference of arrival (TDOA), and $n_x(t)$ and $n_y(t)$ are independent receiver noises for the two sensors. The two receivers are assumed to employ a common timing reference so that the only contribution to $D$ is the propagation delay.

This model neglects several potentially important aspects of real-world collected data sets, such as the Doppler effect, time variation of $D$, cochannel interference, and phase mismatch between the signal components ($s(t)$ components) in each data set. That is, a more complete model might look like

$\displaystyle x(t) = As(t)e^{i2\pi f_{d_x} t} + n_x(t) + I_x(t), \hfill (6)$

$\displaystyle y(t) = Bs(t-D(t))e^{i2\pi f_{d_y} t} + n_y(t) + I_y(t), \hfill (7)$

where $A$ and $B$ are complex scalars (numbers), $f_{d_x}$ and $f_{d_y}$ are the Doppler shifts for the two sensors, and $I_x(t)$ and $I_y(t)$ are one or more cochannel interferers, each of which has its own TDOA, Dopplers, etc.

The noises $n_x(t)$ and $n_y(t)$ are assumed to be stationary, or at least do not possess any cycle frequencies of $s(t)$ that we wish to exploit. They are also assumed to be statistically independent and have mean values of zero.

We have developed multiple TDOA estimation algorithms using the simpler model, as have others, and then we end up applying them to the more complex signal models. This is due to the general degree of difficulty with mathematically deriving algorithms using the more complex models.

### Algorithms

#### Cross Correlation

A simple, and often effective, approach to estimating $D$ is to estimate the cross correlation between $x(t)$ and $y(t)$, and find the location of its maximum. That location is then the TDOA estimate $\hat{D}$.

The cross correlation method works because the autocorrelation for any signal achieves its maximum at a delay (correlation lag) of zero. (It may achieve this same maximum at other delays as well, but will definitely achieve it at $\tau = 0$. Consider $x(t) = \sin(t)$.)

The cross correlation between $x(t)$ and $y(t)$ is

$\displaystyle R_{xy}(t, \tau) = E\left[ x(t+\tau/2) y^*(t-\tau/2)\right] \hfill (8)$

or

$\displaystyle R_{xy}(t, \tau) = E\left[ s(t+\tau/2)s^*(t-\tau/2-D) + s(t+\tau/2)n_y^*(t-\tau/2) \right.$

$\displaystyle \left. + s^*(t-\tau/2-D)n_x(t+\tau/2) + n_x(t+\tau/2)n_y^*(t-\tau/2) \right]. \hfill (9)$

The expectation can be distributed to the four terms. The final three terms are all zero. The cross correlations between the signal and noise are zero because they are independent, and the noises have zero mean values. The final term, “noise times noise,” is zero because the two noises are independent, and have zero mean values. So we end up with

$\displaystyle R_{xy}(t, \tau) = E\left[ s(t+\tau/2)s^*(t-\tau/2-D) \right]. \hfill (10)$

Now, we’d like to express this cross correlation in terms of a function that we know how to maximize, such as the autocorrelation for $s(t)$. But the form isn’t yet familiar; we have defined our functions with symmetrical lag variables $\tau/2$ and $-\tau/2$,

$\displaystyle R_s(t, \tau) = E \left[ s(t + \tau/2) s^*(t-\tau/2) \right].\hfill (11)$

We can force our desired form in the following way. Consider the autocorrelation for two arbitrary time instants $t_1$ and $t_2$:

$\displaystyle R_x(t_1, t_2) = E\left[ x(t_1) x^*(t_2)\right]. \hfill (12)$

Let’s define two new variables that are the difference and the average value of the two time instants $t_1$ and $t_2$,

$\displaystyle u = t_1 - t_2 \hfill (13)$

$\displaystyle v = \frac{t_1 + t_2}{2}. \hfill (14)$

Then $t_1$ and $t_2$ can be expressed in terms of $u$ and $v$ as follows

$\displaystyle t_1 = v + u/2 \hfill (15)$

$\displaystyle t_2 = v - u/2. \hfill (16)$

Therefore the generic autocorrelation for $x(t)$ is re-expressed as

$\displaystyle R_x(t_1, t_2) = R_x(v+u/2, v-u/2) = E[x(v+u/2)x^*(v-u/2)], \hfill (17)$

as desired. For our cross-correlation problem, we identify

$t_1 = t + \tau/2 \hfill (18)$

$t_2 = t - \tau/2 - D. \hfill (19)$

$u = t + \tau/2 - (t - \tau/2 - D) = \tau + D \hfill (20)$

$v = \frac{t + \tau/2 + t - \tau/2 - D}{2} = t - D/2. \hfill (21)$

So we have

$\displaystyle R_{xy}(t, \tau) = E\left[ s(v+u/2) s^*(v-u/2) \right]. \hfill (22)$

$\displaystyle R_{xy}(t, \tau) = R_s(v, u) = \sum_\alpha R_s^\alpha (u) e^{i2\pi \alpha v} \hfill (23)$

$\displaystyle R_{xy}(t, \tau) = \sum_\alpha R_s^\alpha (\tau + D) e^{i2\pi \alpha (t-D/2)}. \hfill (24)$

For the cross-correlation method, let’s focus on the $\alpha = 0$ term here, or

$\displaystyle R_{xy}^0 (\tau) = R_s^0 (\tau + D) e^{0} = R_s^0 (\tau + D). \hfill (25)$

Since the autocorrelation function for $s(t)$ peaks at $\tau = 0$, the cross correlation with $\alpha = 0$ between $x(t)$ and $y(t)$ peaks at $\tau = -D$. (Don’t worry about the sign of $D$ here, it is just an artifact of which receiver we associate with $x(t)$ and which with $y(t)$.)

In practice, we simply form an estimate of the cross correlation and find the location of its maximum.

$\displaystyle \hat{D} = \mbox{\rm arg} \max_{\tau} \left| \hat{R}_{xy} (\tau) \right|. \hfill (26)$

#### Generalized Cross Correlation (GCC)

The cross correlation can be expressed as the inverse Fourier transform of the cross spectrum between $x(t)$ and $y(t)$, which is denoted by $S_{xy}^0 (f)$,

$\displaystyle R_{xy}^0(\tau) = \int S_{xy}(f) e^{i2\pi f \tau} \, df. \hfill (27)$

Notice that if there are frequency bands for which the cross spectrum $S_{xy}(f)$ is small relative to the noise power, the contribution to the inverse transform can be harmful and cause TDOA estimation error. Thus, we may wish to weight the cross spectrum prior to inverse transformation, such as with a general weighting function $W(f)$,

$\displaystyle \hat{D} = \mbox{\rm arg} \max_{D} \int W(f) S_{xy}(f) e^{i2\pi f \tau} \, df . \hfill (28)$

This form can also be derived by modifying our problem statement to include arbitrary linear time-invariant filtering applied to each of $x(t)$ and $y(t)$ prior to cross-correlation-based TDOA estimation. A variety of weighting factors, and their motivations and properties, are discussed in Knapp/Carter’s paper on Generalized Cross Correlation (GCC) (The Literature [R52]).

In this post, we will use $W(f) = S_s^0(f)$, which we will approximate as $W(f) = \hat{S}_x^0(f)$. This is essentially the Eckhart filter ([R52]) when the receiver noises are white and possess the same variances. Other helpful references on the GCC class of estimators were published by Knapp, Carter, Hero, and Kirlin [R53]-[R60].

#### Cyclic Cross Correlation (CCC)

Recall from our development for the cross-correlation method that the time-varying cross correlation between $x(t)$ and $y(t)$ is given by

$\displaystyle R_{xy}(t, \tau) = \sum_\alpha R_s^\alpha (\tau + D) e^{i2\pi \alpha (t-D/2)}. \hfill (29)$

An individual $\alpha \neq 0$ cyclic cross correlation term in this cross correlation is given by

$\displaystyle R_{xy}^\alpha(\tau) = R_s^\alpha(\tau + D) e^{i2\pi \alpha (D/2)}. \hfill (30)$

If the underlying cyclic autocorrelation $R_s^\alpha (\tau)$ achieves its maximum at $\tau = 0$, then the magnitude of the cyclic cross correlation will achieve its maximum at $\tau = -D$. Therefore, the cyclic cross correlation can be used to form an estimate of the delay $D$ in a manner similar to the conventional cross correlation,

$\displaystyle \hat{D} = \mbox{\rm arg} \max_D \left| \hat{R}_{xy}^\alpha (\tau) \right |. \hfill (31)$

For some signals (textbook rectangular-pulse!), the cyclic autocorrelation does not peak at $\tau = 0$, but may exhibit a pair of peaks whose center is indeed at $\tau =0$. For such signals, the search for a single peak can be replaced by a search for a pair of peaks. The TDOA estimate is then the midpoint of the locations of the two peaks.

#### Spectral Coherence Alignment (SPECCOA)

The SPECCOA method arises from a least-squares fitting problem involving the relationship between the cyclic cross correlation and the cyclic autocorrelation. Recall that this relation is given by

$\displaystyle R_{xy}^\alpha(\tau) = R_x^\alpha(\tau + D) e^{i2\pi \alpha (D/2)}. \hfill (32)$

If we perform a least-squares minimization over the error given by the absolute value of the differences between the left and right sides of that equation, we end up with the following estimator

$\displaystyle \hat{D} = \mbox{\rm arg} \max_D \left| \int \hat{S}_{xy}^\alpha (f) \hat{S}_x^\alpha(f)^* e^{i 2 \pi f D} \, df \right|. \hfill (33)$

The SPECCOA estimator is derived in The Literature [R48]-[R49].

#### Spectral Correlation Ratio (SPECCORR)

By taking a system-identification approach, Gardner and Chen (The Literature [R51]) derive yet another cyclostationarity-exploiting TDOA estimator that inverse transforms the ratio of the cyclic cross spectral correlation function and the cyclic spectral correlation function,

$\displaystyle \hat{D} = \mbox{\rm arg} \max_D \left| \int \frac{\hat{S}_{xy}^\alpha (f)} {\hat{S}_x^\alpha(f)} e^{i 2 \pi f D} \, df \right|. \hfill (34)$

#### Spectral Coherence Nulling (SPECCON)

Consider the signal $z(t) = x(t) - y(t-\Delta)$, where $x(t)$ and $y(t)$ are our usual noisy CS signals (4)-(5). If the delay/advance parameter $\Delta$ is equal to $-D$, then the spectral correlation at cycle frequency $\alpha$ for $z(t)$ will be zero, because the spectral correlation for the $s(t)$ components in $x(t)$ and $y(t-(-D)) = s(t + D - D) + n_y(t)$ will be identical. So we can consider minimizing the spectral correlation for $z(t)$ over the parameter $\Delta$. This leads to the SPECCON method (The Literature [R48]).

SPECCON is given by

$\displaystyle \hat{D} = \mbox{\rm arg} \min_\Delta \Re \left\{\int \left[ \vphantom{\int} S_x^\alpha(f) S_y^\alpha(f)^* e^{-i 2 \pi \alpha \Delta} + S_{yx}^\alpha (f)^* S_{xy}^\alpha(f) e^{-i4\pi f \Delta} \right. \right.$

$\displaystyle - (S_x^\alpha(f)^* S_{xy}^\alpha(f) + S_y^\alpha(f)S_{yx}^\alpha(f)^*)e^{-i2\pi(f-\alpha/2)\Delta}$

$\displaystyle \left. \left. \vphantom{\int} - (S_x^\alpha(f) S_{yx}^\alpha(f)^* + S_y^\alpha(f)^*S_{xy}^\alpha(f))e^{-i2\pi(f+\alpha/2)\Delta} \right] \, df \right\}. \hfill (35)$

#### Generalized SPECCOA (GSPECCOA)

GSPECCOA arises from considering all auto- and cross-SCF functions for the two sensor time-series, and jointly fitting them to an ideal SCF for the signal of interest (The Literature [R50]). The resultant algorithm can be expressed as

$\displaystyle \hat{D} = \mbox{\rm arg} \max_D \Re \left\{ \int \hat{S}_{xx}^\alpha (f) \hat{S}_{xy}^\alpha(f)^* e^{i2\pi (f-\alpha/2)D} \, df + \int \hat{S}_{yx}^\alpha(f) \hat{S}_{xx}^\alpha(f)^* e^{i2\pi(f+\alpha/2)D} \, df \right.$

$\displaystyle + \int \hat{S}_{yx}^\alpha(f) \hat{S}_{yy}^\alpha(f)^* e^{i2\pi (f-\alpha/2)D} \, df + \int \hat{S}_{yy}^\alpha(f) \hat{S}_{yx}^\alpha(f)^* e^{i2\pi (f+\alpha/2)D} \, df$

$\displaystyle \left. + \int \hat{S}_{yx}^\alpha(f) \hat{S}_{xy}^\alpha(f)^* e^{i4\pi f D} \, df + \int \hat{S}_{yy}^\alpha(f) \hat{S}_{xx}^\alpha(f)^* e^{i2\pi \alpha D} \, df \right\} . \hfill (36)$

#### Approximate Maximum Likelihood (AML)

I have attempted to derive the weak-signal maximum-likelihood TDOA estimator for the signal model involving gain-and-phase mismatched receivers ((6)-(7) with $I_x(t) = I_y(t) = 0$). The mathematics is rather complex, and tedious, so I’m just going to say that I’m including two approximations to the ML estimator that I call AML1 and AML5. They are similar to GSPECCOA.

### Examples

#### The Case of No Interference

Let’s consider a simulated textbook BPSK signal with square-root raised-cosine pulses. The excess bandwidth parameter (roll-off) is 1.0, the bit rate is $1/T_0 = 1/10$ (10 samples per signalling interval), and the carrier frequency is $f_c = 0.05$. So this signal has five cycle frequencies that we can exploit: the two non-conjugate cycle frequencies of 0 and 1/10, and the three conjugate cycle frequencies of $2f_c = 0.1$, $2f_c - 1/T_0 = 0.0$, and $2_fc + 1/T_0 = 0.2$. A small amount of noise is added so that the inband SNR is approximately 20 dB. The TDOA is set at 10 samples.

In all cases shown here, 65536 samples are processed. All spectral correlation functions are estimated using the frequency-smoothing method. All cross spectral correlations are not smoothed (more on that below), and all auto spectral correlation functions are smoothed with a rectangular smoothing window with width equal to 500 frequency bins. This means the smoothing window width is equal to 500/65536 = 0.0076 (normalized) Hz.

For each simulated algorithm, we plot the TDOA objective function. That is, we plot the function of delay that is either maximized or minimized by the algorithm to find the TDOA estimate. First, consider the conventional case of a non-conjugate cycle frequency of zero:

(Here SPECCORR is coincident with Cross Correlation.) All the algorithms objective functions peak at the correct delay of 10 samples (the true signal TDOA value). The next graph shows the objective functions for the case of the non-conjugate cycle frequency of $\alpha = 1/T_0 = 1/10$:

Once again, all the objective functions peak at the correct TDOA. The shapes of the functions are different from each other, and are also different from those for $\alpha = 0$. Each algorithm provides a slightly different trade-off between the width near the true TDOA and the height at points distant from the TDOA. This is akin to a resolution/sidelobe type of trade-off that is common in spectral analysis and other signal-processing settings.

Here are the objective functions for the three conjugate cycle frequencies:

Several of the algorithms have multi-cycle extensions, where multiple cycle frequencies are exploited simultaneously. Here are the objective functions for those multi-cycle algorithms when they exploit all five cycle frequencies:

#### The Case of Cochannel Interference

Now let’s add a cochannel interferer to our BPSK signal. The interferer will also be a BPSK signal, with slightly different parameters, but will completely overlap the original BPSK signal in both time and frequency–it cannot be removed by any sort of linear time-invariant filtering or time-gating operation.

The interferer has the same pulse function as the BPSK signal of interest, a bit rate of $1/T_1 = 1/9$ and a carrier frequency of $2f_c = 0.045$. These choices imply that the two signals do not share any cycle frequencies except for the (trivial) non-conjugate cycle frequency of zero. The interferer power and the original signal power are the same. This all results in the PSDs shown here:

The interferer’s TDOA is 17 samples. Let’s go through the same sequence of plots that we went through for the case of an interference-free signal. First, the conventional methods and the CSP methods that use only the non-conjugate cycle frequency of zero:

Several of the methods still peak at the TDOA of the BPSK signal, but each also peaks for the TDOA of the interferer. However, methods such as GSPECCOA peak at neither TDOA, but somewhere between the two TDOAs.

Without some further prior information, it would be impossible to reliably associate any of the local maxima with the TDOA for the signal–which peak belongs to which signal?

The next example corresponds to the non-conjugate cycle frequency of 0.1, which is the bit rate for the signal of interest, but is not a cycle frequency for the interferer:

Here are the objective functions for the conjugate cycle frequencies for the BPSK signal of interest:

And here are the objective functions for the case of combining the information in all five cycle frequencies:

All of the methods peak at the correct TDOA, although they also show significant values near the TDOA of the interferer. This is largely due to the inclusion of the non-conjugate cycle frequency of zero, which is sensitive to the presence of all signals. If we eliminate this cycle frequency, we obtain the following objective functions:

Here the SPECCOA and AML5 algorithms do not appear to reflect the presence of the interferer at all.

Finally, if we reverse the roles of the two modulated signals, and use the four cycle frequencies for the interferer instead of those for the original BPSK signal, we obtain the following set of objective functions:

Now the functions peak at the TDOA of the interferer, as expected. This is another illustration of the signal selectivity property of cyclostationary parameters.

#### Captured-Data Example

To reassure us that the various CS-exploiting methods actually work on real-world data, let’s apply them to a captured CDMA-EVDO signal. We process $1.3$ seconds of a downlink EVDO signal sampled at $12.5$ MHz. This signal possesses non-conjugate cycle frequencies that are harmonics of $1.230$ MHz. We induce a TDOA of $10.333$ samples by applying a delay operation to the captured data. The $1.3$-second data record is processed by extracting successive blocks with length $1.3$ ms and applying GSPECCOA for various cycle-frequency selections.

Here are the histograms of the obtained TDOA estimates:

In the histogram plot, the histogram-bin locations are constrained to lie within $[9.6 10.6]$ samples, so all TDOA estimates less than $9.6$ samples are lumped into the bin at $9.6$, and similarly for $10.6$. Clearly the cycle-frequency set that uses all five harmonics provides nearly perfect performance, producing the nearest estimate to $10.3333$ samples (thatis, $10$ samples) that it can for almost all trials.

### Possible Pitfalls

#### Temporal Resolution

I’ve presented the algorithms’ objective functions in terms of samples, and I used a simulated example where the true TDOA is equal to an integer number of samples ($10$). In practice, the TDOA is simply a function of the propagation speed and angle of arrival, and so is continuous, not discrete. But we want to work with sampled data, so we are stuck with discrete-valued TDOA estimates.

The discrete nature of the TDOA estimate in digital signal processing leads to unavoidable errors when applied to actual captured data. Perhaps even more importantly, it leads to the possibility of erroneous performance analyses when Monte Carlo simulations are used. In particular, if the TDOA is selected to be an integer number of samples, the various algorithms can yield perfect estimates (zero errors) over many many Monte Carlo trials.

One way of avoiding this is to interpolate the data prior to applying the TDOA estimate.  It could be interpolated such that the minimum non-zero error is small relative to the desired error in some application. At least in this case, the relevant performance can be captured, and the algorithms correctly ranked.

When interpolation is applied to the EVDO example above, we obtain the following histograms:

As the level of interpolation increases, the TDOA estimates inch closer to the true value of $10.3333$ (which does not lie on any interpolated point, since the interpolation factors are multiples of two).

#### Spectral Smoothing

It is tempting to apply lots of smoothing to the auto and cross spectral correlation functions in a given CSP algorithm, so as to decrease the variance of the estimate, and thereby improve the variability of the TDOA estimates. But there is a catch. We can apply plenty of smoothing to any auto spectral correlation estimate, but the ideal cross spectral correlation functions contain a sinusoidal factor whose period in the frequency variable depends on the unknown TDOA $D$. Thus, these functions can easily be oversmoothed, causing the destructive interference between periods of the sinusoidal factor. I typically smooth the auto estimates and do not smooth the cross estimates at all.

I’ll leave it to you, dear reader, to calculate the cross spectral correlation between $x(t)$ and $y(t)$ in order to see the problematic complex sinusoidal term.