A New Computing Paradigm for Off-Grid Direction of Arrival Estimation Using Compressive Sensing

Department of Electrical Engineering, International Islamic University, Islamabad 44000, Pakistan School of Engineering & Applied Sciences, ISRA University, Islamabad 44000, Pakistan Future Technology Research Center, National Yunlin University of Science and Technology, 123 University Road, Section 3, Douliu, Yunlin 64002, Taiwan Department of Computer and Electrical Engineering, COMSATS University Islamabad, Attock Campus, Attock, Pakistan Department of Electrical Engineering, Air University Islamabad, 44000, Pakistan


Introduction
The direction of arrival (DOA) has been under investigation for a long time [1][2][3], and with the advancement in wireless communications, it has applications in a variety of fields like radar, acoustic signal processing, medical imaging, and seismology. The goal of DOA is to estimate the location of the closely spaced sources in the presence of a noise. It is common to use antenna arrays with different structure [4,5], and there are many algorithms for estimation of the DOA [6,7]. The performance of these algorithms depends upon the number of the samples and signal-to-noise ratio (SNR). These algorithms can be divided in to three main categories [8], conventional beamforming, subspace techniques, and maximum likelihood technique. The most notable algorithms used for DOA are multiple signal classification (MUSIC) algorithm, estimation of signal parameter via a rotational variant technique (ESPIRIT), and CAPON. These algorithms require the sampling to be on the Nyquist rate, i.e., sampling frequency should be at least twice the highest frequency present in the signal.
Compressive sensing (CS) has gained a lot of attention over the years because of its ability to exploit the concept of sparsity [9,10]. The CS has applications in the fields like radars [11], image reconstruction and restoration [12,13], blind source separation [14,15], beamforming, and source localization [16][17][18]. A signal can be reconstructed only from a small number of linear measurements if it is sparse in certain domain. It means that the information rate of the signal is much smaller than the suggested signal bandwidth. As most of the real-time signals are sparse, therefore, the signal can be reconstructed using a system of equations such that the smaller number of samples is used. It is important that it satisfies the restricted isometric property (RIP). A recursive weighted minimum norm with focal underdetermined system solver (FOCUSS) is used to achieve sparsity in the problem of source localization [19]. In [20], a hardware has been developed for spectral estimation using CS framework.
Normally, in the CS, we are required to find a sparse signal using an overcomplete dictionary. This concept can also be applied to source localization problem where spatial sparsity is exploited; the source is sparse in the spatial domain means that if it is not present at every angle, hence, the concept of CS can be applied due to the sparsity of source. The source localization can be done using single or multiple measurements to achieve super resolution. For multiple measurements, the computational complexity increases with the amount of data growing in a system. Different methods are present for solving multiple measurement problems.
One of the main issues with the source localization using CS is to define the resolution of the grid [21]. It is assumed that the location of the sources falls on the resolution of the grid defined. However, this is not always possible or practical. There may be a scenario when the location of the source does not coincide with the resolution of the grid defined. This creates a grid mismatch or off-grid targets. There are different methods available to solve this problem; one of them is iterative grid refinement. In iterative grid refinement, the grid resolution is changed until the dense grid mitigates the grid mismatch. One of the main drawbacks of using iterative grid refinement is that decreasing the grid resolution may not comply with the RIP condition. The other approach to resolve grid mismatch problem is off-grid sparse method [22]. In this process, initial grid resolution is still defined. The DOA of the sources is not restricted to the grid. A bias is added to the signal model using first order approximation of the manifold matrix. The new model may be nonconvex and difficult to solve. Iterative grid refinement is one of the most used methods. However, it requires a method to best select the discretization value for the grid. Therefore, in this paper, a method based on the distribution of the source energy due to grid mismatch is presented to calculate the discretization value for the grid. The main contributions of the paper are summarized as follows.
(i) A novel framework for solving grid mismatch or offgrid target is presented for DOA estimation problem using CS technique (ii) A fitness function is introduced based on the difference of source energy between adjacent grids due to grid mismatch in DOA estimation (iii) An approach for finding the best discretization value using the designed objective function is presented for iterative grid refinement (iv) The proposed scheme is viably tested for multiple sources with different energy and spatial resolutionbased scenarios in DOA estimation The rest of the article is organized as follows. The mathematical background for DOA estimation using MUSIC algorithm is presented in Section 2, while in Section 3, compressive framework for DOA estimation is presented. In Section 4, the proposed algorithm is provided, and results along with necessary discussion are presented in Section 5. The concluding remarks are given in Section 6.

DOA Estimation
In this section, a general model for direction of arrival (DOA) estimation based on a subspace technique for DOA estimation is presented. Let us consider a uniform linear array (ULA) as shown in Figure 1 with N number of antenna elements, and the distance between the antenna element is λ/2, while there are P number of far field sources at different angles, θ l . Then, the received signal at the m th antenna element is given as where s i is the amplitude of the signal that is received. k and d are wave number and distance between the antenna elements. Then, Equation (1) can be written as In the presence of noise, the received signal is updated as where n is the Gaussian noise, A is the steering vector, and y is the received vector. As mentioned earlier, there are number of ways to solve Equation (3) like MUSIC and MVDR. If it satisfies the Nyquist sampling rate, the correlation of the received matrix is given as Then, accordingly, the spatial spectrum is given as where e n is the eigen vector orthogonal to the steering vector. The performance of the MUSIC algorithm degrades with the reduction in the number of the samples and presence of a noise. In the next section, we will look at the CS approach to the DOA estimation.

Compressive Sensing
The concept and mathematical development of compressive sensing for DOA estimation are briefly presented in this section. In CS framework, a sparse representation of the signal can be reconstructed only from a small number of samples. Let us consider a signal s, a discrete signal which is sparse in certain domain sε ℂ N×1 and y is the received signal of dimension M such that yεℂ M×1 . Then, received signal without noise is given as Here, A is a sensing matrix of dimension Aεℂ M×N , where M ≪ N. We count the number of nonzeros elements in a signal s, which is given by ksk 0 also known as l 0 -norm. This leads to a nonlinear programming (NP) hard problem. To solve NP problem, many approximation methods have been developed. One of the methods is to use l 1 or l p relaxation. If the unknown signal s is considered sparse, then the optimization problem can be mathematically represented as considering p = 0, then the above equation will be an NP hard problem. While considering p = 1, we can recast it as an l 1 -norm problem and solve it using Equation (9).
In practical scenario, there is always a noise. Now, if the received signal is contaminated with a noise n, then, (7) is rewritten as and the optimization problem becomes where ε is a parameter that specifies how much noise is allowed. To formulate the problem in the CS framework, consider Figure 2 [23], where the spacing between the antenna elements is d which is λ/2. As the goal is to find the location of the sources, we consider a ULA with narrow band signal for K number of sources and M number of antenna elements. The received signal is given in (1). As seen in Figure 2, to cast this in the sparse representation problem, an overcomplete dictionary of array steering vector A is introduced, where A = ½θ 1 , θ 2 , ⋯, θ N , N is the sampling of the grid. The N will be much higher than K. Therefore, the matrix A is given as One of the main issues with application of CS in the DOA problem is definition of the grid resolution. The resolution depends upon the sampling grid formulation, and the sampling grid is uniform. If the grid size is defined very fine, it increases the computational requirements. If the size of the sampling grid is large, then the resolution decreases, and close targets cannot be detected.

Proposed Methodology
In this section, we propose a methodology for the selection of discretization value for the grid. Considering Figure 3, the upper grid represents the resolution of the grid with discretization value of r = θ n+1 − θ n , while θ L represents the location of the source.
The first step is to estimate the vector s using the overcomplete dictionary defined for the iteration

Wireless Communications and Mobile Computing
Due to grid mismatch, the energy of the source is distributed among the adjacent grids as shown in Figure 4. The discretized grid resolution is 1°defined, and the location of the source is 60.4°. The energy of the source distributed among the adjacent grids is mathematically presented as

Wireless Communications and Mobile Computing
Let i be the iteration index in which the discretization value is r i . The peaks in the vector s are detected, and the difference is taken as in Equation (14).
Then, the process is repeated with i + 1 th iteration with finer grid discretized value and with dictionary defined with new discretization value. The fitness function is calculated as A termination criterion is defined asset as If it satisfies the termination criteria, then the discretization value in the i th iteration is the best value for the grid to discretize at. For the case of multiple sources, Equation (14) can be generalized for sum of the difference of adjacent peaks and sum of individual peaks if there is no adjacent peak. As mentioned, i is the iteration number. In each iteration, the discretization value is reduced. It is selected by the user. In our simulations, we have selected a discretization value of 1, 0.5, 0.1, and 0.01. The main steps involved in the proposed algorithm are given as follows.
The regularization term λ in (13) plays important role for the accuracy of the solution and must be estimated. It is a compromise between finding a solution that is sparse as possible and has lower error as possible. Two methods for estimating the regularization term are L curve and generalized cross validation (GCV). In GCV, it is more convenient as compared to L curve where we must find the corner [24]. It can be computed using the following relations.
The GCV estimate is variant of the above equation which is obtained by applying necessary calculation and results in GCV function [25]. This technique estimates λ by assuming that the optimum value of λ should be chosen to minimize GCV value.

Simulation Results
We start our simulations with MUSIC algorithm. In Figure 5, the results of the two sources based on DOA estimation are presented that are at an angle of 20°and°. The received signal SNR is 20 dB, and the number of antenna elements is 10. We consider the performance of the algorithm for different number of samples. The number of samples is 50 and 3. It is observed that as the number of the samples decreases, the detection performance degrades, and the targets are not distinguishable.
Next, we consider CS techniques for solving the DOA problem. Considering a single sample at T = 1, we create an overcomplete dictionary having a resolution of 1°. The location of the targets can be resolved by solving Equation (9) with the help of linear programming. We use convex optimization toolbox for solving this problem. For simplicity, we consider a noiseless case. It is shown in Figure 6 that two targets with amplitude of 2 and 1 on normal scale at location 20°a nd 23°are resolved.
Next, we consider a scenario in which the target locations are not aligned with the grid resolution. Considering two targets, one is at 40.5°with amplitude 2, and the other is at 43.5°w ith amplitude 1. In Figure 7, it is shown that the four targets are detected. This creates ambiguity about the location of the targets and the number of the targets. However, it is observed  that the amplitude of the received signal is distributed in the adjacent grid. To solve this, the proposed algorithm is used. The fitness function is calculated for the pair of adjacent peaks. Then, the grid resolution is changed and again solved with the new dictionary elements, and the fitness function is calculated. The new fitness function calculated is compared with the previous one. If the fitness function increases, the steps are repeated; otherwise, we conclude that the optimized value has been reached, and the true location of the targets have been estimated as shown in Table 1.
Next, we consider a scenario for three and four numbers of sources in the DOA estimation problem. Initially, the grid resolution is taken to be 1°. It is assumed that the targets are not aligned with the grid resolution which again creates ambiguity about the number of the sources as shown in Figures 8 and 9. Table 1 shows the iterations for solving the grid mismatch problem using the proposed algorithm.
Next, we consider a different case, where finer resolution is required to detect the sources. Two sources are considered with locations at 40.4°and 43.3°and amplitude of 2 and 1 shown in Table 2. Similarly, four off-grid sources are considered in Table 3, where we present the results for a scenario of four sources at locations 30.4°, 43.3°, 60.5°, and 80.7°. The results show that the location of the sources is accurately estimated using the proposed algorithm.
As shown in the tables with the help of the fitness function, best discretization value for the grid is calculated. As mentioned, to address the basis mismatch, off-grid sparse method is also used. The signal model for the off-grid target with bias is given in [26,27]. The minimization problem can be written as where B = ½bðθ 1 Þ, ⋯, bðθ N Þ and bðθ nk Þ is the derivative of a ðθ k Þ with respect to θ nk . Δ = diag ðδÞ = ½δ 1 , ⋯, δ N T . δ is the difference between the nearest grid point and the direction of the k th signal.
In Figure 10, we consider two targets located at 30:5°and 60:5°. The regularization parameter for each SNR level is calculated using the GCV method. The mean square error (MSE) of the DOA estimation of the proposed method is compared with Cramer Rao lower bound (CRLB) for DOA estimation and off-grid method with bias. For a greater SNR-based scenario, the proposed algorithm is reasonable accurate. Additionally, the proposed approach is simpler and requires less computations.
In Figure 11, we compare the iterations for grid refinement with MSE. We consider two targets at 30.5°and 60.5°, where the grid is refined in each iteration. The resolution of the grid for each iteration is 1°, 0.5°, and 0.1°. It is observed that the MSE is the lowest for 0.5°grid resolution. Although the targets are detected for grid resolution of 0.1°, the MSE increases due to the RIP condition.
In Figure 12, we show the overall process flow structured diagram of the proposed study.

Conclusion
Compressive sensing is an interesting technique for finding the location of the sources using sensor array. However, the definition of the grid is a challenge. If the location of the target coincides with the grid resolution, it is detected, whereas if not, then the signal power of the source is distributed among the adjacent grids. This creates ambiguity about the location of the target. In this paper, we presented an iterative grid refinement mechanism based on a fitness function. The fitness function governs the extent of the grid refinement. When the maximum value of the fitness function is reached, we conclude further refinement of the grid is not required, and the ambiguity about the location of the target is removed. Thus, the best discretization value for the grid is calculated. In the future work, this technique can be applied for multiple input multiple output radar system [28].  Grid refinement based on fitness function 1. s i = min ||y -A i s|| 2 + d|s| 1

Data Availability
There is no data associated with this manuscript.

Conflicts of Interest
The authors declare that they have no conflicts of interest.