Computationally Efficient Sources Location Method for Nested Array via Massive Virtual Difference Co-Array

In this paper, we derive the discrete Fourier transform (DFT) method for direction of arrival (DOA) estimation by generating the massive virtual difference co-array with the nested array. By contrast with the spatial smoothing (SS) subspace-based methods for nested array, which halve the array aperture, the proposed method can take full advantage of the total array aperture. Since the conventional DFT method is a non-parametric method and is limited by Rayleigh threshold, we perform the phase rotation operation to obtain the fine DOA estimates. Owing to the full utilization of the array aperture and phase rotation operation, the proposed method can achieve better performance than SS subspace-based methods for far-field sources especially with massive virtual difference co-arrays which possess a large number of virtual sensors. Besides, as the fast Fourier transform (FFT) is attractive in practical implementation, the proposed method lowers the computational cost, as compared with the subspace-based methods. Numerical simulation results validate the superiority of the proposed method in both estimation performance and complexity.


Introduction
Direction of arrival (DOA) estimation is a significant subject in array signal processing and has various engineering applications such as radar, navigation and wireless communications [1][2][3][4][5][6][7][8]. In the past decades, various subspace-based algorithms, e.g., multiple signals classification (MUSIC) and estimation of signal parameters via rotational invariance techniques (ESPRIT), have been proposed for uniform linear arrays (ULAs) and have drawn considerable attention owing to their high resolution and performance [9][10][11]. Nevertheless, the inter-element spacing of ULAs is restricted to half wavelength to avoid spatial aliasing, which limits the array aperture and resultantly affects the estimation performance [12,13]. Besides, the number of resolvable sources of these conventional methods is less than the number of sensors. That is, an M-sensor array can only resolve M − 1 sources at most. However, the problem of multiple sources detection is of great interest in many applications [14].
Recently, non-uniform arrays based on the concept of difference co-array have aroused wide attention due to the improvement of DOA estimation performance. In [15], the minimum redundancy array (MRA) is proposed to maximize the number of virtual sensors. Unfortunately, there is no (a) We extract the DFT method for the DOA estimation problem by constructing a massive virtual difference co-array with nested array. Besides, it is applicable to any sparse arrays which can generate a large virtual ULA, e.g., coprime arrays. (b) The proposed method remarkably reduces the computational complexity since it applies FFT to obtain initial estimates and searches for fine estimates over a small refined region. Besides, it can avoid the process of eigenvalue decomposition (EVD) and has no need to know the number of sources in advance, which are all inevitable in the existing SS methods and CS methods.
. Without loss of generality, we assume that M is an even integer. According to [25], we take the optimal two-level nested array into account with  Assume that there are K far-field uncorrelated narrowband signals impinging on the nested array from different angles located at ( 1, 2 , ) k kK   . The received signal of the nested array can be represented as [ where k P is the average power of the k -th signal and MM M   I denotes the identity matrix. In practice, the covariance matrix of the signal with a finite number of snapshots [41] is: Assume that there are K far-field uncorrelated narrowband signals impinging on the nested array from different angles located at θ k (k = 1, 2 · · · , K). The received signal of the nested array can be represented as [33] x(t) = As(t) + n(t) = K k=1 a(θ k )s(t) + n(t) (1) where A = [a(θ 1 ), a(θ 2 ), · · · , a(θ K )] ∈ C M×K is the manifold matrix and a(θ k ) = [1, e j2πl 2 sin θ k /λ , · · · , e j2πl M sin θ k /λ ] T ∈ C M×1 is steering vector, where l m ∈ L s (m = 1, · · · , M) denotes the position of sensors, is the signal vector and n(t) stands for the additive white Gaussian noise with mean zero and variance σ 2 n , which is uncorrelated with signals. Here, we assume the number of sources K is known.
The covariance matrix of the signal can be represented as [25]: where P k is the average power of the k-th signal and I M ∈ C M×M denotes the identity matrix. In practice, the covariance matrix of the signal with a finite number of snapshots [41] is: where l = 1, 2, · · · , L and L stands for the number of snapshots.

Vectorization of the Covariance Matrix
According to [42], the covariance matrix of the signal with a finite number of snapshots can be vectorized as: where p = [P 1 , P 2 , · · · , P K ] T , U = vec{I M }, A e = A * • A and I M ∈ C M×M is an identity matrix. Note that the vector y in (4) can be treated as an extended received signal observing with a much longer virtual ULA. It has been presented in [25] that a (M 2 is the direction matrix for the virtual array, whose (n, k)-th element is e j2πn sin θ k /λ (n = −(N − 1)/2, · · · , (N − 1)/2). p can be seen as the signal vector, and u ∈ R N×1 is a vector with a one at the ((1 + N)/2)-th position and the other elements are all zeros.
In [25], by applying the SS technique to the observation vector y v , an [(N + 1)/2] × [(N + 1)/2] matrix R ss as the SS covariance matrix can be achieved. However, almost half the aperture of the virtual ULA is dissipated. To fully utilize the aperture provided by a nested array, we apply the DFT method [40] to y v .

Coarse Initial Estimation
Define the normalized DFT matrix F ∈ C N×N as: where the (p, q)-th element of F is W pq N = e −j2πpq/N . Then the DFT of the virtual steering vector is: whose q-th element is [40]: According to [40], if the number of sensors in the virtual consecutive ULA is infinite, i.e., N → ∞ , as: where δ(·) is the delta function, there always exists an integer p k = N sin θ k /2 that makes [ã v (θ k )] p k = √ N and the others are all zeros. Subsequently, the estimate of θ k (k = 1, 2, · · · , K) can be obtained fromã v (θ k ). However, in practice, the number of sensors of an array is finite even in a massive multiple-input multiple-output (MIMO) system [43]. Then, the power of the round{N sin θ k /2}-th DFT point will not be concentrated at this point, which will depress the accuracy but still can be taken as initial estimates.
Therefore, the DFT of the observation vector can be obtained by z = Fy v and the p-th element is presented as: Accordingly, z has K peaks and we can search the spectrum of z for K largest values whose indices are denoted as p ini k , and obtain their corresponding initial DOA estimates, denoted as:

Fine Estimation
The resolution and estimation accuracy of the DFT method is limited by the Rayleigh limit. However, we will show that the accuracy of the DFT method can be improved with the use of the total array aperture and phase rotation operation.
To decrease the estimation error caused by a finite number of sensors, a phase rotation operation is conducted. The phase rotation matrix is defined as [40]: where η ∈ (−π/N, π/N) is related to the phase shifter. The steering vector corrected by rotation matrix can be achieved as: where the p-th element can be denoted as: Similarly, there must exists a corresponding phase shifter between (−π/N, π/N), referred to as the optimal phase shifter η k , which makes: In this case, the rotated steering vectorã ro v (θ k ) would have only one non-zero element. Then the fine estimate of θ k (k = 1, 2, · · · , K) can be obtained as: We can construct the cost function: where f p ini k denotes the p ini k -th column of the DFT matrix F. By searching η over a small sector (−π/N, π/N) and the optimal phase shifter η k can be obtained by solving the K peaks of z ro = FΦy v . Owing to the phase shifter operation, the power of DFT spectrum is centralized around the theoretical DOA of the signal and hence the performance of the DOA estimation can be greatly improved.

Remark 1.
In practice, we usually use FFT rather than DFT to accelerate the calculation and in this case, we need to employ the zero-padding method [36] if the number of sensors of the virtual ULA is not an integer power of 2.
Remark 2. According to the above derivation of the proposed algorithm, it has no need to know the number of sources in advance, which is attractive in practical implementation.

Remark 3.
The DFT based method is also suitable for other arrays that can be transformed to a virtual ULA, such as the coprime array [16].

Detailed Steps
The detailed steps of the proposed method are as follows: Step 1 Compute the covariance matrix and get the observation vector y v in (5); Step 2 Compute the DFT of y v and search the spectrum of z for K values to obtain θ ini k (k = 1, 2, · · · , K); Step 3 Compute z ro and search η over (−π/N, π/N) to obtain the optimal phase shifter η k (k = 1, 2, · · · , K); Step 4 Compute fine DOA estimates according to (16).

Remark 4.
Our method is designed for far-field sources coming from different directions. Our method fails to resolve sources with the same directions but different ranges. Such sources can be resolved in the near-field cases [44][45][46].

Computational Complexity
In this subsection, we mainly discuss the complex multiplication operation. Specifically, computing covariance matrix of the received signal needs the complexity of O(M 2 L). To accelerate computation, we apply FFT with the complexity of O(N log N). Besides, the process of peak locating and shifted phase searching needs O(N) and O(2πK/∆) respectively, where ∆ denotes the search grid. To sum up, the total complexity of the proposed method is O(M 2 L + N log N + N +2πK/∆). It should be noted that small ∆ leads to high accuracy but high computational complexity as well. However, for massive nested arrays which can produce extremely large number of virtual sensors, high accuracy and low complexity can be obtained even with a large search grid ∆.
Specific complexities of the proposed method along with the SS-ESPRIT [25] and the SS-MUSIC [25] are listed in Table 1 for clarity. The comparison of the computational complexity versus number of sensors is illustrated in Figure 2, where K = 3, L = 100 and the search grid ∆ = 0.0001 rad. As the proposed method does not need to perform EVD and searches over a small sector, it shows clearly that its complexity is much lower than the SS methods.

Method Complexity
Proposed

Resolution and DOF
The resolution of the DFT method is limited by the Rayleigh limit, it cannot resolve sources between two adjacent peaks of the spectrum. Two close sources at 1  and 2  can be resolved only

Resolution and DOF
The resolution of the DFT method is limited by the Rayleigh limit, it cannot resolve sources between two adjacent peaks of the spectrum. Two close sources at θ 1 and θ 2 can be resolved only when |sin θ 1 − sin θ 2 | ≥ 4/N holds, where N is the number of sensors of the virtual ULA.
The maximum number of sources that the DFT method can resolve is N/2 . Figure 3 shows the maximum peaks in the positive axis without considering the noise, where the DOAs θ p = sin −1 (4p/N), p = 0, 1, · · · , N/4 correspond to the even indices of the peaks and M = 16, N = 143. We can see that all 36 peaks can be located and there are 35 peaks in the negative axis, all 71 peaks can be located for specific DOAs (correspond to the odd or even integer indices exactly).

Resolution and DOF
The resolution of the DFT method is limited by the Rayleigh limit, it cannot resolve sources between two adjacent peaks of the spectrum. Two close sources at 1  and 2  can be resolved only

•
The proposed method outperforms SS-ESPRIT and SS-MUSIC [25] for distant sources DOA estimation with nested arrays as it involves the full aperture of the virtual ULA and a phase rotation operation is applied for fine estimation. • The proposed method is computationally efficient since it applies DFT (FFT) and searches for peaks over a small sector. Besides, it can avoid performing EVD operation, which is time-consuming and inevitable in SS methods.

•
The proposed algorithm has no need to know the number of sources in advance, which is practical and attractive in the realistic scenario.

Simulation Results
In the simulation section, we only consider the two-level optimal nested array, i.e., two subarrays have the same number of sensors. The root mean square error (RMSE) is used as the performance comparison metric, which is defined as: where Q is the number of Monte Carlo simulations,θ k,q is the estimate of the q-th trial for the k-th theoretical angle θ k . And in this paper, we set Q = 1000.

Example 1.
The DFT based method is a non-parametric method and its resolution is restricted by the Rayleigh limit. Figure 4 shows the RMSE of the proposed method with the increase of the angular separation of two sources, where M = 16, L = 200, SNR = 20 dB, N = 143, and the search grid for fine estimation is ∆ = 0.0001 rad.
As mentioned above, the resolution is determined by the aperture of virtual array as sin −1 (4/N). In this case, the angular resolution is about 1.6 • . It is depicted clearly in Figure 4 that, when the separation is less than 2 • , the proposed algorithm has difficulty to resolve the two sources. But when sources get far, the estimation accuracy improves significantly owing to the phase rotation operation. Considering the large virtual aperture the massive nested arrays can produce, the proposed method can obtain a rather high resolution.


The proposed method outperforms SS-ESPRIT and SS-MUSIC [25] for distant sources DOA estimation with nested arrays as it involves the full aperture of the virtual ULA and a phase rotation operation is applied for fine estimation.  The proposed method is computationally efficient since it applies DFT (FFT) and searches for peaks over a small sector. Besides, it can avoid performing EVD operation, which is timeconsuming and inevitable in SS methods.  The proposed algorithm has no need to know the number of sources in advance, which is practical and attractive in the realistic scenario.

Simulation Results
In the simulation section, we only consider the two-level optimal nested array, i.e., two subarrays have the same number of sensors. The root mean square error (RMSE) is used as the performance comparison metric, which is defined as:  Figure 4 shows Figure 4 that, when the separation is less than 2 , the proposed algorithm has difficulty to resolve the two sources. But when sources get far, the estimation accuracy improves significantly owing to the phase rotation operation. Considering the large virtual aperture the massive nested arrays can produce, the proposed method can obtain a rather high resolution.  [10 ,30 ,50 ]  θ is given in Figure 5, where 16 M  , =200 L , and the search grid is set to be 0.0001 ad r  for both the proposed DFT algorithm and the SS-MUSIC algorithm. Moreover, the Cramer-Rao Bound (CRB) [47,48] of the nested array is also provided in Figure 5 to evaluate   Figure 5, where M = 16, L= 200, and the search grid is set to be ∆ = 0.0001 rad for both the proposed DFT algorithm and the SS-MUSIC algorithm. Moreover, the Cramer-Rao Bound (CRB) [47,48] of the nested array is also provided in Figure 5 to evaluate the theoretical performance of the proposed method. It is shown clearly that the DFT method outperforms the others as it can utilize the whole virtual consecutive array, while the SS methods sacrifice half of the array aperture to perform spatial smoothing technique. the theoretical performance of the proposed method. It is shown clearly that the DFT method outperforms the others as it can utilize the whole virtual consecutive array, while the SS methods sacrifice half of the array aperture to perform spatial smoothing technique.   Figure 7, where the snapshots is 1000 L  and the search grid is 0.0001 ad r 

Example 4. The estimation performance of the proposed algorithm with different number of sensors is illustrated in
. Obviously, the performance of the proposed algorithm improves a lot as the number of sensors increases. Since the massive nested array can produce an extremely large number of virtual sensors, the proposed algorithm can obtain high angular resolution as well as high estimation performance at the same time.  the theoretical performance of the proposed method. It is shown clearly that the DFT method outperforms the others as it can utilize the whole virtual consecutive array, while the SS methods sacrifice half of the array aperture to perform spatial smoothing technique.   Figure 7, where the snapshots is 1000 L  and the search grid is 0.0001 ad r 

Example 4. The estimation performance of the proposed algorithm with different number of sensors is illustrated in
. Obviously, the performance of the proposed algorithm improves a lot as the number of sensors increases. Since the massive nested array can produce an extremely large number of virtual sensors, the proposed algorithm can obtain high angular resolution as well as high estimation performance at the same time.   Figure 7, where the snapshots is L = 1000 and the search grid is ∆ = 0.0001 rad. Obviously, the performance of the proposed algorithm improves a lot as the number of sensors increases. Since the massive nested array can produce an extremely large number of virtual sensors, the proposed algorithm can obtain high angular resolution as well as high estimation performance at the same time.   Figure 8 that owing to the large aperture of its virtual array, the nested array can achieve better estimation performance even with a much smaller number of physical sensors.    Example 5. In this example, we compare the estimation performance of ULA and nested array, where L = 500, ∆ = 0.0001 rad, the total number of physical sensors is M = 120 for ULA, and M = 20 for the nested array. It is shown explicitly in Figure 8 that owing to the large aperture of its virtual array, the nested array can achieve better estimation performance even with a much smaller number of physical sensors.   Figure 8 that owing to the large aperture of its virtual array, the nested array can achieve better estimation performance even with a much smaller number of physical sensors.

Conclusions
In this paper, we derive the DFT method for DOA estimation with a massive virtual difference co-array generated by nested array. By contrast with the SS subspace-based methods for a nested array which halves the array aperture, the proposed method takes full advantage of the total array aperture. Moreover, since the conventional DFT method is a non-parametric method and is limited by Rayleigh threshold, we conduct the phase rotation operation to achieve fine DOA estimates. Specifically, due to the full utilization of the array aperture and phase rotation operation, the proposed method can attain better performance than SS subspace-based methods for far-field sources especially with a massive virtual difference co-array which consists of a large number of virtual sensors. Besides, as the FFT is attractive in practical implementation, the proposed method lowers the computational cost, as compared with the subspace-based methods. Due to the full utilization of the array aperture of the virtual difference co-array, the proposed method outperforms the commonly-used SS methods for far-field sources.

Conclusions
In this paper, we derive the DFT method for DOA estimation with a massive virtual difference co-array generated by nested array. By contrast with the SS subspace-based methods for a nested array which halves the array aperture, the proposed method takes full advantage of the total array aperture. Moreover, since the conventional DFT method is a non-parametric method and is limited by Rayleigh threshold, we conduct the phase rotation operation to achieve fine DOA estimates. Specifically, due to the full utilization of the array aperture and phase rotation operation, the proposed method can attain better performance than SS subspace-based methods for far-field sources especially with a massive virtual difference co-array which consists of a large number of virtual sensors. Besides, as the FFT is attractive in practical implementation, the proposed method lowers the computational cost, as compared with the subspace-based methods. Due to the full utilization of the array aperture of the virtual difference co-array, the proposed method outperforms the commonly-used SS methods for far-field sources.