Underdetermined DOA Estimation for Wideband Signals via Focused Atomic Norm Minimization

In underwater acoustic signal processing, direction of arrival (DOA) estimation can provide important information for target tracking and localization. To address underdetermined wideband signal processing in underwater passive detection system, this paper proposes a novel underdetermined wideband DOA estimation method equipped with the nested array (NA) using focused atomic norm minimization (ANM), where the signal source number detection is accomplished by information theory criteria. In the proposed DOA estimation method, especially, after vectoring the covariance matrix of each frequency bin, each corresponding obtained vector is focused into the predefined frequency bin by focused matrix. Then, the collected averaged vector is considered as virtual array model, whose steering vector exhibits the Vandermonde structure in terms of the obtained virtual array geometries. Further, the new covariance matrix is recovered based on ANM by semi-definite programming (SDP), which utilizes the information of the Toeplitz structure. Finally, the Root-MUSIC algorithm is applied to estimate the DOAs. Simulation results show that the proposed method outperforms other underdetermined DOA estimation methods based on information theory in term of higher estimation accuracy.


Introduction
Direction of arrival (DOA) estimation is an important issue in numerous fields of array signal processing, which has widespread applications in radar, sonar, acoustics, navigation, communications, speech enhancement [1][2][3] and multiple-input multiple-output (MIMO) systems [4]. Along with the recent development of underwater sonar systems, DOA estimation could supply the direction information, which plays various important roles in underwater target tracking [5][6][7] or target localization. Currently, plenty of mature and robust DOA estimation methods, such as estimating signal parameter via rotational invariance techniques (ESPRIT), multiple signal classification (MUSIC), etc. have achieved super-resolution estimation performance for narrowband signals [8]. Further, compared with narrowband signal sources, the wideband counterparts are able to provide more useful and interfering information where the array phase differences lie on two factors, i.e., the DOAs of sources and frequencies. Meanwhile, in passive sonar system, the radiation noise signal [9] from underwater vehicles is also a wideband signal. In light of such characteristics, wideband DOA (WDOA) estimation has attracted more research interests [10][11][12]. Among the classical wideband processing algorithms, the incoherent signal subspace method (ISM) can only deal with uncorrelated signal, norm of a. Symbols · F , Rank(·), and tr(·) stand for the Frobenius norm, rank, and trace of a matrix, respectively. Expectation operator is denoted by E [·]. diag(·) implies that a diagonal matrix is formed from the given vector as the diagonal elements. vec{·} denotes the vectorization operator. ⊗ denotes the Kronecker product and refers Khatri-Rao product. · represents rounding down.

Signal Model and Problem Formulation
Consider a two-level nested arrays of M sensors with the M 1 = M 2 sensors in the inner sub-array, while M 2 = M − M 1 sensors in the outer sub-array, and both sub-arrays are the uniform linear array (ULA). d 1 denotes the separation between the sensors of the inner sub-array, while d 2 represents the counterpart of the outer sub-array, which is equivalent to (M 1 + 1)d 1 . The NA structure is shown in Figure 1, in which the position of the sensors are given by the sets (1) Assume that K wideband stationary signals at the directions of θ 1 , θ 2 , ..., θ K from the far-field impinge on the nested array, where θ k represents the kth signal source. Then, the received signal x m (t) of the mth sensor can be expressed as where the w m (t) is noise part, whose entries are zero-mean Gaussian random variables with σ 2 variance. τ m (θ k ) represents the propagation delay associated with the kth signal impinging and the mth sensor. Wideband signal after time-sampling is split into L non-overlapping segments, and an N-point discrete Fourier transform (DFT) is carried out over each segment data. The output at frequency bin f n is modeled as x n,l = A n (θ)s n,l + w n,l n = 1, 2, · · · , N; l = 1, 2, · · · , L, where s n,l is the signal vector and w n,l represents the additive noise vector in the frequency domain. The array manifold matrix A n (θ) ∈ C M×K of the nth frequency is A n (θ) = a n (θ 1 ), a n (θ 2 ), · · · , a n (θ K ) , whose columns denote array steering vectors for K signal sources of the nth frequency bin, that is where a n (θ k ) ∈ C M , τ m (θ k ) = h m sin(θ k )/c, and c is the signal source propagation velocity.
Then, the covariance matrix R x n of the nth bin data can be expressed as where R s n = E s n,l s H n,l = diag(ρ n ) denotes the only signal covariance matrix at the frequency bin f n , and ρ n = ρ 2 n,1 , ρ 2 n,2 , · · · , ρ 2 n,K T is the source power vector, while R w n = E w n,l w H n,l = σ 2 w I is the noise covariance matrix, and σ 2 w = σ 2 n,1 , σ 2 n,2 , · · · , σ 2 n,M T represents the noise power vector.
In practice, the covariance matrix R x n is estimated from the sample covariance matrixR x n with the L non-overlapping segments, andR x n is given bŷ VectorizingR x n in Equation (7), we have where with b n (θ k ) = a * n (θ k ) ⊗ a n (θ k ), where q n = ρ 2 n,1 , ρ 2 n,2 , · · · , ρ 2 n,K T represents the signal vector, andǏ = vec{I} ∈ C M 2 ×1 . Note that B n (θ) is a tall matrix and can provide useful information for increasing the DOF; in other words, M < Rank B n (θ) < M 2 . According to the concept of difference co-array [19], we can get a new virtual array position of the ULA The goal of the paper is to determine the θ from r n in Equation (8).

Virtual Array Signal
Recall that the vector r n in Equation (8) contains the repeat rows, which do not help for improving the DOF. After removing the repeat rows and sorting its remaining rows of the r n , we obtain a new vector by the following operation z n = G n (θ)q n +w n , where z n ∈ CM andw n = σ 2 wĨM ∈ CM represent the new observation vector and the new noise vector, respectively. G n (θ) = g n (θ 1 ), g n (θ 2 ), · · · , g n ∈ CM ×K is redefined as the filled virtual manifold matrix, which is equivalent to the counterpart of a ULA withM sensors. The virtual array steering vector g n (θ k ) exhibits a Vandermonde structure of sizeM as follows

Focusing on the Virtual Array
It is known that the CSM, which is a classic way of wideband signal processing [31], can lead to less complexity and better DOA estimation performance. The key point of CSM is to design focusing matrices {P n } N n=1 , whose purpose are to transform the manifold matrices at frequencies {f n } N n=1 to that at the pre-selected reference frequency f 0 as follows [32] where G 0 (θ) denotes the virtual array manifold matrix at the referenced frequency f 0 . Note that rotational signal subspace (RSS) method is a well known way of CSM [13], whose focusing matrices {P n } N n=1 are gained through minimizing a Frobenius norm of the obtained virtual array manifold errors where P n should be a unitary matrix, yielding [33] We now define an excessive matrix E = G n (θ)G 0 (θ) H , whose left and right singular vectors are defined as U l n and U r n , respectively. Accordingly, the product of the virtual array y n is defined through multiplying the focusing matrix byz n y n = P nzn = P n G n (θ)q n + P nw After focusing processing, we can reconstruct a new single wideband signal model by averaging the obtained signal of each frequency as follows

Source Number Detection Using Information Theory Criteria
As is known, the key factor for source number detection is the obtained sample covariance matrix. Unfortunately, we cannot get an optimal covariance matrix directly by the new reconstructed single wideband signal in Equation (18). The spatial smoothing method is one of effective solutions for such problem. We can divide the virtual array intoM +1 2 + 1 subarrays, each withM −1 2 − 1 sensors. The output of the n th subarray is defined asȳ n , whose corresponding covariance matrix R n is written as Meanwhile, we can obtain the averaged covariance matrix R, that is Afterwards, applying EVD on the matrix R, we have where U and λ are the eigen vector matrix and eigenvalues, respectively. Here, we have the non-increasing order of such eigenvalues We next introduce the information theory methods: AIC, AIC-C3, and MDL [29,30] for source number detection. These three methods can be organized as follows: where Λ(k) is the maximum likelihood and whereL and k represent the number of snapshots and the number of DOF, respectively. Then, we estimate the signal source number K bŷ

ANM Principle
In the literature [28,34], the ANM principle is that the signal is expressed as a linear combination of a few atoms over a known atom set, and then, the structural information of these atoms, i.e., Vandermonde structure, is utilized for signal reconstruction from measurements.
We assume the signal Φ contains several components, which exhibit the same Vandermonde structure and especially belong to a basic atom set A. Therefore, we can define the atomic norm of Φ over the atom set A, which has the form whose main task is to seek the sparsest decomposition of Φ over A. In other words, if the signal Φ has a special form, i.e., a linear combination form, the signal Φ is considered to be sparse over the atom set of A in Equation (29). Knowing the A, the components of Φ can be retrieved using 1 norm minimization as follows: arg min Equation (30) can be tackled by seeking out the atomic norm Φ A , while the optimal atom set is also obtained for estimating DOAs.
In the presence of noise, the observed measurement Ψ = Φ + w contains the signal and the noise, thus the corresponding atomic norm is written as where ε denotes the noise threshold. The problem in Equation (31) is called the atomic norm minimization, which can be solved via SDP.

DOA Estimation
Here, recall Equation (18); in the case of being noiseless, Equation (18) only contains the signals information, which can be rewritten as , and all K steering vectors {g 0 (θ k )} K k=1 have the same structure, i.e., Vandermonde structure. Furthermore, they belong to a known atom set A v . According to the ANM definition, the signal s can be considered as the sparse signal over the atom set A v .
Afterwards, the set of atoms by all possible G 0 (θ) can be readily formed as follows: Apparently, the noiseless version of Equation (32) ofȳ is viewed as a positive combination of complex sinusoids over the atom set A v . It has been indicated that the atomic norm of s is said to be minimized, when only those atoms corresponding to the true signals {θ k } K k=1 are selected to linearly describe s , which has the form To practically solve Equation (34), an SDP formulation of s A is given as where T(u) is a Toeplitz matrix. According to the authors of [28,35], any PSD Toeplitz matrix allows for Vandermonde decomposition. For any {q k , θ k } K k=1 , K <M , define u = ∑ K n=1 |q k |g 0 (θ k ) and the Toeplitz matrix T(u) is formed from the first row of the vector u = [u 1 , u 2 , · · · , uM], that is The SDP yields the optimal vector u and the Toeplitz matrix T(u), in which the true signals θ are estimated.
In practice, we usually do not have the s directly, and the measurementȳ contains the noise in Equation (18). According to Equation (31), the ANM formulation can be written as Equation (37) makes a regularized denoising formulation to recover s . ŝ A v represents the sparsity enforcing term, and the quadratic term of Equation (37) is used to control the noise. ξ is a regularization parameter.
According to Equation (35), Equation (37) now can be converted to an equivalent SDP formulation as follows: where ξ is a regularization parameter balancing Toeplitz structure and the noise tolerance to the observationȳ, i.e., ȳ −ŝ 2 2 (least-squares term), whose poor choice may increases the probability of a bad solution. Therefore, an optimal section of ξ plays a very important role in guaranteeing signal recover. The SDP formulation in Equations (35) and (38) can be implemented efficiently by the CVX toolboxes.
Remark on targets angle separation condition: Here, targets angle separation means the required angle separation between neighbouring signal source. According to the IV.2 theorem [34], the sufficient separation condition is Note that Equation (39) is only a sufficient condition in ANM, and always much more conservative in most realistic applications. That is to say, if Equation (39) is indeed not satisfied, one may have a chance to get the DOAs.
The overall implementation of our proposed solution is summarized as Algorithm 1.

Algorithm 1
The proposed undetermined wideband DOA estimation method.

Simulation Results
To show the benefits of our proposed method, we used several simulations to testify its better DOA estimation performance, compared to the existing wideband DOA estimation methods: W-SpSF method [18], spatial smoothing MUSIC (SS-MUSIC) [36], and OGSBL [37]. Moreover, to better demonstrate the estimation performance of these four methods, the Cramer-Rao lower bound (CRLB) for wideband signal provided in [18,38] was used as the benchmark, which represents a lower bound for unbiased estimator. In addition, we used the information theory methods, i.e., AIC, MDL, and AIC-C3, to detect the signal source number.

Simulation Setting
We considered a two-level nest array with six physical sensors, whose inner sub-array and outer sub-array have the same amount of sensors, that is to say, M 1 is equal to M 2 = 3. The physical sensors locate at [0, d 1 , 2d 1 , 3d 1 , 7d 1 , 11d 1 ], and d 1 is equal to 1/2λ, where λ represents the wavelength corresponding to the center frequency f 0 = ( f l + f h )/2. In addition, c = 1500 m/s denotes the signal propagation velocity of underwater acoustic scenarios. Therefore, we have d 1 = l/2λ = c/(2 f 0 ). The unchangeable default simulation parameters setting is shown in Table 1, and a co-array of NA with six sensors is displayed in Figure 2.

Source Number Detection
To verify the detection accuracy of wideband signal with NA, we compared the probability of detection using the mentioned information theory method in the overdetermined and underdetermined cases, respectively. Under the overdetermined scenario, we assumed that five wideband signals impinge on the NA from −30 • , −15 • , 0 • , 15 • , 29 • , and seven wideband signal sources are from −45 • , −30 • , −15 • , 0 • , 15 • , 29 • , 44 • for the underdetermined scenario. We varied the SNR from −15 to 5 dB , and fixed snapshots to 30. Figure 3 shows the result of the probability of detection versus SNR. The detection accuracies of AIC are the highest in low SNR, followed by AIC-C3 and MDL, as the case may be. The detection accuracies of all information theory methods for five signals are higher than those of seven signals. Table 2 displays the needed SNRs of the three methods, when the probability of detection are equal to 0.9 and 1, respectively.  We compared the normalized spatial spectrum of all methods under the overdetermined and underdetermined scenarios, respectively. Firstly, five uncorrelated signal sources were at [−30 • , −15 • , 0 • , 15 • , 29 • ] and the SNR and snapshots were set as −5 dB and 30, respectively. In addition, the grid size was 0.5 degree in W-SpSF method. Figure 4 displays the normalized spatial spectrum for the five signal sources with the NA of six sensors. In Figure 4, the dashed lines are used to represent the true signal sources, while the solid lines stand for the estimated DOAs. It can be seen that five wideband signals can be detected successfully by SS-MUSIC, W-SpSF method, OGSBL, and the proposed method in the low SNR. However, there is a sidelobe peak with the performance of W-SpSF, and worse accuracy for some angles are displayed with the performance SS-MUSIC and OGSBL. Therefore, our proposed method has the best detection and estimation performance as the DOA estimations are closest to the true signals.
Then, we compared the detection performance of all methods for the underdetermined case. The seven signals were the same as in Figure 3. The SNR and snapshots were again fixed as −5 dB and 30, respectively. Figure 5 shows the normalized spatial spectrum for the five signal sources of all methods, respectively. It is shown that all methodd can detect seven peaks successfully in the normalized spatial spectrum nearby the true DOAs. The normalized power of SS-MUSIC is worse for some signals, and there is a certain deviation in this method. W-SpSF method has a sidelobe and a certain deviation, as well. Although OGSBL has stronger normalized power, its DOA estimations are farthest from the true signals. Figure 5 shows that the propose method gives the best detectable capability, since it has strong normalized power and its estimated DOAs are closest to the true signals.  Afterwards, to evaluate the estimation accuracy of our proposed method, we give the estimation results of all methods in the light of the root mean square error (RMSE), that is where θ k andθ kq denote the kth true angle and the kth estimated angle in the qth Monte Carlo trial, respectively. Unless otherwise specified, seven wideband signals impinged on the NA from −45 • , −30 • , −15 • , 0 • , 15 • , 29 • , 44 • . We varied SNR from −10 to 20 dB at the step of 3 dB and fixed snapshots to L = 30. Figure 6 shows the results averaged over Q = 100 Monte Carlo trials as a function of SNR. It can be observed that OGSBL and SS-MUSCI have the same performance, whose estimation accuracies do not improve markedly with the SNR increasing. Although the W-SpSF method outperforms the SS-MUSIC method when SNR > −1 dB, and its RMSE curve is gradually close to the CRLB curve, our proposed method's estimation accuracy is the closest to the counterpart of CRLB. The reason is that the method we propose can recover the optimal covariance matrix with Toeplitz structure, which is beneficial to obtain accurate DOAs. Finally, to further illustrate the estimation performance of our proposed method, the RMSE versus the number of snapshots was examined for all methods. In this simulation, the value of snapshots was varied from 20 to 90, while the SNR was fixed to 5 dB. In Figure 7, we note that all methods can obtain more accurate estimation performance by increasing the number of snapshots. Obviously, the RMSE of OGSBL is still the largest among all methods followed the RMSE of SS-MUSIC (Figure 7). Although the W-SpSF RMSE curve tends to decrease most quickly versus the number of snapshots, its performance is worse than the counterpart of the proposed method in all regions. Therefore, under the case of same snapshots, the proposed method exhibits better performance in comparison with SS-MUSIC, OGSBL, and W-SpSF method.

Conclusions
In this paper, a novel wideband DOA estimation technique equipped with NA is proposed, which seeks and utilizes the Toeplitz structure of the reconstructed virtual array by ANM and thus can work well under the underdetermined wideband case. Our method provides an effective solution through collecting the vectorized covariance matrix via sub-band focusing operation to reformulate the single wideband model of virtual array, recovering the optimal Toeplitz covariance matrix and improving the Root-MUSIC by Toeplitz structure enhancement. Numerical simulations showed that the proposed method outperforms the state-of-the-art methods. In the next work, we will continue to study the wideband coherent signals DOA estimation via ANM with sensors position error.