# A Fast Direct Position Determination Algorithm for LFM Signal Based on Spectrum Detection

Authors: Han Ren; Huijie Liu; Rujiang Guo
Journal: Journal of Sensors
 (2022)

---

## Abstract

For the purpose of simultaneously estimating and locating the linear frequency modulation (LFM) emitter with unknown parameters, an innovative approach, i.e., Fast Direct Position Determination (FDPD), is introduced. The proposed algorithm is based on maximum likelihood estimation (MLE) and spectrum detection. To improve the accuracy and overcome the dramatic complexity of plain maximum likelihood formulation, we further derive the objective function equation of Direct Position Determination (DPD) algorithm and present an enhanced strategy to solve the highly nonlinear optimization problem. By combining the two-step localization method, one-step localization method, and short-time Fourier transform (STFT), our approach realizes jointly estimation of the transmitted signal parameters and emitter localization. Simulation results show that the proposed method is superior compared to the existing DPD, and two-step localization algorithms in terms of localization error and computational complexity, especially for low signal-to-noise ratio (SNR).


---

## Body

## 1. Introduction

Localization of emitters is a fundamental demand in many application fields, such as radar [1, 2], sonar [3, 4], global navigation satellite systems (GNSS) [5, 6], wireless sensor network [7, 8], and seismic exploration [9, 10]. Linear frequency modulation (LFM) signals, known as large time-bandwidth product and high resolution signals, are frequently used in localization systems and have attracted increasing attention of many scholars. As stealth technology and low intercept technology get prosperous and mature, passive emitter localization technology has gradually developed from the traditional single station to using several widely distributed receiver stations. This paper considers a passive localization problem for an emitter by multiple diffusely separated receiver stations.In general, the localization methods contain two main categories of two-step localization method and one-step localization method. Algorithms of two-step localization method [11–20] are commonly composed of two steps: first, each receiver station measures and computes the parameters related to the position of the target. The parameters involved in position include time difference of arrival (TDOA) [11–13], angle of arrival (AOA) [14–17], frequency difference of arrival (FDOA) [18–20], etc. The second step is to construct equations using the parameters to estimate the position of the target. Both steps of two-step method might be optimal; given their input, two-step method is still suboptimal. Since the correlation between different receiver stations is ignored, the precision of the parameter measurement in the first step is limited. Moreover, the parameters may be virtually absent in some of the receiver stations due to the signal strength fluctuation and different path attenuation, leading to data association and performance deterioration issues, especially in low signal-to-noise ratio (SNR) scenarios.One-step localization method [21–35] processes the observed signals in the central processor without parameter extraction and then constructs the cost function only related to the position of the emitter. This kind of localization algorithms finally obtain the position of the emitter through the maximum likelihood estimation (MLE), the least squares, the grid search, gradient-based methods, etc. Therefore, it is also referred to direct position determination (DPD) [21–24]. In 2004, Weiss first proposed direct position determination technology [21], analyzed, and compared the DPD performance with AOA for signals with known (DPD-known) and unknown waveforms (DPD-unknown) [22, 23]. Over the past dozen years, multiple DPD algorithms have been presented to enhance performance of passive localization system based on different approaches including the maximum likelihood (ML) [25, 26], the time frequency analysis (DPD-STFT-w) [27], the multiple signal classification (MUSIC) [28, 29], the expectation maximization (EM-DPD) [30, 31], the minimum variance distortionless response (MVDR) [32–34], and the time-varying delay [35]. Compared with the one-step approaches, the DPD algorithms directly estimate the initial position of the target according to the received signals in the central station, avoiding the process of correlation of parameters of different receivers. Moreover, DPD can locate the signals (simultaneous signals with same frequency, etc.) which are difficult to be processed by traditional methods. At the same time, because the cost function of DPD is only related to the position of the emitter, it has higher localization accuracy and better robustness at lower SNR on account of making full use of the prior information of the signals which are coming from the same emitter [23]. However, since the DPD algorithms deal with the original sampled signals rather than the analytic solutions of the emitter, the computational complexity increases rapidly and the estimation performance of waveform parameters deteriorates drastically when the SNR is lower than a threshold point, i.e., the “threshold effect.”This paper addresses passive localization of an LFM emitter with unknown signal parameters. Aiming at the problem of “threshold effect” and high computational complexity of DPD algorithms at low SNR, we proposed a Fast Direct Position Determination (FDPD) algorithm. The proposed algorithm can realize the combination of emitter localization and parameter estimation of the LFM signal. The main procedures of this approach include the following:(i)
Firstly, the signal spectrum is obtained through short-time Fourier transform (STFT) and other operations (linear fitting, the least square method, etc.); then, the SNR is calculated(ii)
Compared the SNR with the threshold point, the cooperative localization of two-step method and one-step method is combined to finally achieve the optimal localization. The performance of the two-step method and one-step method is similar at high SNR, we choose the two-step method owing to its less computational complexity. Chan algorithm [36] and Taylor algorithm [37], which are two representative algorithms of two-step localization method, are selected to complete collaborative localization: Chan algorithm is used to estimate the initial position coordinates, regarding it as the initial value, and Taylor series expansion algorithm is used to obtain the further position coordinates(iii)
At low SNR, the initial frequency, chirp rate, and other signal parameters are estimated by time frequency analysis, thus we can capture the entire parameter set of the LFM profile of the signal and obtain an objective function. Then, a four dimensional search (4D) is carried out at the position coordinates obtained by Chan and Taylor algorithms. To solve the highly nonlinear optimization problem (four-dimensional search), we turn the four-dimensional search into a two-dimensional parameter estimation combined with grid search. The search range is determined by the error analysis of Chan and Taylor algorithms. Simulation results show that the proposed algorithm reduces the computational complexity, and the positioning accuracy is better than Chan, Taylor, DPD, and other traditional algorithmsThe rest of this paper is organized as follows. Section2 introduces the model of the LFM signal for widely separated receiver stations. In Section 3, the localization problem formulated in an ML estimation paradigm is clarified. Section 4 presents the proposed FDPD algorithm in detail. Section 5 provides the numerical simulations. The conclusions are drawn in Section 6.


## 2. Signal Model

Consider a 2D Cartesian coordinate system where a noncooperative stationary emitter transmitting signals. We let the vectorsplxl,yl,l=1,2,⋯,L represent the positions of L widely separated receiver stations. The transmitter position is denoted by the vector of coordinates px,y. As shown in Figure 1, the emitter transmitted an LFM signal st. Then, the signal is intercepted by the receiver stations and transmitted to the center processor.Figure 1
Sketch map of an emitter and L widely distributed receiver stations.The LFM signal is modeled as(1)st;f0,k,Tp=expj2πf0t+12kt2,0<t<Tp,where f0 denotes the initial frequency of the LFM signal, k is the chirp rate, and Tp is the time duration. The signal observed by the lth receiver station is given by
(2)rlt=blαlpst−τlp−t0;f0,k,Tp+nlt.rlt is a time-dependent vector, bl is an unknown complex scalar representing the attenuation coefficient related with regard to the path from the emitter to the receiver stations, and αlp is the lth array response to signal transmitted from position p. st−τlp−t0;f0,k,Tp is the signal waveform, transmitted at time t0 and delayed by τlp. The vector nlt represents the independent and identically distributed zero-mean complex Gaussian white noise with corvariance matrix Rl=σl2I, and each receiving channel is independent of each other. τlp can be calculated by [35]
(3)τlp=x−xl2+y−yl2c,withc denoting the speed of light.After sampling, the continuous signal of Equation (2) can be expressed in a vector form
(4)rlj=blαlpslj+nlj,0≤j≤Ns−1,where
(5)slj≜st−τlp−t0t=jT,rlj=rltt=jT,nlj=nltt=jT.withNs denotes the number of the samples, fs denotes the sampling frequency, B is the signal bandwidth, and T is the time duration of the received signal, thus Ns=T∗fs (fs≥2B).


## 3. Problem Formulation about Direct Position Determination

Derived from [38], the MLE of the unknown parameter vector can be estimated by examining the likelihood ratio, with H1 corresponding to the target-present hypothesis. When the distribution characteristics of noise n¯l are known, MLE can be used [21, 39]. Here, we assume that n¯l is the identically distributed zero-mean complex Gaussian white noise, with correlation matrix Rl=σl2I. In addition, the receiver stations are distributed widely enough to warrant mutual independence of noise vectors,
(6)nl¯⊥nl′¯,l≠l′.According to Equation (4), we can construct the probability density function (PDF) as
(7)Pr¯H1;p,t0=∏l=1LPrl¯H1;p,t0=C1∏l=1Lexp−12rl¯−blsl¯HRl−1rl¯−blsl¯.where C1 is a constant, and .H represents the operation of conjugation transpose. Thus, the likelihood ratio function can be obtained as
(8)Lr¯;p,t0=Pr¯H1;p,t0Pr¯H0=C1C0exp−12rl¯−blsl¯HRl−1rl¯−blsl¯+12r¯lHRl−1r¯l.Let∂Lrl/∂bl=0, bl=sl¯HRl−1sl¯/sl¯HRl−1sl¯, Rl=σ12I, and s¯2=1. Finally, we can get the maximum likelihood function as follows:
(9)Lr¯;p,t0∝∑l=1L1σl2sl¯Hr¯l2=∑l=1L1σl2∑k=0Ns−1ejωkτlps¯∗kr¯lk2.Remarkably, we observe from Equation (9) that the parameters to be estimated in the objective function include the transmitter position, the signal parameters, and the time t0. It is a highly nonlinear optimization problem whose closed-form description of the analytic optimum solution is difficult to obtain. The numerical solution based on grid search for the region of interest is an alterative. However, the high-dimensional grid search will become computationally prohibitive when the needed localization accuracy increases. If the SNR is lower than a threshold point, the estimation errors rise quickly, i.e., the “threshold effect.” Based on the above arguments, a novel algorithm is investigated in the subsequent sections to trade off algorithm performance for implementation complexity.


## 4. Algorithm Description

In this section, we present a computationally efficient approach to solve a highly nonlinear optimization problem in passive localization. First, we employ STFT technique to extract the valid spectrum fragment of the transmitting signals and then further estimate the instantaneous frequency, spatial power, and other signal parameters. After extracting peak points of the power spectrum, aiming at the disadvantages of low precision of two-step method and high complexity of DPD algorithm at low SNR, we adopt corresponding objective functions to extract the emitter positions. The algorithm is described in detail below.


### 4.1. Estimation of Signal Parameters and Spectrum Detection

Under the circumstance of low SNR, the localization performance of DPD algorithm will deteriorate significantly. In order to improve the localization performance, in this subsection, we adopt the time-frequency analysis method to estimate the relevant parameters and detect the spectrum.


#### 4.1.1. Parameter Estimation Based on Time Frequency Analysis

To estimate the LFM signal parameters, the most important step is to estimate instantaneous frequency (IF) from the received signal. The IF expression of the LFM signal is as follows:(10)flfmIFt≜f0+kt,0≤t≤Tp.After sampling,flfmIFt can be written in a vector form as
(11)flfmIF≜flfmIFt0,⋯,flfmIFtα,⋯,flfmIFtNp−1,where 0≤α≤Np−1, and Np is the number of samples of the LFM signal with Np=Tp∗fs.As can be seen from the above equations, IF is formed by linear combination of initial frequencyf0 and chirp rate k. A large number of IF samples can be obtained through estimation, and the corresponding f0 and k can be estimated from the IF samples by means of parameter fitting.For the LFM signal, short-time Fourier transform (STFT) performs analysis on the signal in the time window and can intuitively reflect the trend of the frequency and the change of spectrum along with time, thus the initial frequency and chirp rate of the signal can be estimated efficiently, which has a good effect on the estimation accuracy, especially at low SNR [40–45]. The common STFT is defined as
(12)Hfτ,f=∫sτht−τexpj2πfτdτ,where hτ is the time window function, which is called analysis window as well.Substituting Equation (12) into the expression of the LFM signal [42], then the STFT for the LFM signal st is defined as
(13)STFTst,f=exp−j2πf−f0+kττ×exp−jπkτ2hf−f0−kτ,0≤τ≤Tp,where
(14)hf≜∑αexp−jπkτ2gτexp−j2πfτ.By taking the square of the module ofhf, the spectrum (SP) of the LFM signal can be written in the following form:
(15)Pτ,f≜STFTst,f2=hf−f0−kτ2.The Gaussian window is chosen as the analysis window so thathf2 is a Gaussian function in essence. hf2 is symmetry about the line f=0 and achieves the maximum value at f=0. The above formula shows that the LFM signal’s spectrum STFTst,f obtains the maximum value at frequency point f0+kτ at each discrete time τ. Therefore, the IF of the LFM signal can be acquired by finding the max point at each discrete time τ.Hence, we can get the IF estimation as the following expression:(16)flfmIFτ=argmaxPτ,f.fFigure2 shows the real value and estimated value of the instantaneous frequency after linearly fitting by the least square method when SNR=−5dB. Thus, the estimated values of chirp rate and initial frequency can be calculated from the fitted straight line.Figure 2
The IF estimated value and real value of LFM signal whenSNR=−5dB.


#### 4.1.2. Extraction of the SP of the LFM Signal

The SP of the observed signal of thelth receiver station can be expressed as [38]
(17)Pt,f=Pt,f,t0<t<t0+τlp,nt,f,else,where nt,f represents the SP of noise.After adding the moving window functionhτ, the extracted LFM signal can be written as
(18)Plfmr,p,t0,Tp=1L∑l=1LPlfmlr,p,t0,Tp.The SP of the LFM signal observed and extracted by thelth receiver station in Equation (17) is shown in Figure 3.Figure 3
SP of the observed and extracted LFM signal of thelth base station.On the basis of the above process, the 2D distribution of the chirp signal’s STFT is converted into 1D relationship between frequency and time, which displays a coarse line in the time-frequency plane. It can be seen from Figure3 that when SNR is at -5 dB, the spectral intensity of the LFM signal is much greater than that of the noise. It can be clearly seen from the figure that the time-frequency characteristics of the LFM signal reflect the sum of power spectral in each window. Summing the maximum value of spectrum of each window helps to estimate the parameters of the LFM signal and realize the spectrum detection. The localization algorithms after spectrum detection will be introduced in the subsequent sections.


### 4.2. The Localization Algorithm Based on Spectrum Detection

According to the “threshold effect” of DPD algorithm, the localization accuracy will deteriorate sharply when the SNR is lower than -5 dB [21]. The mean (or root mean square) of miss distance will also increase rapidly. At high SNR, DPD, Chan, Taylor algorithms have similar performance in localization, but when the SNR is lower than the threshold point, DPD methods have been demonstrated to outperform two-step methods. For most of the DPD algorithms require grid search in the final, which increases computation complexity and transmission demand in the case of lacking prior information. Therefore, Chan and Taylor algorithms are first used for coarse location. Then, grid search is carried out on the basis of the results of the coarse location, which can greatly reduce the workload. The complexity analysis will be detailed in the next subsection. The whole FDPD algorithm process is as follows.


#### 4.2.1. Localization at High SNR

According to [46], in passive multistation localization, after obtaining TDOA measurements by Chan algorithm [13, 36], the distance differences between the emitter and the base stations are obtained, thus multiple TDOA measurements form a set of estimated equations about the emitter. We assume x^,y^ is the estimated position of the emitter and xl,yl,l=1,2,⋯,L is the known position of the lth receiver station, then the emitter position can be calculated by the following formula:
(19)x^y^=x21y21x31y31−1×12K2−K1−R212−R21R112K3−K1−R312−R31R1,where
(20)xl1=xl−x1,yl1=yl−y1,Rl1=Rl−R1=xl−x∧2+yl−y∧2−x1−x∧2+y1−y∧2,Kl=xl2+yl2.It has been shown that the iterative algorithms (Standard least squares, hyperbolic least squares, etc.) require a definition of the initial parameters, which can significantly affect solution accuracy and computational time [47]. Therefore, to ensure the solution accuracy, the position calculated by Chan algorithm is selected for initial position x0,y0 of MS for Taylor expansion of the Taylor algorithm [37]. Then, the weighted least square (WLS) method [48, 49] is used to obtain the least squares estimation of Z as follows:
(21)Z=GTφ−1G−1GTφ−1h,where .T represents the transpose operator, and φ is the covariance matrix of the measurement value of TDOA. The value of the TDOA measurements covariance matrix is estimated under the assumption that there are large numbers of samples, or equivalently, long observation time. Therefore, the Cramer-Rao matrix bound (CRMB) of the estimated TDOA vector is used to represent the value of the matrix [50].
(22)φ=2Tp2π∫0Bω2Sω21+SωtrNω−1×trNω−1Npω−1−Npω−111TNPω−1dω−1,Z=P0R1T,h=R21−R2−R1R31−R3−R1⋮Rl1−Rl−R1,G=x1−xR1−x2−xR2y1−yR1−y2−yR2x1−xR1−x3−xR3y1−yR1−y3−yR3⋮⋮x1−xR1−xl−xRly1−yR1−yl−yRl.where B is the bandwidth of the signal and tr∗ is the trace of matrix ∗. Sω is the signal power spectrum, Nω=diagN1ω,N2ω,⋯,NMω is the noise power spectral matrix, Npω is the lower right M−1 by M−1 partition of the matrix Nω, and 1 is a vector of unity which has the same size as Npω.The selection of the initial value of iterative Taylor-series expansion method has a great influence on the location result. If the initial value is not appropriate, the Taylor algorithm will not converge. On the other hand, Chan algorithm has fast operation speed and low computation complexity. Although the Chan algorithm reduces the localization accuracy in poor channel environment, the localization results can still reflect the general characteristics between the position of the radiation source and the measured value, which is conducive to the convergence of Taylor algorithm. Therefore, we choose Chan algorithm to calculate the TDOA measured values, and the localization results of Chan algorithm are taken as the initial value of Taylor series expansion algorithm, then several iterations are carried out.Hence, this paper adopts the collaborative localization of Chan and Taylor algorithms using Equations (19)-(22) at high SNR, which can not only ensure the small error of the initial value of the Taylor expansion, so as to ensure the localization accuracy, but also greatly reduce the amount of calculation.


#### 4.2.2. Localization at Low SNR

When the SNR is low, the DPD algorithm is used on the basis of the location results of Chan and Taylor algorithms. To get more accurate waveform parameters, STFT is used to assist signal parameter estimation. DPD algorithm can obtain better localization accuracy on the basis of more accurate signal parameters estimation. In order to obtain more accurate positions, here, an objective function is constructed for grid search.The IF of the LFM signal is estimated from Equation (16), so as to estimate the initial frequency and chirp rate, and then substituted into the following function, which describes the fitting error per unit time at different positions, window width, and starting time, the function is given by
(23)Fr,p,t0,Tp=∑αflfmIF∧τ−kτ−f02Tp,where 0≤α≤Ns−1, Ns is the number of samples. ∑αflfmIF∧τ−kτ−f02 is referred to the minimum cost function, thus Fr,p,t0,Tp denotes the mean minimum cost function, i.e., the fitting error per unit time at different positions. From Equation (15), 1/L∑l=1L1/σl2sl¯Hrl¯2 represents the total likelihood ratio function of the signal, which reflects the sum of spectral power in all windows. ∑α∑l=1LmaxfPτ,fhτ reflects the summation of the maximum SP within the sliding window. To obtain higher parameter estimation accuracy, the following objective function is defined:
(24)Lr¯;p,t0=1L∑l=1L1σl2sl¯Hr¯l2×∑α∑l=1LmaxfPτ,fhτFr,p,t0,Tp.Compared with the function in Equation (9), the objective function includes more factors of indicating the parameter estimation by means of STFT, which is maximized by Fr,p,t0,Tp. Then through grid search, we can estimate the emitter position p, initial time t0, and signal pulse width Tp more accurately. In most instances, four-dimensional search cannot be considered as a viable solution to solve a highly nonlinear optimization problem. Therefore, we used an efficient decoupled strategy to solve the problem. The basic idea of the strategy is to turn the four-dimensional search into a two-dimensional parameter estimation combined with grid search. Specifically, based on STFT, a two-dimensional estimator is firstly designed to extract the initial frequency f0 and chirp rate k of the LFM signal in the time-frequency domain. Then, a four-dimensional search, including the previously estimated waveform parameters, is carried out, and the position of the emitter, the corresponding transmitted time, time duration, and initial frequency can be obtained by the search.For determining when and which method should be used, that is, to decide whether the emitter signal is in high or low SNR environment, Equation (15) in 4.1.1 can be utilized. By means of maximizing the spectrum at each window, it can be seen that the maximum spectrum of the LFM signal is obviously different from the maximum spectrum of the noise. The spectrum difference between the LFM signal and the noise can be calculated and then compared with the threshold point to determine whether the emitter signal is at high or low SNR. The simulation diagram of the maximum spectrum value changing over time is shown in Figure 4.Figure 4
The maximum SP in each STFT window with time of -15 dB.The general pseudocodes of the proposed FDPD algorithm are shown in Algorithm1.Algorithm 1: Summary of the FDPD Algorithm.
Input: The observed signal r, the base station coordinates, space of grid search (Nx×Ny) for p;Output: Estimated parameters and the position of the radiation source;1: Forl =1 ··· L do2:   Compute the SP forlth observed signal.3:   Compute the maximum SP of thelth observed signal4:   Compare the maximum SP difference between the LFM5:   Signals and noise signal with the threshold value.6: End for7:If the maximum SP difference > threshold then8:   Compute the emitter positionp using Equation (19)–(22).9:End if10:If the maximum SP difference < threshold then11:   Forl =1 ··· L do12:    Compute the emitter position initialp1 using13:    Equation (19) - (22).14:    Determine the search range (Nx×Ny) by the15:    Error analysis of Chan and Taylor algorithm.16:    Fori=1⋯Nx×Ny do17:     Estimate the IF of thelth observed signal using Equation (16).18:     Estimate the initial frequency and chirp rate.19:     Extract the IF and compute the fitting error every time20:     Window per unit time at different position.21:     Compute the objective function value using Equation (24).22:     Find the maximum of the objective function and compute the23:     Emitter positionp.24:    End for25:   End for26: End if


## 4.1. Estimation of Signal Parameters and Spectrum Detection

Under the circumstance of low SNR, the localization performance of DPD algorithm will deteriorate significantly. In order to improve the localization performance, in this subsection, we adopt the time-frequency analysis method to estimate the relevant parameters and detect the spectrum.


### 4.1.1. Parameter Estimation Based on Time Frequency Analysis

To estimate the LFM signal parameters, the most important step is to estimate instantaneous frequency (IF) from the received signal. The IF expression of the LFM signal is as follows:(10)flfmIFt≜f0+kt,0≤t≤Tp.After sampling,flfmIFt can be written in a vector form as
(11)flfmIF≜flfmIFt0,⋯,flfmIFtα,⋯,flfmIFtNp−1,where 0≤α≤Np−1, and Np is the number of samples of the LFM signal with Np=Tp∗fs.As can be seen from the above equations, IF is formed by linear combination of initial frequencyf0 and chirp rate k. A large number of IF samples can be obtained through estimation, and the corresponding f0 and k can be estimated from the IF samples by means of parameter fitting.For the LFM signal, short-time Fourier transform (STFT) performs analysis on the signal in the time window and can intuitively reflect the trend of the frequency and the change of spectrum along with time, thus the initial frequency and chirp rate of the signal can be estimated efficiently, which has a good effect on the estimation accuracy, especially at low SNR [40–45]. The common STFT is defined as
(12)Hfτ,f=∫sτht−τexpj2πfτdτ,where hτ is the time window function, which is called analysis window as well.Substituting Equation (12) into the expression of the LFM signal [42], then the STFT for the LFM signal st is defined as
(13)STFTst,f=exp−j2πf−f0+kττ×exp−jπkτ2hf−f0−kτ,0≤τ≤Tp,where
(14)hf≜∑αexp−jπkτ2gτexp−j2πfτ.By taking the square of the module ofhf, the spectrum (SP) of the LFM signal can be written in the following form:
(15)Pτ,f≜STFTst,f2=hf−f0−kτ2.The Gaussian window is chosen as the analysis window so thathf2 is a Gaussian function in essence. hf2 is symmetry about the line f=0 and achieves the maximum value at f=0. The above formula shows that the LFM signal’s spectrum STFTst,f obtains the maximum value at frequency point f0+kτ at each discrete time τ. Therefore, the IF of the LFM signal can be acquired by finding the max point at each discrete time τ.Hence, we can get the IF estimation as the following expression:(16)flfmIFτ=argmaxPτ,f.fFigure2 shows the real value and estimated value of the instantaneous frequency after linearly fitting by the least square method when SNR=−5dB. Thus, the estimated values of chirp rate and initial frequency can be calculated from the fitted straight line.Figure 2
The IF estimated value and real value of LFM signal whenSNR=−5dB.


### 4.1.2. Extraction of the SP of the LFM Signal

The SP of the observed signal of thelth receiver station can be expressed as [38]
(17)Pt,f=Pt,f,t0<t<t0+τlp,nt,f,else,where nt,f represents the SP of noise.After adding the moving window functionhτ, the extracted LFM signal can be written as
(18)Plfmr,p,t0,Tp=1L∑l=1LPlfmlr,p,t0,Tp.The SP of the LFM signal observed and extracted by thelth receiver station in Equation (17) is shown in Figure 3.Figure 3
SP of the observed and extracted LFM signal of thelth base station.On the basis of the above process, the 2D distribution of the chirp signal’s STFT is converted into 1D relationship between frequency and time, which displays a coarse line in the time-frequency plane. It can be seen from Figure3 that when SNR is at -5 dB, the spectral intensity of the LFM signal is much greater than that of the noise. It can be clearly seen from the figure that the time-frequency characteristics of the LFM signal reflect the sum of power spectral in each window. Summing the maximum value of spectrum of each window helps to estimate the parameters of the LFM signal and realize the spectrum detection. The localization algorithms after spectrum detection will be introduced in the subsequent sections.


## 4.1.1. Parameter Estimation Based on Time Frequency Analysis

To estimate the LFM signal parameters, the most important step is to estimate instantaneous frequency (IF) from the received signal. The IF expression of the LFM signal is as follows:(10)flfmIFt≜f0+kt,0≤t≤Tp.After sampling,flfmIFt can be written in a vector form as
(11)flfmIF≜flfmIFt0,⋯,flfmIFtα,⋯,flfmIFtNp−1,where 0≤α≤Np−1, and Np is the number of samples of the LFM signal with Np=Tp∗fs.As can be seen from the above equations, IF is formed by linear combination of initial frequencyf0 and chirp rate k. A large number of IF samples can be obtained through estimation, and the corresponding f0 and k can be estimated from the IF samples by means of parameter fitting.For the LFM signal, short-time Fourier transform (STFT) performs analysis on the signal in the time window and can intuitively reflect the trend of the frequency and the change of spectrum along with time, thus the initial frequency and chirp rate of the signal can be estimated efficiently, which has a good effect on the estimation accuracy, especially at low SNR [40–45]. The common STFT is defined as
(12)Hfτ,f=∫sτht−τexpj2πfτdτ,where hτ is the time window function, which is called analysis window as well.Substituting Equation (12) into the expression of the LFM signal [42], then the STFT for the LFM signal st is defined as
(13)STFTst,f=exp−j2πf−f0+kττ×exp−jπkτ2hf−f0−kτ,0≤τ≤Tp,where
(14)hf≜∑αexp−jπkτ2gτexp−j2πfτ.By taking the square of the module ofhf, the spectrum (SP) of the LFM signal can be written in the following form:
(15)Pτ,f≜STFTst,f2=hf−f0−kτ2.The Gaussian window is chosen as the analysis window so thathf2 is a Gaussian function in essence. hf2 is symmetry about the line f=0 and achieves the maximum value at f=0. The above formula shows that the LFM signal’s spectrum STFTst,f obtains the maximum value at frequency point f0+kτ at each discrete time τ. Therefore, the IF of the LFM signal can be acquired by finding the max point at each discrete time τ.Hence, we can get the IF estimation as the following expression:(16)flfmIFτ=argmaxPτ,f.fFigure2 shows the real value and estimated value of the instantaneous frequency after linearly fitting by the least square method when SNR=−5dB. Thus, the estimated values of chirp rate and initial frequency can be calculated from the fitted straight line.Figure 2
The IF estimated value and real value of LFM signal whenSNR=−5dB.


## 4.1.2. Extraction of the SP of the LFM Signal

The SP of the observed signal of thelth receiver station can be expressed as [38]
(17)Pt,f=Pt,f,t0<t<t0+τlp,nt,f,else,where nt,f represents the SP of noise.After adding the moving window functionhτ, the extracted LFM signal can be written as
(18)Plfmr,p,t0,Tp=1L∑l=1LPlfmlr,p,t0,Tp.The SP of the LFM signal observed and extracted by thelth receiver station in Equation (17) is shown in Figure 3.Figure 3
SP of the observed and extracted LFM signal of thelth base station.On the basis of the above process, the 2D distribution of the chirp signal’s STFT is converted into 1D relationship between frequency and time, which displays a coarse line in the time-frequency plane. It can be seen from Figure3 that when SNR is at -5 dB, the spectral intensity of the LFM signal is much greater than that of the noise. It can be clearly seen from the figure that the time-frequency characteristics of the LFM signal reflect the sum of power spectral in each window. Summing the maximum value of spectrum of each window helps to estimate the parameters of the LFM signal and realize the spectrum detection. The localization algorithms after spectrum detection will be introduced in the subsequent sections.


## 4.2. The Localization Algorithm Based on Spectrum Detection

According to the “threshold effect” of DPD algorithm, the localization accuracy will deteriorate sharply when the SNR is lower than -5 dB [21]. The mean (or root mean square) of miss distance will also increase rapidly. At high SNR, DPD, Chan, Taylor algorithms have similar performance in localization, but when the SNR is lower than the threshold point, DPD methods have been demonstrated to outperform two-step methods. For most of the DPD algorithms require grid search in the final, which increases computation complexity and transmission demand in the case of lacking prior information. Therefore, Chan and Taylor algorithms are first used for coarse location. Then, grid search is carried out on the basis of the results of the coarse location, which can greatly reduce the workload. The complexity analysis will be detailed in the next subsection. The whole FDPD algorithm process is as follows.


### 4.2.1. Localization at High SNR

According to [46], in passive multistation localization, after obtaining TDOA measurements by Chan algorithm [13, 36], the distance differences between the emitter and the base stations are obtained, thus multiple TDOA measurements form a set of estimated equations about the emitter. We assume x^,y^ is the estimated position of the emitter and xl,yl,l=1,2,⋯,L is the known position of the lth receiver station, then the emitter position can be calculated by the following formula:
(19)x^y^=x21y21x31y31−1×12K2−K1−R212−R21R112K3−K1−R312−R31R1,where
(20)xl1=xl−x1,yl1=yl−y1,Rl1=Rl−R1=xl−x∧2+yl−y∧2−x1−x∧2+y1−y∧2,Kl=xl2+yl2.It has been shown that the iterative algorithms (Standard least squares, hyperbolic least squares, etc.) require a definition of the initial parameters, which can significantly affect solution accuracy and computational time [47]. Therefore, to ensure the solution accuracy, the position calculated by Chan algorithm is selected for initial position x0,y0 of MS for Taylor expansion of the Taylor algorithm [37]. Then, the weighted least square (WLS) method [48, 49] is used to obtain the least squares estimation of Z as follows:
(21)Z=GTφ−1G−1GTφ−1h,where .T represents the transpose operator, and φ is the covariance matrix of the measurement value of TDOA. The value of the TDOA measurements covariance matrix is estimated under the assumption that there are large numbers of samples, or equivalently, long observation time. Therefore, the Cramer-Rao matrix bound (CRMB) of the estimated TDOA vector is used to represent the value of the matrix [50].
(22)φ=2Tp2π∫0Bω2Sω21+SωtrNω−1×trNω−1Npω−1−Npω−111TNPω−1dω−1,Z=P0R1T,h=R21−R2−R1R31−R3−R1⋮Rl1−Rl−R1,G=x1−xR1−x2−xR2y1−yR1−y2−yR2x1−xR1−x3−xR3y1−yR1−y3−yR3⋮⋮x1−xR1−xl−xRly1−yR1−yl−yRl.where B is the bandwidth of the signal and tr∗ is the trace of matrix ∗. Sω is the signal power spectrum, Nω=diagN1ω,N2ω,⋯,NMω is the noise power spectral matrix, Npω is the lower right M−1 by M−1 partition of the matrix Nω, and 1 is a vector of unity which has the same size as Npω.The selection of the initial value of iterative Taylor-series expansion method has a great influence on the location result. If the initial value is not appropriate, the Taylor algorithm will not converge. On the other hand, Chan algorithm has fast operation speed and low computation complexity. Although the Chan algorithm reduces the localization accuracy in poor channel environment, the localization results can still reflect the general characteristics between the position of the radiation source and the measured value, which is conducive to the convergence of Taylor algorithm. Therefore, we choose Chan algorithm to calculate the TDOA measured values, and the localization results of Chan algorithm are taken as the initial value of Taylor series expansion algorithm, then several iterations are carried out.Hence, this paper adopts the collaborative localization of Chan and Taylor algorithms using Equations (19)-(22) at high SNR, which can not only ensure the small error of the initial value of the Taylor expansion, so as to ensure the localization accuracy, but also greatly reduce the amount of calculation.


### 4.2.2. Localization at Low SNR

When the SNR is low, the DPD algorithm is used on the basis of the location results of Chan and Taylor algorithms. To get more accurate waveform parameters, STFT is used to assist signal parameter estimation. DPD algorithm can obtain better localization accuracy on the basis of more accurate signal parameters estimation. In order to obtain more accurate positions, here, an objective function is constructed for grid search.The IF of the LFM signal is estimated from Equation (16), so as to estimate the initial frequency and chirp rate, and then substituted into the following function, which describes the fitting error per unit time at different positions, window width, and starting time, the function is given by
(23)Fr,p,t0,Tp=∑αflfmIF∧τ−kτ−f02Tp,where 0≤α≤Ns−1, Ns is the number of samples. ∑αflfmIF∧τ−kτ−f02 is referred to the minimum cost function, thus Fr,p,t0,Tp denotes the mean minimum cost function, i.e., the fitting error per unit time at different positions. From Equation (15), 1/L∑l=1L1/σl2sl¯Hrl¯2 represents the total likelihood ratio function of the signal, which reflects the sum of spectral power in all windows. ∑α∑l=1LmaxfPτ,fhτ reflects the summation of the maximum SP within the sliding window. To obtain higher parameter estimation accuracy, the following objective function is defined:
(24)Lr¯;p,t0=1L∑l=1L1σl2sl¯Hr¯l2×∑α∑l=1LmaxfPτ,fhτFr,p,t0,Tp.Compared with the function in Equation (9), the objective function includes more factors of indicating the parameter estimation by means of STFT, which is maximized by Fr,p,t0,Tp. Then through grid search, we can estimate the emitter position p, initial time t0, and signal pulse width Tp more accurately. In most instances, four-dimensional search cannot be considered as a viable solution to solve a highly nonlinear optimization problem. Therefore, we used an efficient decoupled strategy to solve the problem. The basic idea of the strategy is to turn the four-dimensional search into a two-dimensional parameter estimation combined with grid search. Specifically, based on STFT, a two-dimensional estimator is firstly designed to extract the initial frequency f0 and chirp rate k of the LFM signal in the time-frequency domain. Then, a four-dimensional search, including the previously estimated waveform parameters, is carried out, and the position of the emitter, the corresponding transmitted time, time duration, and initial frequency can be obtained by the search.For determining when and which method should be used, that is, to decide whether the emitter signal is in high or low SNR environment, Equation (15) in 4.1.1 can be utilized. By means of maximizing the spectrum at each window, it can be seen that the maximum spectrum of the LFM signal is obviously different from the maximum spectrum of the noise. The spectrum difference between the LFM signal and the noise can be calculated and then compared with the threshold point to determine whether the emitter signal is at high or low SNR. The simulation diagram of the maximum spectrum value changing over time is shown in Figure 4.Figure 4
The maximum SP in each STFT window with time of -15 dB.The general pseudocodes of the proposed FDPD algorithm are shown in Algorithm1.Algorithm 1: Summary of the FDPD Algorithm.
Input: The observed signal r, the base station coordinates, space of grid search (Nx×Ny) for p;Output: Estimated parameters and the position of the radiation source;1: Forl =1 ··· L do2:   Compute the SP forlth observed signal.3:   Compute the maximum SP of thelth observed signal4:   Compare the maximum SP difference between the LFM5:   Signals and noise signal with the threshold value.6: End for7:If the maximum SP difference > threshold then8:   Compute the emitter positionp using Equation (19)–(22).9:End if10:If the maximum SP difference < threshold then11:   Forl =1 ··· L do12:    Compute the emitter position initialp1 using13:    Equation (19) - (22).14:    Determine the search range (Nx×Ny) by the15:    Error analysis of Chan and Taylor algorithm.16:    Fori=1⋯Nx×Ny do17:     Estimate the IF of thelth observed signal using Equation (16).18:     Estimate the initial frequency and chirp rate.19:     Extract the IF and compute the fitting error every time20:     Window per unit time at different position.21:     Compute the objective function value using Equation (24).22:     Find the maximum of the objective function and compute the23:     Emitter positionp.24:    End for25:   End for26: End if


## 4.2.1. Localization at High SNR

According to [46], in passive multistation localization, after obtaining TDOA measurements by Chan algorithm [13, 36], the distance differences between the emitter and the base stations are obtained, thus multiple TDOA measurements form a set of estimated equations about the emitter. We assume x^,y^ is the estimated position of the emitter and xl,yl,l=1,2,⋯,L is the known position of the lth receiver station, then the emitter position can be calculated by the following formula:
(19)x^y^=x21y21x31y31−1×12K2−K1−R212−R21R112K3−K1−R312−R31R1,where
(20)xl1=xl−x1,yl1=yl−y1,Rl1=Rl−R1=xl−x∧2+yl−y∧2−x1−x∧2+y1−y∧2,Kl=xl2+yl2.It has been shown that the iterative algorithms (Standard least squares, hyperbolic least squares, etc.) require a definition of the initial parameters, which can significantly affect solution accuracy and computational time [47]. Therefore, to ensure the solution accuracy, the position calculated by Chan algorithm is selected for initial position x0,y0 of MS for Taylor expansion of the Taylor algorithm [37]. Then, the weighted least square (WLS) method [48, 49] is used to obtain the least squares estimation of Z as follows:
(21)Z=GTφ−1G−1GTφ−1h,where .T represents the transpose operator, and φ is the covariance matrix of the measurement value of TDOA. The value of the TDOA measurements covariance matrix is estimated under the assumption that there are large numbers of samples, or equivalently, long observation time. Therefore, the Cramer-Rao matrix bound (CRMB) of the estimated TDOA vector is used to represent the value of the matrix [50].
(22)φ=2Tp2π∫0Bω2Sω21+SωtrNω−1×trNω−1Npω−1−Npω−111TNPω−1dω−1,Z=P0R1T,h=R21−R2−R1R31−R3−R1⋮Rl1−Rl−R1,G=x1−xR1−x2−xR2y1−yR1−y2−yR2x1−xR1−x3−xR3y1−yR1−y3−yR3⋮⋮x1−xR1−xl−xRly1−yR1−yl−yRl.where B is the bandwidth of the signal and tr∗ is the trace of matrix ∗. Sω is the signal power spectrum, Nω=diagN1ω,N2ω,⋯,NMω is the noise power spectral matrix, Npω is the lower right M−1 by M−1 partition of the matrix Nω, and 1 is a vector of unity which has the same size as Npω.The selection of the initial value of iterative Taylor-series expansion method has a great influence on the location result. If the initial value is not appropriate, the Taylor algorithm will not converge. On the other hand, Chan algorithm has fast operation speed and low computation complexity. Although the Chan algorithm reduces the localization accuracy in poor channel environment, the localization results can still reflect the general characteristics between the position of the radiation source and the measured value, which is conducive to the convergence of Taylor algorithm. Therefore, we choose Chan algorithm to calculate the TDOA measured values, and the localization results of Chan algorithm are taken as the initial value of Taylor series expansion algorithm, then several iterations are carried out.Hence, this paper adopts the collaborative localization of Chan and Taylor algorithms using Equations (19)-(22) at high SNR, which can not only ensure the small error of the initial value of the Taylor expansion, so as to ensure the localization accuracy, but also greatly reduce the amount of calculation.


## 4.2.2. Localization at Low SNR

When the SNR is low, the DPD algorithm is used on the basis of the location results of Chan and Taylor algorithms. To get more accurate waveform parameters, STFT is used to assist signal parameter estimation. DPD algorithm can obtain better localization accuracy on the basis of more accurate signal parameters estimation. In order to obtain more accurate positions, here, an objective function is constructed for grid search.The IF of the LFM signal is estimated from Equation (16), so as to estimate the initial frequency and chirp rate, and then substituted into the following function, which describes the fitting error per unit time at different positions, window width, and starting time, the function is given by
(23)Fr,p,t0,Tp=∑αflfmIF∧τ−kτ−f02Tp,where 0≤α≤Ns−1, Ns is the number of samples. ∑αflfmIF∧τ−kτ−f02 is referred to the minimum cost function, thus Fr,p,t0,Tp denotes the mean minimum cost function, i.e., the fitting error per unit time at different positions. From Equation (15), 1/L∑l=1L1/σl2sl¯Hrl¯2 represents the total likelihood ratio function of the signal, which reflects the sum of spectral power in all windows. ∑α∑l=1LmaxfPτ,fhτ reflects the summation of the maximum SP within the sliding window. To obtain higher parameter estimation accuracy, the following objective function is defined:
(24)Lr¯;p,t0=1L∑l=1L1σl2sl¯Hr¯l2×∑α∑l=1LmaxfPτ,fhτFr,p,t0,Tp.Compared with the function in Equation (9), the objective function includes more factors of indicating the parameter estimation by means of STFT, which is maximized by Fr,p,t0,Tp. Then through grid search, we can estimate the emitter position p, initial time t0, and signal pulse width Tp more accurately. In most instances, four-dimensional search cannot be considered as a viable solution to solve a highly nonlinear optimization problem. Therefore, we used an efficient decoupled strategy to solve the problem. The basic idea of the strategy is to turn the four-dimensional search into a two-dimensional parameter estimation combined with grid search. Specifically, based on STFT, a two-dimensional estimator is firstly designed to extract the initial frequency f0 and chirp rate k of the LFM signal in the time-frequency domain. Then, a four-dimensional search, including the previously estimated waveform parameters, is carried out, and the position of the emitter, the corresponding transmitted time, time duration, and initial frequency can be obtained by the search.For determining when and which method should be used, that is, to decide whether the emitter signal is in high or low SNR environment, Equation (15) in 4.1.1 can be utilized. By means of maximizing the spectrum at each window, it can be seen that the maximum spectrum of the LFM signal is obviously different from the maximum spectrum of the noise. The spectrum difference between the LFM signal and the noise can be calculated and then compared with the threshold point to determine whether the emitter signal is at high or low SNR. The simulation diagram of the maximum spectrum value changing over time is shown in Figure 4.Figure 4
The maximum SP in each STFT window with time of -15 dB.The general pseudocodes of the proposed FDPD algorithm are shown in Algorithm1.Algorithm 1: Summary of the FDPD Algorithm.
Input: The observed signal r, the base station coordinates, space of grid search (Nx×Ny) for p;Output: Estimated parameters and the position of the radiation source;1: Forl =1 ··· L do2:   Compute the SP forlth observed signal.3:   Compute the maximum SP of thelth observed signal4:   Compare the maximum SP difference between the LFM5:   Signals and noise signal with the threshold value.6: End for7:If the maximum SP difference > threshold then8:   Compute the emitter positionp using Equation (19)–(22).9:End if10:If the maximum SP difference < threshold then11:   Forl =1 ··· L do12:    Compute the emitter position initialp1 using13:    Equation (19) - (22).14:    Determine the search range (Nx×Ny) by the15:    Error analysis of Chan and Taylor algorithm.16:    Fori=1⋯Nx×Ny do17:     Estimate the IF of thelth observed signal using Equation (16).18:     Estimate the initial frequency and chirp rate.19:     Extract the IF and compute the fitting error every time20:     Window per unit time at different position.21:     Compute the objective function value using Equation (24).22:     Find the maximum of the objective function and compute the23:     Emitter positionp.24:    End for25:   End for26: End if


## 5. Simulation Results

In this section, simulations are brought forward to verify the performance of the proposed algorithm. Four receiver stations used in the simulation are located at (0,0)m (m denotes meter), (6000,0)m, (6000,6000)m, and (2400,6000)m, respectively. The approximate coordinates of the radiation source are (6000,4500)m. We perform extensive Monte Carlo simulations. In the following analysis, the results are collected and processed by averaging M Monte Carlo experiments. Simulation parameters are shown in Table1, where M is the number of Monte Carlo experiments.Table 1
Simulation signal parameters.
LFM signal parameterParameter valuesPulse width20μsBand width20 MHzSampling frequency100 MHzInitial frequency10 MHzChirp rate1 MHz/μsFor each Monte Carlo experiment, we assume that the error between the estimated positionx^l,y^l of the emitter and the actual position x,y is calculated as follows:
(25)pl=xl∧−x2+y∧l−y2,where pi is the estimation of p of the ith Monte Carlo trial.The mean error (ME) of the estimated position is obtained by the statistics of M Monte Carlo experiments. ME is defined as(26)ME=1M∑lMpl.The root mean square error (RMSE) of the estimated parameters is defined as(27)RMSE=1/M∑iMθ∧l−θ2θ×100%.where θ denotes the parameter to be estimated, such as chirp rate and initial frequency.


### 5.1. Analysis of Simulation Results

In Figure3, we extract the 1D frequency-time curve of the LFM signal based on STFT and perform STFT on the transmitted signalrlt, the spectrum of which is symmetry aboutf=0. In order to get the SNR for spectrum detection, according to 1/L∑α∑l=1LmaxfPτ,fhτ in Equation (18), the amount of maximum SP within the sliding window can be viewed as a factor of SP. The simulation of the variation of the maximum SP in each STFT window with time is shown in the following figure.It can be seen from Figure4 that when the SNR is -15 dB, the maximum SP value of the LFM signal is obviously different from the maximum SP value of the noise signal. Then, the SP differential is calculated and compared with the threshold point to determine the signal environment.Figures5 and 6 are simulation comparison diagrams of FDPD and common direct localization algorithms.Figure 5
The MEs of the estimated emitter position for four different algorithms.Figure 6
The RMSEs of the estimated emitter position for four different algorithms.As can be seen from Figures5 and 6, the DPD-known algorithm can acquire the optimal localization performance for the case that the transmitted signal is completely known to the receivers. Therefore, the performance of DPD-known can be viewed as an upper bound for other algorithms. Moreover, in high SNR environment, namely, higher than 0 dB, all the four algorithms achieve superior localization accuracy and perform nearly equally. When the SNR is below the threshold point, the ME and RMSE of these algorithms rise quickly as SNR decreases. As SNR decreases, the accuracy of the DPD-unknown deteriorates rapidly. When the emitter signal is unknown, FDPD and DPD-STFT-w algorithm has better localization accuracy compared with DPD-unknown algorithm, which can ensure the accuracy that the RMSE is less than 100 m when the SNR is above -5 dB. However, considering the extremely high computational complexity of DPD algorithms (DPD-STFT, EM-DPD, etc.), Chan and Taylor algorithms are used for coarse localization, and parameters are estimated on the basis of them, then the objective function is maximized by using the estimated parameters which are more accurate. This will greatly reduce the amount of calculation, so it is called fast DPD (FDPD) algorithm. The complexity analysis of FDPD algorithm and DPD, DPD-known and DPD-STFT-w algorithms will be described in detail in the next section.


### 5.2. Analysis of Complexity

The complexity of Chan and Taylor algorithms isOn2L−1, where n is the characteristic dimension. In this paper, n is equals to 2. L denotes the number of receiver stations. When L is 4, the complexity is O12, which is relatively low under the simulation background of this paper.The complexity of DPD and DPD-STFT-w algorithm is composed of the complexity of two-dimensional search, complexity of maximum likelihood estimation, and least square iterative fitting. The complexity of maximum likelihood estimation least square iterative fitting isOfL,Ns and OLNsNs+1, respectively, both of which increase with the number of receiver stations L and the number of samples Ns. Hence, the overall computational complexity of DPD algorithm is ONx1Ny1NpNτLNsNs+1+OfL,Ns, where Nx1 and Ny1 denote the number of meshes on the x-axis and y-axis, respectively. Compared with DPD algorithm, DPD-STFT-w algorithm needs to take STFT of the observed signal from every receiver station, whose computation complexity is LNSτ/Δτ. Therefore, the overall computational complexity of DPD-STFT-w algorithm is OLNsτ/Δτ+Nx1Ny1NpNτLNsNs+1+OfL,Ns, where τ is the window width of STFT, Δτ is the gliding step length of STFT, Np is the total number of samples of the LFM signal, and ONto is the computational complexity of the estimation of t0.For the complexity of two-dimensional search is reduced significantly, then the complexity of FDPD algorithm isOLNsτ/Δτ+Nx2Ny2NpNτLNsNs+1+OfL,Ns, where Nx2,Ny2 denotes the number of meshes on the x-axis and y-axis, respectively. It can be seen from the complexity of the three localization algorithms that the complexity is determined by two-dimensional grid search, thus the proposed FDPD method reduces the meshes on the axes through combining parameter estimation, time frequency analysis, STFT, and two-step method with DPD algorithm.The complexity of these algorithms is summarized in Table2.Table 2
Computation complexity of the three different algorithms.
AlgorithmTotal operationsMean runtime (seconds)Chan and TaylorOLNsNs+10.05DPDONx1Ny1NpNτLNsNs+1+OfL,Ns105.11DPD-STFT-wOLNsτΔτ+Nx1Ny1NpNτLNsNs+1+OfL,Ns145.80F-DPDOLNsτΔτ+Nx2Ny2NpNτLNsNs+1+OfL,Ns3.77It can be seen from Table2 that the complexity of FDPD algorithm has been greatly reduced under the condition of similar positioning accuracy. The combination of coarse grid and fine grid is considered in our FDPD algorithm to realize the balance between the number of grid points and the localization accuracy. Therefore, compared with DPD, DPD-STFT-w, and most of the direct localization algorithms, it can save more than 90% of the time.


## 5.1. Analysis of Simulation Results

In Figure3, we extract the 1D frequency-time curve of the LFM signal based on STFT and perform STFT on the transmitted signalrlt, the spectrum of which is symmetry aboutf=0. In order to get the SNR for spectrum detection, according to 1/L∑α∑l=1LmaxfPτ,fhτ in Equation (18), the amount of maximum SP within the sliding window can be viewed as a factor of SP. The simulation of the variation of the maximum SP in each STFT window with time is shown in the following figure.It can be seen from Figure4 that when the SNR is -15 dB, the maximum SP value of the LFM signal is obviously different from the maximum SP value of the noise signal. Then, the SP differential is calculated and compared with the threshold point to determine the signal environment.Figures5 and 6 are simulation comparison diagrams of FDPD and common direct localization algorithms.Figure 5
The MEs of the estimated emitter position for four different algorithms.Figure 6
The RMSEs of the estimated emitter position for four different algorithms.As can be seen from Figures5 and 6, the DPD-known algorithm can acquire the optimal localization performance for the case that the transmitted signal is completely known to the receivers. Therefore, the performance of DPD-known can be viewed as an upper bound for other algorithms. Moreover, in high SNR environment, namely, higher than 0 dB, all the four algorithms achieve superior localization accuracy and perform nearly equally. When the SNR is below the threshold point, the ME and RMSE of these algorithms rise quickly as SNR decreases. As SNR decreases, the accuracy of the DPD-unknown deteriorates rapidly. When the emitter signal is unknown, FDPD and DPD-STFT-w algorithm has better localization accuracy compared with DPD-unknown algorithm, which can ensure the accuracy that the RMSE is less than 100 m when the SNR is above -5 dB. However, considering the extremely high computational complexity of DPD algorithms (DPD-STFT, EM-DPD, etc.), Chan and Taylor algorithms are used for coarse localization, and parameters are estimated on the basis of them, then the objective function is maximized by using the estimated parameters which are more accurate. This will greatly reduce the amount of calculation, so it is called fast DPD (FDPD) algorithm. The complexity analysis of FDPD algorithm and DPD, DPD-known and DPD-STFT-w algorithms will be described in detail in the next section.


## 5.2. Analysis of Complexity

The complexity of Chan and Taylor algorithms isOn2L−1, where n is the characteristic dimension. In this paper, n is equals to 2. L denotes the number of receiver stations. When L is 4, the complexity is O12, which is relatively low under the simulation background of this paper.The complexity of DPD and DPD-STFT-w algorithm is composed of the complexity of two-dimensional search, complexity of maximum likelihood estimation, and least square iterative fitting. The complexity of maximum likelihood estimation least square iterative fitting isOfL,Ns and OLNsNs+1, respectively, both of which increase with the number of receiver stations L and the number of samples Ns. Hence, the overall computational complexity of DPD algorithm is ONx1Ny1NpNτLNsNs+1+OfL,Ns, where Nx1 and Ny1 denote the number of meshes on the x-axis and y-axis, respectively. Compared with DPD algorithm, DPD-STFT-w algorithm needs to take STFT of the observed signal from every receiver station, whose computation complexity is LNSτ/Δτ. Therefore, the overall computational complexity of DPD-STFT-w algorithm is OLNsτ/Δτ+Nx1Ny1NpNτLNsNs+1+OfL,Ns, where τ is the window width of STFT, Δτ is the gliding step length of STFT, Np is the total number of samples of the LFM signal, and ONto is the computational complexity of the estimation of t0.For the complexity of two-dimensional search is reduced significantly, then the complexity of FDPD algorithm isOLNsτ/Δτ+Nx2Ny2NpNτLNsNs+1+OfL,Ns, where Nx2,Ny2 denotes the number of meshes on the x-axis and y-axis, respectively. It can be seen from the complexity of the three localization algorithms that the complexity is determined by two-dimensional grid search, thus the proposed FDPD method reduces the meshes on the axes through combining parameter estimation, time frequency analysis, STFT, and two-step method with DPD algorithm.The complexity of these algorithms is summarized in Table2.Table 2
Computation complexity of the three different algorithms.
AlgorithmTotal operationsMean runtime (seconds)Chan and TaylorOLNsNs+10.05DPDONx1Ny1NpNτLNsNs+1+OfL,Ns105.11DPD-STFT-wOLNsτΔτ+Nx1Ny1NpNτLNsNs+1+OfL,Ns145.80F-DPDOLNsτΔτ+Nx2Ny2NpNτLNsNs+1+OfL,Ns3.77It can be seen from Table2 that the complexity of FDPD algorithm has been greatly reduced under the condition of similar positioning accuracy. The combination of coarse grid and fine grid is considered in our FDPD algorithm to realize the balance between the number of grid points and the localization accuracy. Therefore, compared with DPD, DPD-STFT-w, and most of the direct localization algorithms, it can save more than 90% of the time.


## 6. Conclusions

In this paper, we consider the passive localization of an LFM signal emitter with unknown signal parameters. The problem is formulated as a highly nonlinear optimization problem at the beginning; then, the FDPD approach for jointly estimating the position of the emitter and its transmitted signal parameters is presented. The FDPD method combines two-step method (Chan algorithm, Taylor algorithm) and one-step method (DPD algorithm) with short-time Fourier transform to ensure the localization accuracy and greatly reduces the computational complexity compared with DPD algorithm. Spectrum detection is aimed at improving the “threshold effect” of DPD algorithm. When the SNR is higher than the threshold point, it can achieve almost the same localization performance as the DPD algorithm when the signal parameters are known. When the SNR is lower than the threshold point, the computational complexity can be saved more than90%. Simulation results show that the proposed method is obviously superior to the traditional localization algorithms. However, the FDPD method is currently only supported in terms of several types of signals and some specific conditions. Further work is underway to consider the boundary conditions and the scenario with multiple emitters transmitting complex signals, including NLFM, BPSK, and QPSK.



---

*Source: 2222247-2022-06-18.xml*