Non-Iterative Subspace-Based DOA Estimation in the Presence of Nonuniform Noise
M. Esfandiari S. A. Vorobyov S. Aliban M. Karimi




 



IEEE Signal Processing Letters
Shell et al.: Bare Demo of IEEEtran.cls for IEEE Journals




Non-Iterative Subspace-Based DOA Estimation in the Presence of Nonuniform Noise
    Majdoddin Esfandiari, 
	    Sergiy A. Vorobyov, Fellow, IEEE,
        Simin Alibani, 
        and Mahmood Karimi 
M. Esfandiari and S. A. Vorobyov are with Dept. Signal Processing and Acoustics, Aalto University, PO Box 15400, 00076 Aalto, Finland. S. Alibani and M. Karimi are with School of Electrical and Computer Engineering, Shiraz University, Shiraz, Iran.
Emails: majdoddin.esfandiari@aalto.fi; sergiy.vorobyov@aalto.fi; alibani.simin@gmail.com; karimi@shirazu.ac.ir.
Corresponding author is S. A. Vorobyov.

    December 30, 2023
=======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




The uniform white noise assumption is one of the basic assumptions in most of the existing directional-of-arrival (DOA) estimation methods. 
In many applications, however, the non-uniform white noise model is more adequate. Then the noise variances at different sensors have to be also estimated as nuisance parameters while estimating DOAs. In this letter, different from the existing iterative methods that address the problem of non-uniform noise, a non-iterative two-phase subspace-based DOA estimation method is proposed. The first phase of the method is based on estimating the noise subspace via eigendecomposition (ED) of some properly designed matrix and it avoids estimating the noise covariance matrix. In the second phase, the results achieved in the first phase are used to estimate the noise covariance matrix, followed by estimating the noise subspace via generalized ED. Since the proposed method estimates DOAs in a non-iterative manner, it is computationally more efficient and has no convergence issues as compared to the existing methods. Simulation results demonstrate better performance of the proposed method as compared to other existing state-of-the-art methods.




Array processing, direction-of-arrival (DOA) estimation, subspace based methods, nonuniform noise, spectral analysis.







§ INTRODUCTION

Direction of Arrival (DOA) and spectral estimation are the fundamental problems in array processing and spectral analysis with many applications in radar, sonar, navigation and communication systems, as well as acoustic tracking to mention just a few <cit.>–<cit.>. There exist several DOA estimation approaches. Among the most notable are the maximum likelihood (ML), beamforming-based, parametric subspace-based, and sparse representation-based approaches <cit.>–<cit.>. Parametric subspace-based DOA estimation methods, such as multiple signal classification (MUSIC) <cit.>, <cit.> and estimation of signal
parameters via rotational invariance technique (ESPRIT) <cit.> are well known to provide high accuracy and high resolution for DOA estimation with low computational complexity in comparison to the methods such as ML <cit.>. In addition to the traditional far-field narrow band signal assumption, the assumption of uncorrelated sources is also critical for the former methods. A fundamental assumption that applies to all the aforementioned methods is however the presence of spatially uniform white noise. Under this assumption, the analytic concentration of the ML function with respect to the noise variance single parameter becomes possible, while parametric subspace-based methods are just built on this assumption since it enables separation of signal and noise subspaces <cit.>, <cit.>. 

In diverse practical scenarios, the spatially uniform white noise assumption may be violated. Indeed, the sensor noise may be  non-uniform <cit.>, <cit.>–<cit.>, spatially correlated <cit.>, <cit.>, or block-correlated <cit.>. Spatially white non-uniform noise arises when sensor noise powers are non-identical across the array, and leads to diagonal noise covariance matrix with not-identical entries. 
To overcome the problem of performance degradation in the presence of non-uniform noise, a variety of DOA estimation algorithms and techniques have been proposed. In <cit.> and <cit.>, two deterministic and stochastic ML estimators are respectively proposed based on iterative procedures. These two estimators suffer from high computational complexity. Thus, two iterative subspace estimation algorithms with lower complexity, called iterative maximum likelihood subspace estimation (IMLSE) and iterative least squares subspace estimation (ILSSE), respectively, based on ML and least squares (LS) have been proposed in <cit.> for estimating signal subspace and noise covariance matrix. These algorithms then use spectral MUSIC method for performing the DOA estimation. 

In this letter, we propose a new subspace-based method for DOA estimation in spatially non-uniform noise, which is non-iterative, thus leading to lower computational complexity and avoiding any convergence issues. The method has two phases. In the first phase, the noise subspace is initially estimated via eigendecomposition (ED) of some properly designed matrix without knowing the noise covariance matrix. In the second phase, the noise covariance matrix is then estimated by exploiting the results of the first phase and then the generalized ED is applied to the output array covariance matrix and noise covariance matrix. Simulation results demonstrate the efficiency and superiority of the proposed method in terms of both performance and complexity over the existing methods. 




§ SIGNAL MODEL

Consider an array of M sensors receiving L (L<M is known) independent narrowband signals impinging from the sources in far-field. The signal observed at time instance t by the array is given as

𝐱(t)=𝐀(θ)𝐬(t)+𝐧(t)
   
where 𝐀 (θ) ≜ [𝐚(θ_1), ⋯, 𝐚(θ_L)] is the M× L matrix whose columns are the signal steering vectors 𝐚(θ_i), i=1,⋯, L, θ≜ [θ_1, ⋯, θ_L]^T is the vector of unknown source DOA's, 𝐬(t) is the L× 1 vector of source signals, 𝐧(t) is the M× 1 vector of zero-mean spatially and temporally white Gaussian noise, N is the number of snapshots, t is the discrete time index, and [ · ]^T denotes the transpose. 

Using (<ref>), the array output covariance matrix can be expressed as

𝐑 ≜E{𝐱(t) 𝐱^H(t)} = 𝐀 (θ) 𝐏 𝐀^H (θ) + 𝐐
  
where 𝐏≜ E{𝐬(t)𝐬^H(t)} and 𝐐≜ E{𝐧(t)𝐧^H(t)} are respectively the L× L signal and M× M noise covariance matrices, and E{·} and (·)^H  denote the expectation and Hermitian transpose operators, respectively. For uncorrelated sources, 𝐏 is a diagonal matrix. Under the non-uniform uncorrelated noise assumption, 𝐐 is also a diagonal matrix of the form

𝐐=diag { σ_1^2,⋯, σ_M^2 } 
  
where σ_m^2, m=1, ⋯, M are the sensor noise variances and diag{·} denotes a diagonal matrix.




§ NEW PROPOSED METHOD

In parametric subspace-based DOA estimation, the noise subspace needs to be first estimated. Under the non-uniform noise, also the noise covariance matrix possibly needs to be estimated. Then the orthogonality of the noise subspace basis vectors and the source steering vectors can be used to estimate the source DOA's.  



To estimate the noise subspace, recall that in (<ref>), 𝐀^H(θ) is an L× M full row-rank matrix, and there are M-L orthonormal vectors 𝐮_l, l=1,⋯, M-L satisfying the following homogenous equation

𝐀^H(θ)𝐮_l=0
      
where 0 is the vector of zeros. Multiplying both sides of (<ref>) on the right by 𝐮_l and using (<ref>), we obtain 

𝐑𝐮_l = 𝐀 (θ) 𝐏 𝐀^H (θ) 𝐮_l + 𝐐 𝐮_l = 𝐐 𝐮_l,   l=1,⋯, M-L.

  
According to (<ref>), M-L vectors 𝐮_l, l=1,⋯, M-L which span the range space of the noise subspace are the solutions of the generalized ED problem for the matrices 𝐑 and 𝐐, while all M-L eigenvalues being equal to one. However, since the matrix 𝐐 is unknown, 𝐮_l, l=1,⋯, M-L cannot be computed as simple as in the uniform noise case. 

We observe, however, that the array output covariance matrix 𝐑 can be written as the following sum of two matrices

𝐑 = 𝐑_1+𝐑_2
   
where

[𝐑_1]_i,j = 
[𝐑]_i,j,       i≠j 

0,       i=j

 
and

𝐑_2    = diag { [𝐑]_1,1, ⋯, [𝐑]_M,M } 

   =diag { ∑_k=1^L s_k + σ_1^2, ⋯, ∑_k=1^L s_k +σ_M^2 } 
 
with s_k being the received power of the kth source. 

Substituting (<ref>), (<ref>), and (<ref>) into (<ref>), we obtain

𝐑_1𝐮_l = (𝐐 - 𝐑_2) 𝐮_l = - ( ∑_k=1^L s_k ) 𝐮_l .
 
It can be seen from (<ref>) that 𝐮_l, l=1,⋯, M-L can be obtained by applying ED to the matrix 𝐑_1 only. As a matter of fact, since adding the scaled identity matrix of the form b ·𝐈 to 𝐑_1 does not alter eigenvectors and only shifts eigenvalues of 𝐑_1, the noise subspace basis vectors 𝐮_l, l=1,⋯, M-L can be computed by applying ED to every matrix whose diagonal elements are identical and off-diagonal elements are equal to the off-diagonal elements of 𝐑. Thus, after applying ED to 𝐑_1, the noise subspace basis vectors 𝐮_l, l=1,⋯, M-L, or in matrix form 𝐔≜ [𝐮_1,⋯,𝐮_M-L], can be obtained even without the need to estimate 𝐐. This novel result is stated and proved in the following lemma. 

Lemma 1: The noise subspace basis vectors 𝐮_l, l=1,⋯, M-L are the M-L eigenvectors of 𝐑_1 whose corresponding eigenvalues are the smallest.

Proof: Assume that there exists an M × 1 vector 𝐝 that satisfies the following conditions

    𝐀(θ)𝐏𝐀^H(θ)𝐝≠0  
 
    𝐑_1𝐝=λ𝐝 .
 
Adding 𝐑_2 𝐝 to both sides of (<ref>) and using (<ref>), we obtain

𝐑_1𝐝+𝐑_2𝐝=𝐑𝐝=λ𝐝+𝐑_2𝐝
.  

Inserting (<ref>) into (<ref>), and rearranging the terms, we have

𝐀(θ)𝐏𝐀^H(θ)𝐝=λ𝐝+(𝐑_2-𝐐)𝐝 .
 
Moreover, substituting (<ref>) and (<ref>) into (<ref>), we obtain

𝐀(θ)𝐏𝐀^H(θ)𝐝 = ( λ+ ∑_k=1^L  s_k ) 𝐝 .
 

It can be seen from (<ref>) that 𝐝 is an eigenvector of the matrix 𝐀 (θ) 𝐏𝐀^H (θ) while its corresponding eigenvalue is equal to λ + ∑_k=1^L s_k. Since the matrix 𝐀 (θ) 𝐏𝐀^H (θ) has L positive eigenvalues and the condition (<ref>) has to be satisfied, it can be concluded that

λ+ ∑_k=1^L s_k > 0   ⇒  λ>-∑_k=1^L s_k .
 
In other words, (<ref>) indicates that -∑_k=1^L s_k is the lower bound on the smallest eigenvalue of 𝐑_1. Thus, the noise subspace basis 𝐮_l, l=1,⋯, M-L is composed of M-L eigenvectors of 𝐑_1 with the smallest eigenvalues. □

Knowing 𝐔, the spectral-MUSIC method, for example, can be used for the source DOA's estimation by finding the locations of L peaks in the pseudo-spectrum

𝐒(θ)=1/𝐚^H(θ)𝐔𝐔^H𝐚(θ) .
 
However, the estimate of 𝐔 can be further improved using the results of the initial 𝐔 estimation and	exploiting (<ref>). Indeed, if there exists an estimate of 𝐐, (<ref>) can be solved by applying generalized ED to the matrices 𝐑 and 𝐐. Then more accurate noise subspace basis vectors 𝐮_l, l=1,⋯, M-L can be found as stated in the following lemma.

Lemma 2: The noise subspace basis vectors 𝐮_l, l=1,⋯, M-L are the M-L eigenvectors, obtained by applying generalized eigendecomposion to the matrices 𝐑 and 𝐐 whose corresponding eigenvalues are the smallest.

Proof: Similar to the proof of Lemma 1, assume that there exists an M × 1 vector 𝐝 that satisfies (<ref>) and 𝐑𝐝=λ𝐐𝐝. 

Inserting (<ref>) into 𝐑𝐝=λ𝐐𝐝, and rearranging the terms, we have

𝐀(θ)𝐏𝐀^H(θ)𝐝=(λ-1)𝐐𝐝 .
 
Multiplying both sides of (<ref>) on the left by 𝐝^H, we obtain

𝐝^H𝐀(θ)𝐏𝐀^H(θ)𝐝=𝐝̅^H𝐏𝐝̅=(λ-1)𝐝^H𝐐𝐝 
 
where 𝐝̅=𝐀^H(θ)𝐝. Since both 𝐏 and 𝐐 are positive definite matrices and the condition (<ref>) has to be satisfied, it can be concluded that

λ- 1 > 0   ⇒  λ>1 .
 
In other words, (<ref>) indicates that 1 is the lower bound on the smallest eigenvalue of generalized ED of 𝐑 and 𝐐. Thus, the noise subspace basis 𝐮_l, l=1,⋯, M-L is composed of M-L eigenvectors with the smallest eigenvalues. □

To apply Lemma 2, we first need to find an estimate of 𝐐. Towards this end, first we write 𝐐 as the following sum of two diagonal matrices

𝐐 = σ^2 𝐈 + 𝐐_nun
 
where σ^2 represents the common part of sensor noise powers, which is computed later, and 𝐐_ nun is a diagonal matrix whose diagonal elements, except for one of them, are nonzero. The place of this zero element is the place of the smallest diagonal element of 𝐑. As a result, the rank of 𝐐_ nun is M-1, and  

𝐞^T_k 𝐐_nun = 0 
 

where 𝐞_k is the M × 1 unit vector such that 

[𝐞_k]_i = 
0,       i≠k 

1,       i=k



and k is the index of the smallest diagonal element of 𝐑. 

Multiplying both sides of (<ref>) by 𝐞_k^T on the left and using (<ref>) and (<ref>), we obtain

𝐞_k^T 𝐑 𝐮_l = 𝐞_k^T (σ^2 𝐈 + 𝐐_nun) = σ^2 𝐞_k^T 𝐮_l .
 
Equation (<ref>) can be written for all vectors 𝐮_l, l=1,⋯, M-L in the following matrix-vector form

𝐞_k^T 𝐑 𝐔 = σ^2 𝐞_k^T 𝐔
 
where 𝐔 is composed of 𝐮_l, l=1,⋯, M-L obtained by Algorithm 1. Consequently, σ^2 can be computed as

σ^2 = 𝐞_k^T 𝐑 𝐔 𝐔^H 𝐞_k/𝐞_k^T 𝐔 𝐔^H 𝐞_k .
 

The only issue remaining is the construction of the matrix 𝐐_ nun. Let us set the nonzero diagonal elements of 𝐐_ nun as the differences of the corresponding elements in 𝐑 with the smallest diagonal element of 𝐑, that is,

𝐐_nun = diag { [𝐑]_1,1 - c, ⋯, [𝐑]_M,M - c }
  
where c is the smallest diagonal element of 𝐑.

Finally, the matrix 𝐐 can be formed by utilizing (<ref>), (<ref>), and (<ref>). With the matrices 𝐑 and 𝐐, the noise subspace basis 𝐮_l, l=1,⋯, M-L can be re-estimated as stated in (<ref>) and Lemma 2. Then the new 𝐔 can be formed and the source DOA's can be estimated by finding, for example, the locations of L peaks in (<ref>). The corresponding algorithm for DOA estimation in non-uniform noise is summarized in Algorithm 1, where the sample data covariance matrix 𝐑̂ is used as an estimate of the array output covariance matrix 𝐑. Steps 1 and 2 represent the first phase of the algorithm that can be followed by step 7 directly. Steps 3–6 represent the second correction phase of the algorithm.



Complexity analysis: For the proposed method, the ED of 𝐑_1 or the generalized ED of 𝐑 and 𝐐 are involved. The corresponding complexity is O(M^3) <cit.>. It is equivalent to the complexity in each iteration of IMLSE or ILSSE <cit.>. The difference arises from the fact that IMLSE and ILSSE are iterative methods (the number of iterations can be comparable to M in many scenarios to converge to their best result).



§ SIMULATION RESULTS

A ULA with M=8 omnidirectional sensors, which are separated by half wavelength, is considered. Two far-field uncorrelated narrow-band signals (with equal powers) impinge on the array simultaneously from θ_1=-3^∘ and θ_2=6^∘, respectively.  

The worst noise power ratio (WNPR), and the signal-to-noise ratio (SNR) are defined as
WNPR = σ^2_ maxσ^2_ min and 
SNR = σ^2_ sM∑_m=1^M (σ^2_m)^-1, respectively, where σ^2_ max and σ^2_ min are the maximum and minimum sensors noise powers, respectively, and σ^2_ s is the signal power. The number of snapshots (N) and the number of Monte Carlo trials (K) are set to 500 and 5000, respectively. The root mean squared error (RMSE) of DOA's estimation is defined as

RMSE = √(1/KL∑_k=1^K ∑_l=1^L (θ̂_k,l-θ_l)^2)

where θ̂_k,l denotes the lth DOA estimate in the kth trial. 

To validate the performance of the proposed methods, two examples are considered, and the results are compared to the performance of the standard spectral-MUSIC and well as IMLSE and ILSSE methods both after the first iteration only and also after convergence is achieved.[The performance of ILSSE is plotted only for the SNRs when ILSSE converges for all trials.] The number of sources is assumed to be known for all methods tested.

Example 1: The non-uniform noise covariance matrix is fixed in all simulation runs and is given as 𝐐=diag{1, 1, 1, 1, 1, 20, 30, 50}, resulting in WNPR=50. Fig. <ref> shows the RMSEs for the methods tested versus SNR. The Cramer-Rao bound (CRB) <cit.> is also shown. It can be seen from the figure that the proposed method after both phases demonstrates the best threshold behavior, despite it is non-iterative and has lower computational complexity. Although the IMLSE method after the first iteration only possesses as low computational complexity as the proposed method, its DOA estimation accuracy is very poor. Thus, multiple iterations have to run for it to achieve its best result.


  

Example 2: To provide more comprehensive insights into the performance of the methods tested, the noise covariance matrix is chosen in this example multiple times randomly with maximum WNPR of 30. Then the RMSE results are also averaged over 50 different realizations of the noise covariance matrix 𝐐 for which of each 5000 Monte Carlo trials are averaged. The average RMSEs for the methods tested are shown in Fig. <ref>.[Since the curves in Fig. 2 are resulted also from averaging over different realization of the nuisance parameters, i.e., 𝐐's with different WNPRs, the CRB is not applicable in the setup and is not shown. Indeed, fixed parameters of interest as well as nuisance parameters have to be assumed for CRB.] It can be seen from the figure that the second phase of the proposed method shows the best performance and improves the threshold behavior by about 2 dB as compared to the next best performing method that is the IMLSE method. 

Furthermore, Table I shows the average run time of the methods tested for SNRs=-5, 0, 5, 10 , and 15 dB for the setup of Example 2. The simulation is performed on a PC running an Intel(R) Xeon(R) 3.40GHz CPU. It can be observed that the proposed method is superior in terms of the required time which is reduced by orders of magnitude compared to the existing methods. The average number of iterations for IMLSE and ILSSE are about 27 and 52, respectively, for SNR=5 dB, for example. Inspecting Fig. 2 and Table I together, it can be seen that even the first phase of the proposed method leads to superior performance with lower complexity compared to the competitive methods after the first iteration.






§ CONCLUSION

A novel computationally efficient non-iterative two-phase subspace-based parametric method for DOA estimation in the presence of unknown spatially non-uniform noise has been proposed. The noise subspace estimation problem is converted in the first phase of the proposed method to the problem of finding eigenvectors of a properly designed matrix so that the noise covariance matrix estimation is evoided. In the second phase, the covariance matrix of the non-uniform noise is first estimate based on the results of the first phase and then it is used for finding the noise subspace more accurately by means of generalized ED of this matrix and the data covariance matrix. It is of importance that the proposed method has low computational complexity and is non-iterative, and thus, has no issues with convergence. Moreover, it is superior to the existing iterative state-of-art methods in both the performance and especially the computational cost. 





22
r14a
H.L. Van Trees, Optimum Array Processing, Wiley, NY, 2002.

r14b
A. M. Zoubir, M. Viberg, R. Chellappa, and S. Theodoridis, Eds., Array and Statistical Signal Processing, Academic Press Library in Signal Processing, Academic Press, vol. 3, 2014.

r14
B. Friedlander and T. Tuncer, Classical and Modern Direction-of-Arrival Estimation.  New York: Academic Press, 2009.

r15a
P. Stoica and R. Moses, Introduction to Spectral Analysis, Prentice Hall, 1997.

r15
Y. Wang, J. Li and P. Stoica, Spectral Analysis of Signals. Morgan and Claypool, 2005.

r1
P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramer-Rao bound,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 5, pp. 720–741, May 1989.

r2
M. Pesavento and A. Gershman, “Maximum-likelihood direction-of-arrival estimation in the presence of unknown nonuniform noise,”  IEEE Transactions on Signal Processing, vol. 49, no. 7, pp. 1310–1324, July 2001.

r6
R.O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Transactions on  Antennas and Propagation, vol. 34, no. 3, pp. 276–280, Mar. 1986.

r7
R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 7, pp. 984-995, July 1989.

r8
D. Malioutov, M. Cetin, and A.S. Willsky, “ A sparse signal reconstruction perspective for source localization with sensor arrays,”
IEEE Transactions on Signal Processing, vol. 53, no. 8, pp. 3010–3022, Aug. 2005.

Mahdi
M. Shaghaghi and S. A. Vorobyov, “Subspace leakage analysis and improved DOA estimation with small sample size,” IEEE Transactions on Signal Processing, vol. 63, no. 12, pp. 3251–3265, June 2015.

Matt
M.W. Morency, S.A. Vorobyov, and G. Leus, “Joint detection and localization of an unknown number of sources using algebraic structure of the noise subspace,” IEEE Transactions on Signal Processing, vol. 66, no. 17, pp. 4685–4700, Sept. 2018. 

r3
C. Chen, F. Lorenzelli, R. Hudson, and K. Yao, “Stochastic maximum-likelihood DOA estimation in the presence of unknown nonuniform noise,” IEEE Transactions on Signal Processing, vol. 56, no. 7, pp. 3038–3044, July 2008.

r5
B. Liao, S. Chan, L. Huang, and C. Guo, “Iterative methods for subspace and DOA estimation in nonuniform noise", IEEE Transactions on Signal Processing, vol. 64, no. 12, pp. 3008–3020, June 2016.

r11
B. Liao, C. Guo, L. Huang,  and J. Wen, “Matrix completion based direction-of-arrival estimation in nonuniform noise,” in Proc. Intern. Digital Signal Processing Conf. 2016, pp. 66–69.

r12
B. Liao, G. S. Liao, and J. Wen, “A method for DOA estimation in the presence of unknown nonuniform noise,” Journal of Electromagnetic Waves and Applications, vol. 22, no. 14-15, pp. 2113–2123, 2008.

r12a
B. Liao and S. C. Chan, “A simple method for DOA estimation in the presence of unknown nonuniform noise,” in  Proc. IEEE Conf. Acoustics, Speech and Signal Processing, 2015.

r9
M. Agrawal and S. Prasad, “A modified likelihood function approach to DOA estimation in the presence of unknown spatially correlated Gaussian noise using a uniform linear array,” IEEE Transactions on Signal Processing, vol. 48, no. 10, pp. 2743–2749, Oct. 2000.

r10
M. Malek-Mohammadi, M. Jansson, A. Owrang, A. Koochakzadeh, and M. Babaie-Zadeh, “DOA estimation in partially correlated noise using low-rank/sparse matrix decomposition,” in Proc. IEEE Workshop Sensor Array and Multichannel Signal Processing, 2014.

r2a
S. A. Vorobyov, A. B. Gershman, and K. M. Wong, “Maximum likelihood direction of arrival estimation in unknown noise fields using sparse sensor arrays,”  IEEE Transactions on Signal Processing, vol. 53, no. 1, pp. 34–43, Jan. 2005.

r20
G. H. Gloub and C. F. Van Loan, Matrix Computations.  Baltimore: John Hopkins University Press, 1996.




