DOA estimation of LFM signal based on Krylov subspace method

Direction of arrival estimation of LFM signal is an essential task in radar, sonar, acoustics and biomedical. In this paper, a short time Fourier transform multi-step knowledge aided iterative generalized minimum residual (STFT-MS-KAI-GMRES) approach is presented to amend the angle measurement of this signal. A three stage algorithm is proposed. First, the process is initiated with formulating an estimation algorithm for the carrier frequency and chirp rate, followed by calculation of STFT of the output of array element; this yields a spatial time–frequency distribution matrix. Next, the Krylov subspace-based estimation algorithm is formulated in the presence of MS-KAI-ESPRIT algorithm. If the number of antennas increases, the accuracy of the algorithm will increase, but we will incur more communication costs. Results are presented showing attainment of the CRLB by STFT-MS-KAI-GMRES the for an adequately large signal to noise ratio. An important feature of the method presented in the current study is the low computational complexity that has higher suitability for practical applications.


Introduction
One of the significant matters in array signal processing is estimation of direction of arrival (DOA), attracting considerable attention [1,2]. Some characteristics of the Linear Frequency Modulation (LFM) signal as a chirp signal include having a better detection performance and transferring a high amount of information and. Because of these characteristics, it has a wide application in radar, solar, and other detection devices [3][4][5]. Therefore, one of the focus points of the research works is the DOA estimation of LFM signals [6][7][8]. There exist ideal energy collection properties in the time-frequency plane in LFM signals so the time-frequency analysis systems are superior in addressing the LFM signals [9]. The spatial time-frequency distribution (STFD) concept was firstly developed by researchers in 1999, and they used the bilinear time-frequency distribution in the area of the DOA estimation [10]. Besides, some DOA estimation algorithms were developed based on the time-frequency systems [11][12][13][14]. In [11], DOA estimation was proceeded using Wigner-Ville distribution (WVD), and a satisfactory estimation performance was obtained. Despite owning the B Hamidreza Bakhshi bakhshi@shahed.ac.ir 1 Department of Electrical Engineering, Shahed University, Tehran, Iran optimal time-frequency resolution by WVD, it has a fairly high computational complexity, and its cross-terms have a serious impact on the accuracy of the estimation concerning multiple signals. Development of such algorithms is further limited by these factors. DOA estimation of LFM signals was studied in [12,20] using fractional Fourier transform (FrFT). A new array data model was presented in this work in the fractional Fourier transform domain (FrFD). For obtaining the DOA estimation, they employed the traditional Eigen subspace approaches. These researchers used FrFT instead of WVD to solve the cross-term interference problem. Nevertheless, when the new steering vector matrix is constructed, the FrFT-based algorithms are required to search the energy collection point in two dimension that creates a high level of computational complexity. In [13] and [14], a DOA estimation algorithm is built based on short-time Fourier transform (STFT) that is effective in reducing cost of algorithm implementation. Moreover, an effective approach is proposed in these works for precise selection single-source auto-term location for construction of the STFD matrixes of each source, and a satisfactory estimation performance was achieved. The STFT-MS-KAI-GMRES is used in this work, which has a low computational complexity and high resolution in low SNR, in comparison with other methods, like MS-KAI-CG [15,19]. An ESPRIT-based method called MS-KAI-ESPRIT was developed in previous studies [16], which was presented for improving the sample covariance matrix estimation in the finite sample size region for stationary signals. It is possible to apply the STFT-MS-KAI-GMRES algorithm for estimating the DOA of multiple LFM signals. Cross-term interference is not generated in this algorithm, and its estimation precision is higher in comparison with the traditional methods. Additionally, it preserves the computational advantage of the traditional ones.
The present study proposes a new method to estimate DOA and the chirp parameters. We develop the STFT-MS-KAI-GMRES algorithm with the ability of obtaining considerable gains in RMSE or probability of resolution performance compared to other methods, and it is a good option for usages with short data records in large-scale antenna systems for wireless communications, radar, and other large sensor arrays.
The major advantage of this method is the computational load reduction. Thus, it can be realized usefully for real time usages. Besides, it is possible to use the presented chirp rate estimation mechanism in low SNR scenarios.
This study is presented as following: the signal model is explained in the second section, and an estimation algorithm is formulated for the carrier frequency and chirp rate. The third section presents the estimation algorithm for the carrier frequency and chirp rate. Then, the STFT of the output of the element is obtained. By selection of the time-frequency points in the energy collection area of LFM signals, the spatial time-frequency matrix can be gained. Lastly, the Krylov subspace-based estimation algorithm (with effective extraction of prior knowledge of LFM signal) is developed in the presence of MS-KAI-ESPRIT algorithm. MLS criterion is used for obtaining the closed-form solution. According to the simulation results in IV, the present algorithm provides a satisfactory performance. The conclusion and advantage of the method proposed in this paper are discussed in the last section.

System model
Consider P far-field LFM signal radiating from source location θ [θ 1 , θ 2 , . . . , θ P ] T . The LFM signal can be modeled as: where (.) T denotes transpose, b is magnitude of signal, k represents chirp rate, and f c indicates the carrier frequency. It is assumed that s(t) is normalized to its magnitude. Let s 1 (t), which is the sum of s(t) and s(t − τ ) (its timedelayed form): where τ is an arbitrary delay parameter. The function B(t) has roots as Eq. (3): so we have Eq. (4): by using Eq. (4), the chirp rate is estimated as follows: For calculating the carrier frequency estimation (like the chirp rate estimation process mentioned above but with a slight difference), by aggregation s(t) and base-band signal and calculating the roots of the new function, the estimation of carrier frequency will be obtained: Also, the space between the sensors is a distance (d ≤ λ c 2 ), where λ c represents the signal wavelength, and without generality loss, there is −π 2 ≤ θ 1 ≤ · · · ≤ θ P ≤ π 2 . Modelling of the i th data snapshot of the M-dimensional array is as follows:

x(i) A(θ)s(i) + n(i)
where is data vector of LFM signals that chirp rate and carrier frequency estimated in Eqs. (5), (6), and n(i) ∈ C M×1 is vector of noises which are additive white Gaussian noise with zero mean and varianc σ 2 n and C is a complex set of numbers. The Vandermonde matrix A(θ) [a(θ 1 ) . . . a(θ P )] ∈ C M×P called the array steering, possesses the array steering vectors a(θ j ) equivalent to the nth source that is given as follows: where n 1 . . . P.
Given modelling n(i) and s(i) uncorrelated linearly independent variables, calculation of M × M signal covariance matrix is as follows [21]: Due to unknown state of the true signal covariance matrix, the following formula is used for its estimation: which estimation accuracy depends on N (number of snapshots).

The STFT-based spatial time-frequency distribution
One of the oldest and most well-known time-frequency analysis approaches is STFT. It was developed in 1946 by Gabor [17]. Its underlying concept is cutting out one sequence of time domain signal by the use of window function and proceeding Fourier transform for obtaining the signal's local time-frequency properties. STFT is defined as follows: where h(t) represents the window function. STFT is firstly applied on the nth element output, with obtaining: and S x(n) (t, f ) is presented as follows: which S X (t, f ), S S (t, f ) and F(t, f ) respectively are STFT for x n (t), signal and noise. We can get the STFD based on STFT via the array output's correlation operation, so we have: where D xx (t, f ) is STFD matrix of the array output and A (θ) is array steering matrix. The mathematical expectation ("Appendix") of STFD matrix can be calculated by:

STFT-based MS-KAI-GMRES algorithm
MS-KAI-GMRES algorithm is regarded as Krylov subspace algorithms. The critical point in the suggested algorithm is modification of the sample data STFD matrix estimation at each iteration, which is achieved through gradual incorporation of the knowledge presented by the newer Vandermonde matrixes that contain the newer estimates from the previous iteration. We can calculate refined estimates of the projection matrixes of the noise and signal subspaces according to these updated Vandermonde matrixes. A a scaled version of the undesirable terms is derived fromD for computing modified STFD matrix. Based on [14], Eq. (16) is obtained by generalizing Eq. (10):

( As(i) + n(i))( As(i) + n(i)) H
The signal and the noise elements are represented by the first two terms ofD in Eq. (16). A means steering matrix, N denotes number of snapshots, and n and s represent noise and signal vector. In Eq. (16), the last two terms are undesirable by-products that can be considered as estimates for the correlation between the noise and signal vectors. The studied system model is based on zero-mean noise vectors and independent of the signal vectors. Hence, the noise and signal components are not correlated. Consequently, for an adequately large number of samples N, the last two terms of Eq. (16) are zero. However, in practice (e.g., in radar applications) there are a limited number of available samples. In these cases, the last two terms in Eq. (16) could take significant values, creating the deviation of the estimates of the noise and signal subspaces from the true noise and signal subspaces.

end if end for end for
Then, the GMRES algorithm is used for estimating the DOAs [18]. The major point in the presented algorithm is that STFD matrix is modified at each iteration. Table 1 presents the algorithm's steps. The superscript (.) (1) denotes the estimation task that is performed in the first step. A process including n 1 : P iterations is initiated using the DOA estimates by formation of the Vandermonde matrix. Then, it is followed by estimation of the amplitudes of the sources and the square norm of the differences between the vector containing estimates and observation vector and the available DOAs is reduced [21]. It should be noted that the calculation of carrier frequency and chirp rate is done regarding Eqs. (5) and (6). Formulation of this problem is as follows: Using the least squares method, Eq. (17) is minimized, and the solution is presented as: Then, we estimate the noise component as the difference between the observations made by the array and the estimated signal, as follows: When the noise and signal vectors are estimated, we can compute the third term in Eq. (16): (20) whereÂ,ŝ,n andQ A respectively mean: estimation of steering vector, estimation of signal vector, estimation of noise vector and estimation projection matrix of the signal subspace: is an estimate of the signal subspace's projection matrix, and: whereQ ⊥ A of is an estimation of the noise subspace's projection matrix. Then, we calculate the modified data STFD matrixD (n+1) as part of the process of n 1 : P iterations, which is gained by computation of a scaled version of the estimated terms from the primary sample STFD matrix as follows: where the superscript (.) n represents the nth iteration. The reliability or scaling factor μ is incrementally increased from 0 to 1, which leads to modified data STFD matrixes. They give origin to newly estimated DOAs That is indicated also by the superscript (.) n+1 . Here, the GMRES algorithm is used that was already explained. As a result of the Eigen value decomposition (EVD) of the modified STFD matrix Eq. (16): where the matrixÛ n ∈ C M×M−P indicates the noise subspace and the matrixÛ s ∈ C M×P indicates the signal subspace. The P largest and the M − P smallest Eigen values are observed in the diagonal matrixesˆ s andˆ n . Then, using steering vectors of newer estimated DOAs, a new Vandermonde matrix is built. This new matrix is used for computing the newer estimates of the projection matrixes of the noiseQ ⊥ (B n +1) and signalQ (B n +1) subspaces. Then, the newer estimates of the projection matrixes are employed for computing the initial sample matrixD, and the number of sources and sensors, the stochastic maximum likelihood objective function U (n+1) (μ) for values of μ at the nth iteration, as in the following: ln(.) and det(.) are described in Table 1 (part B).
The unavailable DOA estimates with a higher likelihood at each iteration are selected by the preceding computation. Next, it is possible to determine the set of estimated DOAs equivalent to the optimal value of µ, minimizing Eq. (25) also at each nth iteration. In the end, the set of the estimates gained at the Pth iteration forms the presented MS-KAI-GMRES algorithm output (Table 1).

Experiment setups and simulator details
In this paper, all simulations are performed by a system with windows 10, RAM 8GB, Core i7 and a 64-bit processor. This section uses Matlab 2018b software to simulate some subspace-based methods (MUSIC, ESPRIT, two-step ESPRIT and multi-step ESPRIT).
To evaluate the performance of each method, the criterion of mean squares error is used.   Fig. 2. For the purpose of having the same basis of comparison for the RMSE and the resolution probability of all mentioned algorithms, it is assumed that there is no previous knowledge, i.e., despite availability of known DOAs, MS-KAI-GMRES is computed as they are not available.

Simulation results
In this paper, it considered that there is a ULA with M 10 sensors, P 2, inter-element spacing λ C 2 and assuming the number of snapshots as N 1000. It is supposed there is not a prior knowledge of the last true DOAs (10.2°, 12.6°, 15°and 17.4°). The proposed MS-GMRES algorithm's performance is higher than GMRES, given [18] in the range between − 5 and 20 dB. Figure 2 indicates the RMSE in dB versus SNR. Here, the term CRLB represents the square root of the deterministic Cramer-Rao bound. The definition of the RMSE is as follows: where L is the trial number. According to the results, MS-KAI-GMRES outperforms in low SNR.
In Table 2, the effect of increasing the number of antennas can be seen.

Conclusion
The proposed STFT-MS-KAI-GMRES algorithm, gradually uses the knowledge of source signals that is achieved online. It also exploits the STFD matrix structure and its perturbations. Firstly, it initiates with estimating the parameters of the signal's carrier frequency and chirp rate by the method  f summing the main signal with delayed format and signal summing with the base-band format, respectively. Then, the nonstationary signal (LFM) is taken to the time-frequency domain, and by solving the problem of low SNR and few snapshots, and applying the GMRES method to estimate the STFD matrix, the DOA of LFM signal is estimated. As mentioned, the GMRES method is one of Krylov's subspace methods, and the computational load of the proposed method is good compared to its high accuracy. The excellent performance can justify the pretty high computational cost that is linked with additional matrix multiplications, the iterations required for progressive incorporation of the knowledge achieved online as newer estimates, and the increment used for reducing the unwanted by-products. The low computational complexity is the main characteristic of the method proposed in this paper that helps its feasibility for real-time applications. Additionally, it is possible to employ the mechanism presented for chirp rate estimation under low SNR. In the case of a moving target with different speeds, the scenario of estimating the chirp rate and carrier frequency as well as the target angle is challenging and can be considered as future research.