Fast TDOA and FDOA Estimation for Coherent Pulse Signals

Abstract

Time difference of arrival (TDOA) and frequency difference of arrival (FDOA) have been widely used for localizing temporally continuous signals passively. Temporal sparsity of pulse signals makes their TDOA and FDOA estimation processes much different, and computational complexity is a major concern in this area. This paper addresses the problem of fast TDOA and FDOA estimation of pulse signals and focuses mainly on narrowband coherent pulses. By decoupling the effects of TDOA and FDOA in the cost function of localization approximately, we propose a fast coarse TDOA and FDOA estimation method. The estimates are then refined with the cross-ambiguity function (CAF) algorithm within a small TDOA and FDOA neighborhood. In the simulations, the proposed method is demonstrated to have satisfying TDOA and FDOA estimation precisions, and it exceeds existing counterparts largely in computational efficiency.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In the field of signal processing, passive localization has been the subject of extensive research for decades [1, 2]. Time difference of arrival (TDOA) and frequency difference of arrival (FDOA) are two of the major intermediate parameters used for source localization [3–8]. By formulating a cost function of the observations with respect to all unknown parameters, the maximum likelihood method can be used to estimate the location-dependent parameters of TDOA and FDOA. However, as some redundant parameters are also included in the cost function (e.g., signal waveforms), and they usually have extremely high dimensions, the maximum likelihood algorithm can hardly be implemented via exhaustive searching. Therefore, the cross-ambiguity function (CAF) algorithm is often used as an alternate to save computational loads during TDOA and FDOA estimation [9–12], which can be combined with discrete Fourier transform (DFT) to further improve computational efficiencies [13].

Existing TDOA and FDOA estimation methods are mostly proposed for continuous signals [6–12]. Continuous signals typically have durations on the order of a few milliseconds to tens of milliseconds in the fields of communications and acoustics. During such a short period of time, only negligible position changes are introduced between emitters and receivers. Therefore, it is reasonable to treat the TDOA and FDOA as static parameters.

When a pulse-radiating emitter, like radar, is the one that needs to be localized, the duration of each pulse signal is often on the order of microseconds [14]. According to the CRLBs given in [15], such short signals may result in very low FDOA estimation precisions, which hardly meet the requirement of high-precision emitter localization. Therefore, multiple pulse signals of the same emitter should be utilized jointly to improve TDOA and FDOA estimation precisions. As the pulse repetition intervals (PRIs) of many radars are usually on the order of submilliseconds [14], if the number of pulses is not too many, the total time span of observed signals will be on the order of a few milliseconds to tens of milliseconds. Such a time span is comparative with that of continuous communication and acoustic signals. Vehicle-mounted, shipborne, or even airborne emitters and receivers move very slightly during such a short time span; thus the localization problem of pulse-radiating sources can also be addressed approximately in a static scenario. The CAF-based TDOA and FDOA estimation method for continuous signals can be extended to pulse signals with moderate modifications. However, continuous signals generally have modulation-dependent nonstationary waveforms, and no special temporal structures can be exploited to facilitate the TDOA and FDOA estimation process.

Pulse signals are temporally sparse; in the recent years, some research has focused on their TDOA and FDOA estimation. In [16], the authors introduce the Chirp Z- transform interpolation technique to improve estimation efficiency. In [17], the fact that FDOA is the time derivative of TDOA is exploited to eliminate FDOA ambiguity. In [18], a Keystone transform-based method is proposed to estimate TDOA and FDOA of broadband signals. More recently, Ma et al. made deep insights into the TDOA and FDOA estimation problem of interleaved pulse trains [19]. These works effectively use the temporal sparsity of pulse signals to facilitate the TDOA and FDOA estimation process. Besides temporal sparsity, pulse signals also have some other characteristics that are potentially helpful for parameter estimation )e.g., interpulse coherence). In many scenarios, radar emits coherent pulse signals in the same direction [20, 21], so that their echoes at the radar receiver will be phase coherent and can be accumulated to improve signal detection and parameter estimation performance. However, whether and how the interpulse coherence of pulse signals can be utilized in passive systems to improve TDOA and FDOA estimation performance is still an unanswered question. This paper takes a step further along this direction; we make use of the interpulse coherence to accelerate the process of TDOA and FDOA estimation for pulse-radiating emitters.

In this paper, we mainly address the TDOA and FDOA estimation problem of coherent pulse signals, based on the assumption that the observation time is short and the target-receiver geometry changes negligibly; thus the TDOA and FDOA remain constant during the observation time. We propose utilizing subspace decomposition to recover the signal waveforms on the two receivers based on the coherence between different pulses, which helps to suppress additive noises on data samplings. The subspace decomposition procedure also roughly decouples the TDOA and FDOA in the observation model, so that they can be estimated independently via two separated numerical procedures. If the coupling between the TDOA and the FDOA needs to be considered more rigorously, the coarse TDOA and FDOA estimates can be further refined via CAF searching in a small neighborhood. Different from the majority of existing CAF-based TDOA and FDOA estimation methods [8–12], the CAF of the proposed method is calculated based on the recovered signal waveforms on the two receivers, so it avoids repeated CAF calculations between multiple pulse pairs collected by the two receivers; thus its computational burden is significantly reduced.

The rest of the paper mainly consists of four parts. Section 2 formulates the problem of TDOA and FDOA estimation for coherent pulse trains. In Section 3, a fast TDOA and FDOA estimation method for coherent pulses is proposed. In Section 4, simulations are carried out to demonstrate the TDOA and FDOA estimation performance of the proposed method. Section 5 summarizes the whole paper.

2. Problem Formulation for TDOA and FDOA Estimation of Pulse Signals

Pulse signals of emitters like radar usually do not contain much complicated modulations [14]; they usually contain consistent modulations, for example, linear frequency modulated (LFM) signals with the same modulation parameters. Therefore, they can be modeled as coherent signals [20, 21]. Further assume that the signal bandwidth is very small when compared with the carrier frequency; then the signals can be approximated as narrowband ones, and the Doppler frequency shifts of the signals at different stations remain unchanged during the observation time.

The signals collected by the two receivers can be formulated as$$\begin{array}{c}\text{(1)}& x_{1,k}t=a_k{e}^{j{\phi }_k}st+u_{1,k}t,x_{2,k}t=b_k{e}^{j{\phi }_k}{e}^{j{\varphi }_0}{e}^{j2\pi f_{d}t}st-t_d+u_{2,k}t,\\ t\in T_k^1,T_k^2;k=1,\dots ,K,\end{array}$$where $x_{1,k}t$ and $x_{2,k}t$ are the observation data of the $k$ th pulse at time $t$ on the two receivers, $K$ is the number of pulses, and $st$ is the waveform of the first pulse on the first receiver. The frequency offset of the pulse signals at the first receiver is not explicitly indicated in (1); that is because only the frequency difference between the two receivers is considered. $a_k,b_k\in ℝ$ denote the relative amplitudes of the $k$ th pulse signal received by the two receivers with respect to that of the first pulse received by the first receiver, which satisfies $a_1=1$. $t_d$ and $f_d$ are the constant TDOA and FDOA of the pulses between the two receivers. The Gaussian white noises at the two receivers are denoted by $u_{1,k}t$ and $u_{2,k}t$, respectively, which satisfy $E{u_{1,k}t}^2=E{u_{2,k}t}^2={\sigma }^2$ with the noise variance ${\sigma }^2$ assumed to be known. When the noise variance is unknown, similar conclusions can be obtained as these for known ${\sigma }^2$ following the guidelines of analyses in [15]. $T_k^1$ and $T_k^2$ denote the starting and ending time instants of the $k$ th pulse, respectively. ${\varphi }_0$ stands for the initial phase offset between the two receivers caused by the asynchronization between different sampling clocks, and ${\phi }_k\in \mathrm{0,2}\pi$ represents the additional phase shift of the $k$ th pulse signal relative to the first one satisfying ${\phi }_1=0$.

The duty cycle of most pulse signals is very small and interpulse gaps should be skipped over to suppress noises and only observations within pulses should be retained for TDOA and FDOA estimation. Denote the sampling interval by $T$; then the discrete samples at the two receivers are given by$$\begin{array}{c}\text{(2)}& x_{1,k}n=a_k{e}^{j{\phi }_k}sn+u_{1,k}n,x_{2,k}n=b_k{e}^{j{\phi }_k}{e}^{j{\varphi }_0}{e}^{j{\varphi }_{k}v}{e}^{jvn}sn-\tau +u_{2,k}n,\\ n=1,\dots ,N_0;k=1,\dots ,K,\end{array}$$where $\tau =t_d/T$ and $v=2\pi f_{d}T$ are TDOA- and FDOA-dependent parameters, respectively. For notational simplicity, these two parameters are directly referred to as TDOA and FDOA in the rest of this paper unless otherwise stated. These two parameters will be estimated intermediately before obtaining the desired TDOA and FDOA estimates. ${\varphi }_{k}v=2\pi f_d{T}_k^1\triangleq v·{\overline{n}}_k$ represents the additional phase shift caused by the Doppler frequency shift at the beginning of the $k$ th pulse, where ${\overline{n}}_k=T_k^1/T$ and ${\overline{n}}_1=0$. In the formulations, we assume that the sampling interval is very small with respect to the pulse width, so that the nonalignment between the pulse starting edges with the nearest sampling instants does not significantly disturb the coherence of different pulses.

The samples of the $k$ th pulse at the two receivers can be written in a vector form as follows:$$\begin{array}{c}\text{(3)}& {\mathbf{x}}_{1,k}=a_k{e}^{j{\phi }_k}s+u_{1,k}\mathbf{,}{\mathbf{x}}_{2,k}=b_k{e}^{j{\phi }_k}{e}^{j{\varphi }_0}{e}^{j{\varphi }_{k}v}{\mathbf{D}}_{v}s-\tau +u_{2,k},\\ k=1,\dots ,K,\end{array}$$where ${\mathbf{x}}_{i,k}={x_{i,k}1,\dots ,x_{i,k}{N}_0}^T$ and $u_{i,k}={u_{i,k}1,\dots ,u_{i,k}{N}_0}^T$ for $i=\mathrm{1,2}$, $s={s1,\dots ,s{N}_0}^T$, $s-\tau ={s1-\tau ,\dots ,s{N}_0-\tau }^T$, ${\mathbf{D}}_v=\text{diag}\exp jvl$, with $\mathbf{l}={\mathrm{0,1},\dots ,N_0-1}^T$, and $N_0$ being the sample number within each pulse. In practical applications, the pulses of two receivers can be aligned first according to their arriving time, and only the residual TDOA requires precise estimation. Therefore, it is reasonable to assume that the TDOA is smaller than the sampling instant and much smaller than the pulse width; thus the edge effect that the signals at the two receivers are not accurately aligned due to time delay can be ignored without introducing significant model errors, and the delayed signal $s-\tau$ can be approximated by a spectrally group-delayed form of the original signal [15], which is given by$$\begin{array}{}\text{(4)}& s-\tau =F^H{D}_{\tau }\mathrm{Fs,}\end{array}$$where $\mathbf{F}=1/\sqrt{N_0}\exp -j2\pi /N_{0}l{l}^T$ and ${\mathbf{D}}_{\tau }=\text{diag}\exp -j2\pi \tau /N_0\mathbf{l}$. The superscripts $g^T$, $g^{\ast }$, and $g^H$ denote transpose, conjugate, and conjugate transpose operators, respectively. In (4), the time-delayed signal $s-\tau$ is interpreted as a variant of $s$ after a series of transformations, where $s$ is first transformed to the spectral domain by left-multiplying the Fourier matrix $\mathbf{F}$; then each spectral component is group-delayed by a TDOA $\tau$ and finally transformed back to the temporal domain by left-multiplying the conjugate transposed Fourier matrix ${\mathbf{F}}^H$. One can demonstrate (4) by left-multiplying both sides with $\mathbf{F}$; then it becomes more apparent that each spectral component of $s-\tau$ is a phase-shifted replica of the corresponding component of $s$, and the shifted phases are time-delay and frequency dependent and they form ${\mathbf{D}}_{\tau }$.

Unknown parameters contained in the above observation model include amplitudes $\mathbf{a}={a_2,\dots ,a_K}^T$ and $\mathbf{b}={b_1,b_2,\dots ,b_K}^T$ and phases $\phi ={{\phi }_2,\dots ,{\phi }_K}^T$ for different pulses, and $\eta ={s_r^T,s_i^T}^T$ and $\theta ={{\varphi }_0,\tau ,v}^T$ shared by all pulses, where subscripts $g_r$ and $g_i$ denote the real and imaginary parts of a variable, respectively. The whole parameter set in the observation model is$$\begin{array}{}\text{(5)}& \xi ={{\eta }^T,a^T,b^T,{\phi }^T,{\theta }^T}^T.\end{array}$$

Denote ${\mathbf{x}}_k={{\mathbf{x}}_{1,k}^T,{\mathbf{x}}_{2,k}^T}^T$, $\mathbf{x}={{\mathbf{x}}_1^T,\dots ,{\mathbf{x}}_K^T}^T$; then $\mathbf{x}$ has a complex Gaussian distribution in the case of white Gaussian noise, which is given by$$\begin{array}{}\text{(6)}& \mathbf{x\sim }ℂ\mathbb{N}\mathbf{\mu },{\sigma }^2{\mathbf{I}}_{2K{N}_0},\end{array}$$where$$\begin{array}{c}\text{(7)}& \mu ={{\mu }_1^T,\dots ,{\mu }_K^T}^T,{\mu }_k=\begin{array}{c}{a}_k{e}^{j{\phi }_k}s\\ b_k{e}^{j{\phi }_k}{e}^{j{\varphi }_0}{e}^{j{\varphi }_{k}v}{D}_v{F}^H{D}_{\tau }\mathrm{Fs}\end{array}\triangleq \begin{array}{c}{\mu }_{k,1}\\ {\mu }_{k,2}\end{array}.\end{array}$$

The task of TDOA and FDOA estimation is to estimate $\tau$ and $v$ based on the observations of the two sensors, that is, $x_{\mathrm{1,1}},x_{\mathrm{1,2}},\dots ,x_{1,K}$ and $x_{\mathrm{2,1}},x_{\mathrm{2,2}},\dots ,x_{2,K}$.

3. Fast TDOA and FDOA Estimation for Coherent Pulse Trains

Based on the distribution function in (6), the maximum likelihood algorithm can be applied to obtain optimal TDOA and FDOA estimates via exhaustive searching in the unknown parameter space. However, as the dimension of $\xi$, that is, $2N_0+3K+1$, is usually very high, the computational complexity of this exhaustive searching procedure will be prohibitive in most practical applications. The CAF technique has been introduced to significantly improve the computational efficiency of the maximum likelihood method [9–12]; it only takes $\tau$ and $v$ as unknown parameters. Unfortunately, after the simplification, the computational load of the TDOA and FDOA estimation process is still very heavy when the required precisions of the estimates are high [13]. In this section, we propose a fast method to realize TDOA and FDOA estimation for coherent pulse trains. The method consists of a coarse parameter estimation procedure and a refinement one. The coarse-to-fine strategy has previously been used for radar target detection [22], and it can also be combined with DFT techniques [14] to further improve its computational efficiency.

3.1. Coarse TDOA and FDOA Estimation

In order to estimate the TDOA and FDOA of coherent pulse trains, the observation data at the two receivers within the $K$ pulses are rewritten in a matrix form as follows:$$\begin{array}{c}\text{(8)}& {\mathbf{X}}_1={\mathbf{x}}_{\mathrm{1,1}},\dots ,{\mathbf{x}}_{1,K}=s\times a_1{e}^{j{\phi }_1},\dots ,a_K{e}^{j{\phi }_K}+u_{\mathrm{1,1}},\dots ,u_{1,K}\triangleq s{\alpha }_1^H+{\mathbf{U}}_1,{\mathbf{X}}_2={\mathbf{x}}_{\mathrm{2,1}},\dots ,{\mathbf{x}}_{2,K}={\mathbf{D}}_v{\mathbf{F}}^H{\mathbf{D}}_{\tau }\mathbf{Fs}\times e^{j{\varphi }_0}{b}_1{e}^{j{\phi }_1}{e}^{j{\varphi }_{1}v},\dots ,b_K{e}^{j{\phi }_K}{e}^{j{\varphi }_{K}v}\\ +u_{\mathrm{2,1}},\dots ,u_{2,K}\triangleq s{}^{\prime }{\alpha }_2^H+{\mathbf{U}}_2,\end{array}$$where ${\alpha }_1={a_1{e}^{j{\phi }_1},\dots ,a_K{e}^{j{\phi }_K}}^H$, ${\mathbf{U}}_1=u_{\mathrm{1,1}},\dots ,u_{1,K}$, $s{}^{\prime }={\mathbf{D}}_v{\mathbf{F}}^H{\mathbf{D}}_{\tau }\mathbf{Fs}$, ${\alpha }_2=e^{-j{\varphi }_0}{b_1{e}^{j{\phi }_1}{e}^{j{\varphi }_{1}v},\dots ,b_K{e}^{j{\phi }_K}{e}^{j{\varphi }_{K}v}}^H$, ${\mathbf{U}}_2=u_{\mathrm{2,1}},\dots ,u_{2,K}$.

The likelihood function of the observed data matrix is$$\begin{array}{}\text{(9)}& p{\mathbf{X}}_1,{\mathbf{X}}_2|\mathbf{\xi }={\pi {\sigma }^2}^{-2N_{0}K}\exp -\frac{1}{{\sigma }^2}{{\mathbf{X}}_1-s{\alpha }_1^H}_2^2+{{\mathbf{X}}_2-s{}^{\prime }{\alpha }_2^H}_2^2.\end{array}$$

The maximum of the likelihood function is equivalent to the minimum of the following cost function when ${\sigma }^2$ is known:$$\begin{array}{}\text{(10)}& W\mathbf{\xi }={{\mathbf{X}}_1-s{\alpha }_1^H}_2^2+{{\mathbf{X}}_2-s{}^{\prime }{\alpha }_2^H}_2^2.\end{array}$$

If ${\sigma }^2$ is unknown, its maximum likelihood estimate can be calculated and substituted into (9), and the same cost function as the one in (10) can also be obtained. Moreover, if the noise variances at the two receivers are different, one can introduce their diverse values to reformulate (9) and then obtain a cost function accordingly similar to that in (10), and then continue with the coarse parameter estimation procedure.

In a noiseless case, the minima of $W_1\xi ={{\mathbf{X}}_1-s{\alpha }_1^H}_2^2$ and $W_2\xi ={{\mathbf{X}}_2-s^{\prime }{\alpha }_2^H}_2^2$ have consistent locations corresponding to the true values of the unknown parameters. In practical applications, as the observations are noise contaminated, each of the two separated subfunctions has a Gaussian distribution centralized at the minimum. We propose to treat ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$ as independent observations at first, and recover $s$, $s^{\prime }$, ${\alpha }_1$, and ${\alpha }_2$ coarsely based on them, and then exploit the $\tau$- and $v$-dependent relations between them to estimate TDOA and FDOA.

After ignoring the relations between ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$, it becomes infeasible to obtain unique estimates of $s$, $s^{\prime }$, ${\alpha }_1$, and ${\alpha }_2$; that is because $s{\alpha }_1^H=s{e}^{j{\phi }_0}{{\alpha }_1{e}^{j{\phi }_0}}^H$ and $s^{\prime }{\alpha }_2^H=s{}^{\prime }{e}^{j{\phi }_0}{{\alpha }_2{e}^{j{\phi }_0}}^H$ for any ${\phi }_0$. Therefore, we substitute them with ${\mathbf{r}}_1$, ${\mathbf{r}}_2$, ${\mathbf{w}}_1$, and ${\mathbf{w}}_2$ to distinguish the intermediate estimates from the original $\tau$- and $v$-dependent vectors.

In the two cost subfunctions, the minuends of $s{\alpha }_1^H$ and $s{}^{\prime }{\alpha }_2^H$ are both rank-1, so maximum likelihood estimates of ${\mathbf{r}}_1$, ${\mathbf{r}}_2$, ${\mathbf{w}}_1$, and ${\mathbf{w}}_2$ can be obtained according to the principal components of ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$; that is,$$\begin{array}{}\text{(11)}& X_1\stackrel{{\lambda }_{\max }}{⟶}{\widehat{r}}_1,{\widehat{w}}_1,\text{(12)}& {\mathbf{X}}_2\stackrel{{\lambda }_{\max }}{⟶}{\widehat{r}}_2,{\widehat{w}}_2.\end{array}$$

In the above estimation formula, ${\lambda }_{\max }$ represents the largest eigenvalue of the corresponding matrix, and $\stackrel{{\lambda }_{\max }}{⟶}$ indicates that eigendecomposition is performed on the matrix to obtain the left and right eigenvectors associated with the largest eigenvalue. The eigendecomposition process can be realized with the svd function in MATLAB or via a two-step numerical procedure, which first decomposes the covariance matrices ${\mathbf{X}}_1^H{\mathbf{X}}_1$ and ${\mathbf{X}}_2^H{\mathbf{X}}_2$ to estimate ${\widehat{w}}_1$ and ${\widehat{w}}_2$, and then estimate ${\widehat{r}}_1$ and ${\widehat{r}}_2$ via ${\widehat{r}}_1={\mathbf{X}}_1{\widehat{w}}_1$ and ${\widehat{r}}_2={\mathbf{X}}_2{\widehat{w}}_2$. The eigendecomposition process is scale and phase invariant; that is, ${\widehat{r}}_1{\widehat{w}}_1^H=\eta e^{j{\phi }_0}{\widehat{r}}_1{{\eta }^{-1}{e}^{j{\phi }_0}{\widehat{w}}_1}^H$ and ${\widehat{r}}_2{\widehat{w}}_2^H=\eta e^{j{\phi }_0}{\widehat{r}}_2{{\eta }^{-1}{e}^{j{\phi }_0}{\widehat{w}}_2}^H$ for any $\eta$ and ${\phi }_0$. In order to guarantee the uniqueness of the estimates, we include some additional constraints for the eigenvectors; that is, ${\widehat{w}}_{1}1={\widehat{w}}_{2}1=1$ and ${{\widehat{r}}_1}_2={{\widehat{r}}_2}_2$, with ${\widehat{w}}_{1}1$ and ${\widehat{w}}_{2}1$ representing the first element of ${\widehat{w}}_1$ and ${\widehat{w}}_2$, respectively.

During the eigendecomposition process of ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$, the element-wise phase and amplitude relations between $s$ & $s{}^{\prime }$ and ${\alpha }_1$ & ${\alpha }_2$ may have not been accurately retained in the estimates of ${\mathbf{r}}_1$ & ${\mathbf{r}}_2$ and ${\mathbf{w}}_1$ & ${\mathbf{w}}_2$, but the phase relations along the vectors are retained, which are $\tau$- and $v$-dependent and can be exploited for TDOA and FDOA estimation. According to the phase relations between $s$ and $s{}^{\prime }$, and that between ${\alpha }_1$ and ${\alpha }_2$, the FDOA and TDOA estimates can be obtained based on ${\widehat{r}}_1$, ${\widehat{r}}_2$, ${\widehat{w}}_1$, and ${\widehat{w}}_2$ as follows:$$\begin{array}{}\text{(13)}& \widehat{v}=\arg \underset{v}{\min }{\mathbf{h}v-\exp j\times \mathrm{angle}{\widehat{w}}_1\odot {\widehat{w}}_2^{\ast }}_2^2,\text{(14)}& \widehat{\tau }=\arg \underset{\tau ,{\phi }_0}{\min }{e^{j{\phi }_0}\mathbf{q}\tau -\exp j\times \mathrm{angle}\mathbf{F}{\mathbf{D}}_{\widehat{v}}^{\ast }{\widehat{\mathbf{r}}}_2\odot {\mathbf{F}{\widehat{\mathbf{r}}}_1}^{\ast }}_2^2,\end{array}$$where $\mathbf{h}v={e^{j{\varphi }_{1}v},\dots ,e^{j{\varphi }_{K}v}}^T$, the definition of ${\varphi }_{k}v$ is the same as that in (2), $\mathbf{q}\tau =\exp -j2\pi \tau /N_0\mathbf{l}$, and $\odot$represents element-wise multiplication between vectors; $\mathrm{angle}g$ indicates the argument of a complex number. A phase shift is added in both (13) and (14) to compensate for the phase mismatching between the two items in each cost function. A phase of $e^{-j{\varphi }_{1}v}$ should be added in (14) to shift the phase of the first element of $\mathbf{h}v$ to 0 because ${\widehat{w}}_{1}1={\widehat{w}}_{2}1=1$, but as ${\varphi }_{1}v=0$, according to (2), this item is left out from (13). The initial phase of the subtrahend in (14) is not predictable; thus an unknown phase of $e^{j{\phi }_0}$ is added in (14), and the phase is optimized jointly with $\tau$. The estimation equations in (13) and (14) decouple originally interconnected TDOA and FDOA variables, and they can be realized fast via DFT. Therefore, coarse TDOA and FDOA estimates can be obtained from (13) and (14) via computationally very efficient algorithms.

The coarse TDOA/FDOA estimation algorithm and the refinement procedure afterwards both work based on the assumption of coherent pulse signals. If the signals are only partially coherent or even incoherent, the result of the eigendecomposition process in (11) and (12) will be unpredictable, one will not be able to extract TDOA/FDOA-dependent measurements easily by skipping over the large amount of waveform-dependent redundant variables, and the coarse TDOA/FDOA estimation procedure fails accordingly. Fortunately, the result of the eigendecomposition process contains clues about the coherence of pulse signals. Only one of the eigenvalues is significantly larger than zero in the case of completely coherent signals, while more than one eigenvalue is large if the assumption of coherence is deviated. In the latter case, one may have to retrogress to the original CAF-based method for TDOA/FDOA estimation, or single out coherent pulse signals via an eigen-value-checking criterion to implement coarse TDOA/FDOA estimation first, and then take all pulses into consideration to refine the estimates.

3.2. TDOA and FDOA Refinement

As the relations between ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$ are artificially ignored when recovering $s$, $s{}^{\prime }$, ${\alpha }_1$, and ${\alpha }_2$, the estimates of $\tau$ and $v$ are probably suboptimal and we call them coarse estimates. The estimates can be refined by taking the $\tau$– and $v$–dependent relations between ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$, together with those between ${\mathbf{w}}_1$ and ${\mathbf{w}}_2$, into consideration. In order to improve parameter estimation accuracy, an additional two-dimensional searching procedure in the TDOA-FDOA domain is performed in a small neighborhood of the coarse TDOA and FDOA estimates, and an interpolation procedure is implemented around the peak of the resulting CAF to refine the parameter estimates.

In order to save computational load, the recovered signal waveforms of ${\widehat{r}}_1$ and ${\widehat{r}}_2$ are used to replace the original pulse observations when calculating the CAF; that is,$$\begin{array}{}\text{(15)}& \text{CAF}\tau ,v={\widehat{r}}_1^H{\mathbf{F}}^H{\mathbf{D}}_{\tau }^H\mathbf{F}{\mathbf{D}}_v^H{\widehat{r}}_2\times h^{H}v\times \exp j\times \mathrm{angle}{\widehat{w}}_1\odot {\widehat{w}}_2^{\ast }.\end{array}$$

The cost function in (14) exploits jointly the observations at the two receivers, which are used for independent parameter estimation in (12) and (13), and takes into account the coupling between TDOA and FDOA, so it is expected to obtain higher parameter estimation accuracies. From the perspective of computational efficiency, the proposed method needs to compute the CAF only within a small neighborhood of the coarse estimates obtained by (12) and (13), and the calculations at each CAF grid contain only a correlation between ${\widehat{r}}_1$ and a time-delayed and frequency-shifted replica of ${\widehat{r}}_2$. On the contrary, original CAF-based TDOA and FDOA estimation methods have no coarse estimates, and they compute CAF between all the $K$ pulse signal pairs at the two receivers. Therefore, the computational load of the proposed method is expected to be largely cheaper than existing CAF-based methods. The cross-correlation between ${\widehat{r}}_1$ and ${\widehat{r}}_2$ can be speeded by combining the CAF calculations with the DFT technique to further improve the overall computational efficiency [13].

Another major difference between the proposed subspace-based estimation method and the original CAF-based one is that, when estimating FDOA using (12), unambiguous estimates can be obtained only within a scope depending on the pulse repetition intervals (PRIs). For example, when the pulse train has a constant interval of 1 ms, the equation can only estimate unambiguous FDOA in a range of 1 KHz. If the candidate FDOA exceeds this range, all ambiguous FDOA estimates in (12) should be retained for further identification. The local FDOA estimates are then substituted into (14) to obtain multiple pairs of coarse TDOA-FDOA estimates. The ambiguity is finally eliminated to yield a global optimal estimate according to (15). This process increases the computational load of the proposed method by a multiple equaling the number of unambiguous FDOA estimates. However, in despite of the existence of local minima, as the pulse number is large in most practical applications, the overall computational complexity of the proposed method is still significantly lower than that of the original CAF-based method. In practice, the value of the FDOA is constrained by various factors, such as relative target-receiver speed and observation geometry, making its range much limited; thus the number of unambiguous FDOA estimates and TDOA-FDOA pairs will be very small. In addition, the ambiguity effect will be further reduced significantly when the pulse intervals are not constant.

3.3. Computational Complexity Analysis

The proposed method mainly includes three steps. First, eigendecomposing ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$ to estimate ${\mathbf{r}}_1$, ${\mathbf{r}}_2$, ${\mathbf{w}}_1$, and ${\mathbf{w}}_2$ according to (11) and (12), which requires $O{K}^3+2K{N}_0$ complex multiplications. Second, estimating $\tau$ and $v$ coarsely according to (13) and (14), which requires only some numerical calculations and the computational complexity is negligible. Third, refining the estimates of $\tau$ and $v$ by calculating the CAF according to (15), which requires $2N_0^2$ complex multiplications at each TDOA-FDOA grid and can be accelerated via FFT. As the number of samplings in each pulse is generally much larger than that of pulses, that is, $N_0\gg K$, the computational complexity of the proposed method majorly lies in the CAF calculation procedure.

4. Simulations and Analyses

In this part, we carry out simulations to demonstrate the performance of the proposed TDOA and FDOA estimation method with respect to various factors, such as pulse number, SNR, pulse width, and PRI. Theoretical analyses on how these factors affect the performance are provided in another paper by the same authors [23], and the results are included in this part for theoretical verification. Assume that the pulse signals are linear frequency modulated (LFM) with a bandwidth of 1 MHz, and the signal carrier frequency is 1 GHz. Therefore, the signals are approximately narrowband. The received signals are downconverted to a low intermediate frequency and then sampled with a frequency of 10 MHz. The time delay of the signals at the two receivers is 0, and the frequency shift is 1 kHz.

In the simulations, the original CAF-based TDOA and FDOA estimation method [10] is carried out for performance comparison. The TDOA and FDOA searching steps for calculating the CAF are set to be equal to their corresponding CRLBs [23], and the search ranges are centered at the true TDOA and FDOA values and extended by 10 steps on both sides. The CRLBs used for setting the searching grids are calculated according to the signal waveforms for convenience, and they can be substituted with their estimates obtained in the proposed method in practical applications. Numerical interpolation [24] is then implemented near the peak of the CAF to obtain final TDOA and FDOA estimates. Two implementations of the proposed method are included for performance evaluation, named Subspace method and Subspace-CAF method. The TDOA and FDOA estimates of the Subspace method are directly obtained from (12) and (13), while those of the Subspace-CAF method are obtained from (14). In the Subspace-CAF method, as coarse TDOA and FDOA estimates have been obtained in the first stage, the search ranges of the subsequent CAF-based refinement procedure are reduced to the CRLB on both sides of the coarse estimates, and the search step is reduced to 1/10 times of the CRLB. This parameter setting ensures that the number of CAF search grids in the proposed method is equal to that of the original CAF-based method. The number of simulations in each scenario is 10 000.

In the first group of simulations, we fix the pulse width at 30us, the SNR at 5 dB, and the PRI at 100us and increase the pulse number from 5 to 80. The TDOA and FDOA estimation RMSE of the three methods are obtained and shown in Figures 1(a) and 1(b). The results show that the TDOA and FDOA estimates obtained by the Subspace method are close to their true values; they have high precisions but still deviate from the theoretical lower bounds with significant margins. By refining the TDOA and FDOA estimates via local CAF searching, their precisions can be further improved to approach the CRLB well. The precision improvement of Subspace-CAF over Subspace is rooted in the joint exploitation of the observations at the two receivers; the coupling between TDOA and FDOA is taken into account in this procedure and higher parameter estimation accuracies are thus obtained. When compared with the original CAF-based method, the Subspace-CAF method has a slightly deteriorated TDOA estimation precision when the number of pulses is small, but its FDOA estimation accuracy is higher than that of the CAF-based method.

[figure(s) omitted; refer to PDF]

Figure 1(c) shows the average time of the three methods in implementing a single TDOA and FDOA estimation simulation. It can be seen that the computational loads of the two proposed methods remain stable when pulse number increases, and the CAF-based refinement procedure aggravates the computational complexity of the method by about 4 times. On the contrary, the computational complexity of the traditional CAF-based method increases superlinearly with the number of pulses. When the pulse number increases to 80, its computational load is about 7 times that of the Subspace-CAF method, which has comparable parameter estimation precisions. This is mainly because the CAF computing process of the Subspace-CAF method only calculates the cross-correlation between the two recovered signal waveforms, while that of the original CAF-based method calculates the cross-correlation between all pulse pairs at the two receivers. This result verifies the significant advantages of the proposed method in computational efficiency.

Based on the above simulations, we then fix the number of coherent pulses at 40, and vary the SNR of the pulse signals from −5 dB to 35 dB. The TDOA and FDOA estimation RMSE of the three methods, together with the corresponding CRLB, are obtained and shown in Figure 2. Similar to the results in Figures 1(a) and 1(b), the Subspace method is able to obtain high-precision TDOA and FDOA estimates, but there is a certain gap between their precisions and the CRLB. The parameter estimation precisions of the Subspace-CAF method and the traditional CAF-based method approach the CRLB better. The Subspace-CAF method has a slightly lower TDOA estimation accuracy than that of the CAF method when the SNR is lower than 5 dB, and its FDOA estimation accuracy is slightly higher. The computational efficiencies of the three methods can be deduced from the results in Figure 1(c) by fixing the pulse number at 40, and they do not vary significantly with SNR, which again demonstrates the advantage of the proposed methods in computational efficiency. In this and the following simulations, the results on computational efficiencies of different methods can be inferred from Figure 1(c), so they are excluded to avoid redundancy.

[figure(s) omitted; refer to PDF]

Then we fix the number of coherent pulses at 40 and the SNR on both receivers at 5 dB and then vary the pulse width from 5 us to 30 us. The TDOA and FDOA estimation accuracies of the three methods are shown in Figures 3(a) and 3(b). The comparisons of the TDOA and FDOA estimation accuracies of the three methods are similar to those in Figures 1 and 2.

[figure(s) omitted; refer to PDF]

Finally, we fix the number of coherent pulses at 40, the SNR on both receivers at 5 dB, and the pulse width at 30us and increase the PRI from 0.1 ms to 10 ms. The TDOA and FDOA estimation RMSE of the three methods are shown in Figures 4(a) and 4(b). The comparisons of TDOA and FDOA estimation accuracies of the three methods are similar to the results in Figures 1–3.

[figure(s) omitted; refer to PDF]

5. Conclusions

In this paper, a fast TDOA and FDOA estimation method is proposed for pulse signals. It decouples TDOA and FDOA in the cost function approximately to speed up the parameter estimation process and then put the coarse estimates into the original cost function to realize refined TDOA and FDOA estimation. Simulation results show that the proposed method exceeds the traditional CAF-based method largely in computational efficiency, and the refined TDOA and FDOA estimates have satisfying precisions that are comparable with those of the traditional CAF-based method. Specifically, when the overall SNR of the observed data is low in cases of small pulse numbers and narrow pulse widths, the TDOA estimation accuracy of the proposed method is slightly inferior to that of the traditional CAF-based method, while its FDOA estimation accuracy is slightly superior to that of the latter method.

References

[1] F. Ma, Z. M. Guo, F. Guo, "Direct position determination for wideband sources using fast approximation," IEEE Transactions on Vehicular Technology, vol. 68 no. 8, pp. 8216-8221, DOI: 10.1109/tvt.2019.2921981, 2019.

[2] F. Ma, Z. M. Guo, F. Guo, "Direct position determination in asynchronous sensor networks," IEEE Transactions on Vehicular Technology, vol. 68 no. 9, pp. 8790-8803, DOI: 10.1109/tvt.2019.2928638, 2019.

[3] K. Deergha Rao, D. C. Reddy, "A new method for finding electromagnetic emitter location," IEEE Transactions on Aerospace and Electronic Systems, vol. 30 no. 4, pp. 1081-1085, DOI: 10.1109/7.328756, 1994.

[4] K. Becker, "Passive localization of frequency-agile radars from angle and frequency measurements," IEEE Transactions on Aerospace and Electronic Systems, vol. 35 no. 4, pp. 1129-1144, DOI: 10.1109/7.805432, 1999.

[5] K. Becker, "An efficient method of passive emitter location," IEEE Transactions on Aerospace and Electronic Systems, vol. 28 no. 4, pp. 1091-1104, DOI: 10.1109/7.165371, 1992.

[6] K. C. Ho, Y. Chan, "Geolocation of a known altitude object from TDOA and FDOA measurements," IEEE Transactions on Aerospace and Electronic Systems, vol. 33 no. 3, pp. 770-783, DOI: 10.1109/7.599239, 1997.

[7] K. C. Ho, W. Xu, "An accurate algebraic solution for moving source location using TDOA and FDOA measurements," IEEE Transactions on Signal Processing, vol. 52 no. 9, pp. 2453-2463, DOI: 10.1109/tsp.2004.831921, 2004.

[8] K. C. Ho, X. Lu, L. Kovavisaruch, "Source localization using TDOA and FDOA measurements in the presence of receiver location errors: analysis and solution," IEEE Transactions on Signal Processing, vol. 55 no. 2, pp. 684-696, DOI: 10.1109/tsp.2006.885744, 2007.

[9] S. Liu, T. Shan, R. Tao, Y. D. Zhang, G. Zhang, F. Zhang, Y. Wang, "Sparse discrete fractional Fourier transform and its applications," IEEE Transactions on Signal Processing, vol. 62 no. 24, pp. 6582-6595, DOI: 10.1109/tsp.2014.2366719, 2014.

[10] S. Stein, "Algorithms for ambiguity function processing," IEEE Transactions on Acoustics, Speech, & Signal Processing, vol. 29 no. 3, pp. 588-599, DOI: 10.1109/tassp.1981.1163621, 1981.

[11] R. J. Ulman, E. Geraniotis, "Wideband TDOA/FDOA processing using summation of short-time CAF’s," IEEE Transactions on Signal Processing, vol. 47 no. 12, pp. 3193-3200, DOI: 10.1109/78.806065, 1999.

[12] A. Ramachandra, "Cross Ambiguity Function for Emitter Location, Master Thesis," 2008.

[13] K. P. Bentz, "Computation of the Cross Ambiguity Function Using Perfect Reconstruction DFT Filter banks Doctor’s Dissertation," 2007.

[14] M. I. Skolnik, "Introduction to Radar Systems," 2001.

[15] E. Angel, E. Angel, "Joint TDOA and FDOA estimation: a conditional bound and its use for optimally weighted localization," IEEE Transactions on Signal Processing, vol. 59 no. 4, pp. 1612-1623, DOI: 10.1109/tsp.2010.2103069, 2011.

[16] T. Shan, S. Liu, Y. D. Zhang, M. G. Amin, R. Tao, Y. Feng, "Efficient architecture and hardware implementation of coherent integration processor for digital video broadcast‐based passive bistatic radar," IET Radar, Sonar & Navigation, vol. 10 no. 1, pp. 97-106, DOI: 10.1049/iet-rsn.2015.0006, 2016.

[17] S. Yao, Q. He, X. Ouyang, C. Xia, "A novel method for unambiguous DFO estimation of radar signals in LEO dual-satellite geolocation system," Journal of Astronautics (in Chinese), vol. 39 no. 11, pp. 1275-1283, 2018.

[18] X. Xiao, G. Fucheng, D. Feng, "Low‐complexity methods for joint delay and Doppler estimation of unknown wideband signals," IET Radar, Sonar & Navigation, vol. 12 no. 4, pp. 398-406, DOI: 10.1049/iet-rsn.2017.0368, 2018.

[19] F. Ma, Z. M. Liu, F. Guo, D. Yang, L. Yang, "Joint TDOA and FDOA estimation for interleaved pulse trains from multiple pulse radiation sources," IEEE Transactions on Aerospace and Electronic Systems, vol. 56 no. 5, pp. 4099-4111, DOI: 10.1109/taes.2020.2987050, 2020.

[20] M. Pourhomayoun, M. L. Fowler, "Cramer-Rao lower bound for frequency estimation for coherent pulse train with unknown pulse," IEEE Transactions on Aerospace and Electronic Systems, vol. 50 no. 2, pp. 1304-1312, DOI: 10.1109/taes.2014.130024, 2014.

[21] G. B. Jordan, "Comparison of two major classes of coherent pulsed radar systems," IEEE Transactions on Aerospace and Electronic Systems, vol. AES-11 no. 3, pp. 363-371, DOI: 10.1109/taes.1975.308087, 1975.

[22] S. Liu, H. Zhang, T. Shan, Y. Huang, "Efficient radar detection of weak manoeuvring targets using a coarse‐to‐fine strategy," IET Radar, Sonar & Navigation, vol. 15 no. 2, pp. 181-193, DOI: 10.1049/rsn2.12028, 2021.

[23] L. Yan, "Performance analysis of TDOA and FDOA estimation for pulse signals," International Journal of Antennas and Propagation, vol. 2022,DOI: 10.1155/2022/7672417, 2022.

[24] R. Tao, E. Q. Chen, W. Q. Zhang, "Two-stage method for joint time delay and Doppler shift estimation," IET Radar, Sonar & Navigation, vol. 2 no. 1, pp. 71-77, DOI: 10.1049/iet-rsn:20060014, 2008.