CSP-Based Time-Difference-of-Arrival Estimation – Cyclostationary Signal Processing

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.

tdoa_physical_setup — Figure 1. Illustration of the geometric relationship between a transmitter and two receivers in the context of time-delay estimation (or time-difference-of-arrival estimation).

Consider the diagram shown in Figure 1. 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=d$ 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 captured data. Consider this post as a starting point to a study or development effort rather than as a definitive performance characterization.

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 captured 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, Doppler shift, etc.

The noises $n_x(t)$ and $n_y(t)$ are assumed to be stationary, or at least do not possess any of the 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 data adhering to the more complex signal models. This is due to the generally high 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)$

$\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)$

This leads to

$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 (no CSP):

single_alpha_0.0 — Figure 2. Time-delay objective functions for the case of a single non-conjugate cycle frequency that is equal to zero. The GCC and Cross Correlation methods are conventional (inherently non-CSP) methods.

(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$ :

single_alpha_0.1 — Figure 3. Time-delay objective functions for the case of a single non-conjugate cycle frequency that is equal to the bit rate of the BPSK signal ( $0.1$ ). Since we are exploiting cyclostationarity here, the conventional methods’ objective functions are not applicable and so are not shown.

Figure 3. Time-delay objective functions for the case of a single non-conjugate cycle frequency that is equal to the bit rate of the BPSK signal ( $0.1$ ). Since we are exploiting cyclostationarity here, the conventional methods’ objective functions are not applicable and so are not shown.

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:

single_beta_0.0 — Figure 4. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to zero. This cycle frequency corresponds to $2f_c - f_{bit}$ .

Figure 4. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to zero. This cycle frequency corresponds to $2f_c - f_{bit}$ .

single_beta_0.1 — Figure 5. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.1$ . This cycle frequency corresponds to $2f_c$ .

Figure 5. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.1$ . This cycle frequency corresponds to $2f_c$ .

single_beta_0.2.jpg — Figure 6. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.2$ . This cycle frequency corresponds to $2f_c + f_{bit}$ .

Figure 6. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.2$ . This cycle frequency corresponds to $2f_c + f_{bit}$ .

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. If it could, we could identify each signal’s band, filter them to isolate them, and then apply the cross-correlation (non-CSP) methods and we’d be doing very well indeed.

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 $f_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:

psds_for_tdoa — Figure 8. Power spectrum estimates for the noisy received data, which is the combination of the two BPSK signals and noise, as well as for each noise-free BPSK signal for illustration purposes.

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 (stationary-signal processing):

single_interf_alpha_0.0_interf — Figure 9. Time-delay objective functions for the case of a single non-conjugate cycle frequency that is equal to zero as applied to the case of a BPSK signal and a cochannel BPSK interferer. The GCC and Cross Correlation methods are conventional (inherently non-CSP) methods.

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 BPSK signal of interest, but is not a cycle frequency for the BPSK interferer:

single_interf_alpha_0.1_interf — Figure 10. Time-delay objective functions for the case of a single non-conjugate cycle frequency that is equal to $0.1$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The bit rate of the BPSK signal of interest is equal to the employed cycle frequency ( $0.1$ ).

Figure 10. Time-delay objective functions for the case of a single non-conjugate cycle frequency that is equal to $0.1$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The bit rate of the BPSK signal of interest is equal to the employed cycle frequency ( $0.1$ ).

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

single_interf_beta_0.0_interf — Figure 11. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The cycle frequency is equal to $2f_c - f_{bit}$ where $f_c$ is the carrier frequency of the BPSK signal of interest and $f_{bit}$ is the bit rate of the BPSK signal of interest.

Figure 11. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The cycle frequency is equal to $2f_c - f_{bit}$ where $f_c$ is the carrier frequency of the BPSK signal of interest and $f_{bit}$ is the bit rate of the BPSK signal of interest.

single_interf_beta_0.1_interf — Figure 12. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.1$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The cycle frequency is equal to $2f_c$ where $f_c$ is the carrier frequency of the BPSK signal of interest.

Figure 12. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.1$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The cycle frequency is equal to $2f_c$ where $f_c$ is the carrier frequency of the BPSK signal of interest.

single_interf_beta_0.2_interf — Figure 13. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.2$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The cycle frequency is equal to $2f_c + f_{bit}$ where $f_c$ is the carrier frequency of the BPSK signal of interest and $f_{bit}$ is the bit rate of the BPSK signal of interest.

Figure 13. Time-delay objective functions for the case of a single conjugate cycle frequency that is equal to $0.2$ as applied to the case of a BPSK signal and a cochannel BPSK interferer. The cycle frequency is equal to $2f_c + f_{bit}$ where $f_c$ is the carrier frequency of the BPSK signal of interest and $f_{bit}$ is the bit rate of the BPSK signal of interest.

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

multi_interf — Figure 14. Time-delay objective functions for the case of simultaneous exploitation of all five of the BPSK signal’s cycle frequencies (two non-conjugate and three conjugate) as applied to the noisy combination of the BPSK signal and BPSK interferer. Only those time-delay estimation algorithms that inherently support multiple cycle frequencies are shown.

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:

just_four_interf — Figure 15. Time-delay objective functions for the case of simultaneous exploitation of four of the five cycle frequencies for the BPSK signal of interest (the non-conjugate cycle frequency of $0$ is omitted) as applied to the noisy combination of the BPSK signal and BPSK interferer. Only those time-delay estimation algorithms that inherently support multiple cycle frequencies are shown.

Figure 15. Time-delay objective functions for the case of simultaneous exploitation of four of the five cycle frequencies for the BPSK signal of interest (the non-conjugate cycle frequency of $0$ is omitted) as applied to the noisy combination of the BPSK signal and BPSK interferer. Only those time-delay estimation algorithms that inherently support multiple cycle frequencies are shown.

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:

just_four_snoi — Figure 16. Time-delay objective functions for the case of simultaneous exploitation of four of the five cycle frequencies for the BPSK interferer (the non-conjugate cycle frequency of $0$ is omitted) as applied to the noisy combination of the BPSK signal and BPSK interferer. Only those time-delay estimation algorithms that inherently support multiple cycle frequencies are shown.

Figure 16. Time-delay objective functions for the case of simultaneous exploitation of four of the five cycle frequencies for the BPSK interferer (the non-conjugate cycle frequency of $0$ is omitted) as applied to the noisy combination of the BPSK signal and BPSK interferer. Only those time-delay estimation algorithms that inherently support multiple cycle frequencies are shown.

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 Fourier-based 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:

tdoa_evdo — Figure 17. Histograms of GSPECCOA TDOA estimates as applied to a captured CDMA EVDO signal with artificially induced TDOA of $10.333$ samples. The dominant non-conjugate cycle frequencies exhibited by this signal are multiples of $1.2300$ MHz.

Figure 17. Histograms of GSPECCOA TDOA estimates as applied to a captured CDMA EVDO signal with artificially induced TDOA of $10.333$ samples. The dominant non-conjugate cycle frequencies exhibited by this signal are multiples of $1.2300$ MHz.

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 (that is, $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 actual TDOA is simply a function of the propagation speed and angle of arrival, and so is a continuous-valued variable, not a discrete one. 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. (My algorithm is perfect! Let’s publish!)

One way of avoiding this misleading result 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:

tdoa_evdo_interp — Figure 18. Histograms of GSPECCOA TDOA estimates as applied to a captured CDMA EVDO signal with artificially induced TDOA of $10.333$ samples. The dominant non-conjugate cycle frequencies exhibited by this signal are multiples of $1.2300$ MHz. Here the first five harmonics of $1.2300$ MHz are used in GSPECCOA. The two data records are interpolated (sample rate is increased) prior to application of GSPECCOA, causing the nearest objective-function point to get closer and closer to the true time delay.

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.

Author: Chad Spooner

I'm a signal processing researcher specializing in cyclostationary signal processing (CSP) for communication signals. I hope to use this blog to help others with their cyclo-projects and to learn more about how CSP is being used and extended worldwide. View all posts by Chad Spooner

2 thoughts on “CSP-Based Time-Difference-of-Arrival Estimation”

Andrey says:

March 18, 2021 at 1:02 pm

Hi! Thank you for your blog.
Could you explain in details or simple words CyclicMUSIC method for DOA estimation? Origin of this method is WILLIAM A. GARDNER, 1988, “Simplification of MUSIC and ESPRIT by Exploitation of Cyclostationarity”.

1. Chad Spooner says:
  
  March 20, 2021 at 11:54 am
  
  Hey Andrey! Thanks for visiting the CSP Blog and the comment.
  
  I have a draft post on cyclic MUSIC, but I’m not sure when I can get it completed. I have too many draft posts! And a day job that gets in the way of the CSP Blog.
  
  In MUSIC, the basic idea is that the singular vectors for the array-data correlation matrix can be decomposed into two subspaces: a signal subspace and a noise subspace. The signal subspace is spanned by the steering vectors corresponding to the signals impinging on the array. The steering vectors are the collection of array-response vectors corresponding to all the plane waves that can impinge upon the array. The noise and signal subspaces are orthogonal, so you can find which steering vectors correspond to the signals by looking for those steering vectors that are orthogonal to the noise subspace. So you have to reliably find the noise subspace, which is hard in practice (in my experience) because the dimensions of the signal and noise subspaces are not known in advance unless you know the number of impinging signals in advance.
  
  In cyclic MUSIC, the correlation matrix for the array data is replaced by either the non-conjugate cyclic correlation matrix or the conjugate cyclic correlation matrix. The structure of the matrices is similar–there is a signal subspace and a noise subspace. But the signal subspace dimension now reflects only the number of impinging signals that exhibit the chosen cycle frequency used in the cyclic correlation matrix estimate.
  
  I think the best expositions of cyclic MUSIC are by its originator Stephan Schell (full disclaimer: he’s a friend and former colleague). Check out The Literature [R104, R105, R106].