Statistics on noise covariance matrix for covariance fitting-based compressive sensing direction-of-arrival estimation algorithm: For use with optimization via regularization

A direction-of-arrival (DOA) estimation algorithm based on a sparse representation of sensor measurements with a basis composed of samples from the array manifold (Das and Gerstoft, 2017; Das, 2017; Edelmann et al., 2013; Edelmann and Gaumond, 2011; Gaumond and Edelmann, 2013; Hawes et al., 2017; Liu et al., 2017; Xenaki and Gerstoft, 2016; Yardim et al., 2014; Shen et al., 2017a; Shen et al., 2017b; Tian et al., 2017; Wu et al., 2018; Xenaki and Gerstoft, 2014, 2015) is considered in this paper.

Many studies have been conducted on DOA estimation based on compressive sensing (CS) (Malioutov et al., 2005, 2011; Zheng and Kaveh, 2013). These algorithms have been widely employed due to their advantages over the conventional DOA estimation algorithms such as the conventional beamforming algorithm, the MUSIC algorithm, and the Capon beamforming algorithm:

• CS-based DOA estimation methods (Malioutov et al., 2005) show less performance degradation due to incorrect estimation of the number of incident signals in comparison with the MUSIC algorithm.

• CS-based DOA estimation methods (Malioutov et al., 2005, 2011; Zheng and Kaveh, 2013) are superior to the conventional beamforming algorithm in terms of the resolution and the main beamwidth.

• CS-based DOA estimation method (Malioutov et al., 2005) shows less performance degradation than the MUSIC algorithm for coherent incident signals.

A data fitting algorithm using single snapshot measurement has been proposed to recover the signal based on sparse measurement (Malioutov et al., 2005). In the algorithm, the measurement vector, whose entry is associated with each sensor, is recovered via sparse recovery.

Single measurement formulation has been extended to the multiple measurement formulation when multiple measurements are available (Malioutov et al., 2005). In a data fitting algorithm using multiple measurements, the data matrix, whose rows are associated with the sensors and whose columns are associated with the snapshots, should be recovered. The computational complexity of the multiple snapshot algorithm can be extremely large when the number of the snapshots is large. To alleviate this problem, the singular value decomposition (SVD) has been applied to the multiple snapshot algorithm (Malioutov et al., 2005).

Instead of recovering the data vector or the data matrix, the signal covariance matrix can be recovered using the information of the array covariance matrix and the array manifold (Malioutov et al., 2011; Zheng and Kaveh, 2013). In the covariance fitting algorithm, the computational complexity does not greatly increase with an increase in the number of snapshots since the number of entries to be recovered is independent of the number of snapshots.

Both the data fitting algorithm and the covariance fitting algorithm can be implemented via regularization, and the parameter of the regularization should be determined. How the regularization parameter in the data fitting algorithm can be chosen has been presented in Malioutov et al. (2005). A regularization parameter selection scheme for the covariance fitting algorithm has been presented in Zheng and Kaveh (2013). However, the scheme in Zheng and Kaveh (2013) is only valid when the number of snapshots is very large.

In this paper, we propose how the regularization parameter for the covariance fitting algorithm can be chosen. The proposed scheme can be applied not only for a large number of snapshots but also for a small number of snapshots. That is, the proposed scheme works irrespective of the number of the snapshots, which will be shown in the numerical results.

The performance of the sparse spectrum fitting (SpSF) algorithm is highly dependent on the regularization parameter in the cost function of the SpSF algorithm. Therefore, the regularization parameter should be properly chosen for successful implementation of the algorithm. The scheme presented in the previous study (Zheng and Kaveh, 2013) is only applicable when the number of snapshots is very large, which implies that the scheme is not generally used irrespectively of the number of snapshots. The scheme, which can be applied for any number of the snapshots, will be proposed in this paper.

The data matrix $x(t)$ received from M sensor array can be defined as

where L is the number of snapshots and $Nθ$ is the number of search angles for θ. The entries in Eq. (4) are complex Gaussian random variables whose real and imaginary parts are zero-mean Gaussian distributed with variance $σ2/2$⁠. The steering matrix A can be defined as follows:

Due to $Nθ≫M$⁠, Eq. (1) is an underdetermined system. $Nθ$ is much greater than P, which is the number of incident signals.

When considering multiple snapshots, matrices can be defined as follows:

Using these matrices, Eq. (1) can be extend to

The sparse S can be obtained from Eq. (10) via sparse recovery and we can estimate the incident angles from the indices of S where the non-zero entries are located.

The array covariance matrix can be written as

where $Rs$ is the signal covariance matrix. Note that $X$ is given in Eq. (10). Equation (11) is an array covariance matrix obtained from the ensemble average. In practice, it is impossible to obtain the ensemble average, and it is replaced by the array covariance matrix obtained from the time average

If the number of snapshots is particularly large, the array covariance matrices from the time average and that from the ensemble average are identical. In Eq. (12), $R̂s$ is the estimate of the signal covariance matrix $Rs$ and $E$ is the estimate of the noise covariance matrix.

When the incident signals are uncorrelated, Eq. (12) can be written as

The variables in (13) are defined as follows:

$vec(•)$ is the operator that creates a long column vector by connecting the columns of the matrix in order. $p°$ is a column vector composed of the diagonal components of the covariance matrix of the signal covariance matrix $[p°=diag(Rs)]$⁠. Since the signals are assumed to be uncorrelated signals, all the components except the diagonal entries of the signal covariance matrix are zero. Thus, we only have to estimate the diagonal components. Therefore, each entry of $p°$ implies the power for each signal. The cost function for determination of $p°$ can be defined as

If the number of snapshots is close to infinity, the analytic value of $λ$ is given by (Zheng and Kaveh, 2013)

$p$ can be estimated via sparse recovery by adopting Eq. (20) in Eq. (19).

For correlated incident signals, the off-diagonal entries of the signal covariance matrix $Rs$ in Eq. (12) as well as the diagonal entries are non-zero. Unlike the uncorrelated case where only the diagonal entries of the covariance matrix should be estimated, all the entries of the signal covariance matrix should be estimated. The cost function for estimating $Rs$ can be written as

If the number of snapshots is close to infinity, the $λ$ value in Eq. (20) can be used.

When the number of snapshots is finite, the optimal regularization parameter $λ$ of Eqs. (19) and (21) should be dependent on the number of snapshots. In this section, for finite snapshots, the criterion for optimal regularization parameter is presented. For uncorrelated incident signals, the criterion for the optimal regularization parameter is obtained from Eq. (13),

When the incident signals are correlated, from Eq. (12) it can be easily shown that

In order to estimate the DOAs through Eqs. (19) and (21), we need to estimate the signal covariance matrix in the viewpoint of fitting the signal covariance matrix to the array covariance matrix, which accounts for the former term of Eq. (21), and in the view point of the spatial sparsity of the signal covariance matrix, which accounts for the latter term of Eq. (21), for various $λ$ values. For the determination of an optimal value of the regularization parameter $λ$⁠, the optimization in Eqs. (19) or (21), depending on whether the incident signals are correlated or uncorrelated, should be repeated for different positive values of $λ$⁠. The recovered signal covariance matrix satisfying Eq. (22) for uncorrelated incident signals, or Eq. (23) for correlated incident signals, is chosen to be the optimal signal covariance matrix. The regularization parameter associated with the optimal signal covariance matrix is the optimal regularization parameter. The DOAs can be estimated from the indices of the nonzero entries of the optimal signal covariance matrix.

To implement the above scheme, the mean of $‖E‖f2$⁠, which is equal to $‖Ev‖22$⁠, should be available. The mean of $‖E‖f2$ can be derived as follows. $‖E‖f2$ can be expressed as

Using Eq. (24), the square of the diagonal components, $(Eii)2$⁠, can be explicitly written, in terms of an additive complex Gaussian random noise on the sensor arrays

In Eq. (25), the mean of $|ni(tl)|4$ can be written as

The mean of $|ni(tl)|2|ni(tl′)|2$ can be written as

where the first equality is due to the fact that the noises at different sampling times are uncorrelated. Substituting Eqs. (26) and (27) in Eq. (25) results in

$|Eij|2$ in Eq. (24) can be more explicitly written as

The expectation of the second summation in Eq. (29) is identically zero since the noises at different sampling times and the noises at different sensors are uncorrelated. Therefore, the expectation of $|Eij|2$ is reduced to

Substituting Eqs. (28) and (30) in Eq. (24) results in

Equation (31) is used in Eq. (23) for the determination of the optimal $λ$ value, and $Rs(λ)$ in Eq. (23) is estimated by optimizing the cost function in Eq. (21) for different $λ$ values. Our derivation presented in this section is different from Eq. (20) in the previous study (Zheng and Kaveh, 2013) in that we explicitly take into account the number of snapshots in our derivation since the covariance matrix in Eq. (12) is obtained from finite L snapshots of data received on the array elements. See Fig. 1.

$R̂$⁠, which is an estimate of array covariance matrix, is obtained from the time average of the received array data. Once discrete search angles are determined, a matrix $A$ whose each column corresponds to an array vector for each search angle is constructed. For each regularization parameter, the signal covariance matrix $Rs(λ)$ can be obtained via sparse recovery algorithm. Choose $Rs(λ)$ for which $‖R−ARs(λ)AH‖f2$ is very close to $‖E‖f2$⁠.

Note that the matrix $A$ is dependent on how the discrete search angles are selected since the array vectors for all the search angles constitute the columns of $A$⁠. Finally, the directions-of-arrival are estimated from the nonzero indices of the diagonal elements of the finally chosen R_s since diagonal elements of the signal covariance matrix can be interpreted as power spectrum for search angles.

In this section, numerical results are presented to validate the propose scheme. We tabulate the various parameters for the simulation in Table I.

There are three different simulations in this section. The sensor array geometry is Uniform Linear Array (ULA). Candidate regularization parameters are logarithmically equally spaced between 10⁰ and 10². The number of candidate regularization parameters is 20. Figures 2 and 3 show the spectra of the SpSF algorithm when the number of sensor is 8 for various number of snapshots. Figure 4 shows the performances of the SpSF with the optimal regularization parameter and other regularization parameters when the number of sensors is 5 and the number of snapshots is 600. Figures 5–7 show the statistical performance of the SpSF with the optimal regularization parameter and other regularization parameters when the number of sensors is 8 and the numbers of snapshots are 400, 600, and 800.

Figures 2 and 3 show the spectra of the SpSF algorithm when the number of snapshots is changed from 100 to 10 000 for signal-to-noise ratio (SNR) $=0$ dB and SNR $=5$ dB, respectively. When the number of snapshots is small, it is shown that there are spurious peaks around −20 dB. In case of the spectrum of the conventional DOA estimation algorithm, the difference of the noise floor to the peak is less than 40 dB. In the case of the SpSF, the difference between the noise floor and the peak is more than 100 dB. Therefore, it is shown that the spurious peaks which occur when the number of snapshots is small and do not have a significant influence on the DOA estimation. As shown in Figs. 2(c) and 3(c), the number of spurious peaks decreases and the sharp spectrum is generated as the number of snapshots increases. It can be also confirmed that the SpSF algorithm for the optimal $λ$ accurately estimates the true DOAs.

Figure 4(a) shows the result of using the optimal regularization parameter obtained from the method presented in the paper. Figures 4(b) and 4(c) show the result of using other regularization parameters.

In Fig. 4(a), the optimal regularization parameter is used and this parameter improves the DOA estimation accuracy of the SpSF. Therefore, two sharp peaks are observed and the difference between peak values and noise floor is greater than 100 dB.

In Fig. 4(b), an arbitrary regularization parameter smaller than the optimal regularization parameter is used, and this parameter cannot improve the DOA estimation accuracy of the SpSF as much as the optimal regularization parameter. Therefore, the difference between peak values and noise floor is smaller than the difference value in Fig. 4(a). In this case, if there are additional obstacles in DOA estimation, the performance of DOA estimation seems to be greatly degraded.

In Fig. 4(c), an arbitrary regularization parameter greater than the optimal regularization parameter is used, and the performance with this regularization parameter is inferior to that with the optimal regularization parameter. A large regularization parameter in the cost function corresponds to an emphasis on sparsity rather than model fitting. Due to low SNR and a few numbers of sensors, the performance has deteriorated. In this case, if there are additional obstacles in DOA estimation, the performance of DOA estimation seems to be greatly degraded.

In Figs. 5–7, the statistical performance of the SpSF with respect to SNR is shown. In case of the optimal regularization parameter, the value determined from the proposed scheme is used. In case of the other two values, it is arbitrary values other than the best regularization parameter.

Figures 5–7 show the mean square error (MSE) with respect to SNR when the optimal regularization parameter is selected and the MSE with respect to the SNR when the 16th and 19th values of candidate regularization parameters are selected. Two arbitrary non-optimal values are chosen to compare the performance for these values with that for the optimal regularization parameter. In Figs. 5–7, the performance with optimal regularization parameter is much better than that with non-optimal regularization parameter.

Figures 5(a), 6(a), and 7(a) show the MSE of the SpSF for the first source. The SNR varies from −15 dB to 5 dB. It is shown in Figs. 5(a), 6(a), and 7(a) that the DOA estimation performance of SpSF for the optimal regularization parameter is superior to the DOA estimation performance of the SpSF for two non-optimal regularization parameters. The same is true for Figs. 5(b), 6(b), and 7(b), where the MSEs of the second incident signal are illustrated.

In this paper, a scheme for determining a regularization parameter for the SpSF algorithm for a general number of snapshots, has been presented. The insight on the meaning of the regularization parameter in the formulation is a relative weight between the fitting of the signal covariance matrix to array covariance matrix and the sparsity of the signal covariance matrix. Our derivation exploits the statistical property of complex Gaussian noise incident on array sensors. The fact that additive complex Gaussian noises on different sensors or different snapshots are uncorrelated is also used in the derivation. The feasibility of the proposed scheme is shown in the numerical results.

This work was supported by the Agency for Defense Development (ADD) in Korea under the Contract No. UD160015DD.