Semiparametric Stochastic CRB for DOA Estimation in Elliptical Data Model
^{†}^{†}thanks: The work of Stefano Fortunati has been partially supported by the Air Force Office of Scientific Research under award number FA95501710065.
Abstract
This paper aims at presenting a numerical investigation of the statistical efficiency of the MUSIC (with different covariance matrix estimates) and the IAAAPES Direction of Arrivals (DOAs) estimation algorithms under a general Complex Elliptically Symmetric (CES) distributed measurement model. Specifically, the density generator of the CESdistributed data snapshots is considered as an additional, infinitedimensional, nuisance parameter. To assess the efficiency in the considered semiparametric setting, the Semiparametric Stochastic CramérRao Bound (SSCRB) is adopted as lower bound for the Mean Square Error (MSE) of the DOA estimators.
I Introduction
Estimation of the Direction of Arrival (DOA) of a certain number of sources using an array of active or passive sensors is a standard signal processing problem. There is a huge literature on this topic covering different aspects: parametric and nonparametric estimation methods, deterministic and random signal models and so on (see e.g. [1, 2] and the references therein). Along with the statistical signal models and the related estimation methods, considerable efforts have been dedicated to the derivation of suitable performance bounds for DOA estimation (see e.g. the standard references [3, 4] and the overview provided in [5]). However, most of the works on array processing and, in particular, the ones on lower bounds, assume a Gaussian model for the collected data snapshots. A valuable exception is represented by the paper [6], where the signal model is statistically characterized using the set of the Complex Elliptically Symmetric (CES) distributions [7]. The CES class is a wide family of nonGaussian distributions that encompasses the Gaussian, the Generalized Gaussian, the distribution and all the Compound Gaussian distributions as special cases. As detailed in [7] and [8], many experimental evidences have shown their ability to characterize the heavytailed behaviour of real datasets collected in different applications such as radar/sonar, indoor/outdoor wireless communications or seismic data processing. Note that, in all these examples, the DOA estimation is a key aspect.
The main goal of this paper is then to provide a numerical investigation of the statistical efficiency of some of the most widely used DOA estimation algorithms by dropping the classical Gaussianity assumption in favour of a more general model in which i) the collected snapshots are assumed to be CESdistributed with unknown density generator and ii) the signal and disturbance components are uncorrelated. To handle the additional unknown infinitedimensional parameter, that is the desity generator characterizing the actual CES data distribution, we exploit our recent findings on the Semiparametric Stochastic CramérRao Bound (SSCRB) [9].
To make the paper as selfcontained as possible, Section II provides a very short introduction on CES distribution. The assumed measurement model and the related SSCRB on the estimation of the DOAs of narrowband sources are introduced in Section III. Sections IV and V discuss the DOA estimation algorithms while their efficiency with respect to (w.r.t.) the SSRCB is investigated in Section VI. Some concluding remark is collected in Section VII.
Ii Preliminaries: essentials on CES distributions
This section presents some basic properties of CES distributed random vectors that will be used in the remaining of the paper. For a complete and insightful discussion on CES distributions, we refer the reader to [7] and references therein.
An dimensional CESdistributed random vector is fully characterized by its mean vector , its scatter matrix and its density generator , where is a suitable set of functions. Under the absolutely continuous assumption, i.e. when the scatter matrix has full rank, the pdf of a CESdistributed vector is:
(1) 
Moreover, satisfies the circularity property, i.e. . Any CESdistributed vector admits the following representation:
(2) 
where is a complex random vector uniformly distributed on the unit complex sphere and is the socalled 2ndorder modular variate, such that (s.t.):
(3) 
whose pdf is given by:
(4) 
From (2) and by exploiting the properties of [7, Lemma 1], we have that the covariance matrix of the CESdistributed vector is .
In order to remove the wellknown scale ambiguity, we impose a constraint on the functional form of the density generator . Following the same procedure adopted in [10, 9], we assume that is parameterized in order to satisfy the constraint:
(5) 
As a consequence of (5), the scatter matrix equates the covariance matrix of [7, Sec. III.C]. For further reference, we define the set as the set of all the density generators satisfying the constraint in (5). Moreover, all the expectation operator w.r.t. the “constrained” pdf of the secondorder modular variate in (3) will be indicated as .
Iii The measurement model
Suppose to have a uniformly linear array (ULA) of omnidirectional sensors and narrowband sources characterized by spatial frequencies . For the ULA configuration, the spatial frequency and the (conic) angle of arrival is linked by where is the spacing between the sensor and is the wavelength of the transmitted signal. The adoption of , instead of , as direction parameter allows us to discard the nonlinearity due to the function and the dependence on the “systemdependent” parameters and . For a ULA, the steering vector can be expressed as . Furthermore, by defining as the vector collecting all the source spatial frequencies, we can define the steering matrix as the matrix whose th column is the steering vector related to the th spatial frequency. In the rest of the paper, we assume to have a set of zeromean, independent and identically CESdistributed (i.i.d.) data snapshots , s. t.:
(6) 
where is the true^{1}^{1}1Note that we use the subscript to distinguish between the true density generator and a generic function in . This notation will be adopted for any other vector, matrix or function in the paper., but generally unknown, density generator. The scatter matrix is:
(7) 
where , is the source covariance matrix, is the identity matrix of dimension and is the “noise” power. The true, but again generally unknown, parameter vector , is defined as:
(8) 
and the vector is the dimensional real vector such that:
(9) 
where the operator selects all the entries strictly below the main diagonal of taken in the same columnwise order as the ordinary operator [11, Sec. 2.4] while is a column vector collecting the diagonal elements of . Remark 1: It is immediate to verify that, when the density generator is , i.e. when the snapshot is Gaussiandistributed, the signal model in (6) is the classical random signal model used e.g. in [3] and [4]:
(10) 
where is the, zero mean, circular Gaussian signal random vector whose covariance matrix is and is the white Gaussian measurement noise. This observation motivates the characterization of the parameter in the general model (6) as noise power.
Iiia The Stochastic CRB
As largely discussed in the array processing literature, we are generally interested in the estimation of , that is the vector of the spatial frequencies, while the other two terms in the parameter vector in (8), i.e. the signal covariance and the noise power have to be considered as nuisance parameters. A CRB for the estimation of in the presence of the (finitedimensional) nuisance parameter vectors and has been discussed in [4] under the Gaussian model assumption discussed in Remark 1. This bound, called the Stochastic CRB (SCRB), is given by:
(11) 
(12) 
where is the Hadamard product, where and
(13) 
IiiB The Semiparametric Stochastic CRB
Following the theoretical results obtained in our recent work [9], we can take these two steps further:

We drop the Gaussianity assumption in favour of the more general CES model in (6),

By relying on the semiparametric framework^{2}^{2}2The reader that is not familiar with the semiparametric theory may have a look at the books [12] and [13] and to the wide statistical literature available on this topic. Moreover, we may suggest the reader to look into our recent works [14, 15, 9] where the semiparmaetric nature of the CES distributions has been analysed., we consider as additional nuisace parameter the density generator itself.
Roughly speaking, the semiparametric framework addressed in [14, 15, 9] allow us to derive a lower bound on the performance of any estimator of the vector of the spatial frequencies when the signal covariance , the noise power and even the density generator are unknown. The only assumption used to derive this bound, that we call Semiparametric SCRB (SSCRB), is that the data snapshots are CESdistributed as in (6). As proved in [9], the SSCRB can be expressed as:
(14) 
where the matrix is the one given in (12), is the 2ndorder modular variate defined in (3) and where is the actual density generator of the data snapshots.
Iv MUSIC with robust scatter matrix estimators
After having introduced the measurement model and the related SSCRB, we now present the MUSIC estimation algorithm whose efficiency, w.r.t. the SSCRB, will be investigated by simulations in Section VI. The MUSIC estimator of the vector of spatial frequencies is given by (see e.g. [16]):
(15) 
where is the steering vector and are the eigenvectors corresponding to the smallest eigenvalues of the estimated data covariance matrix . It is worth underling that the MUSIC algorithm does not require the apriori knowledge of the functional form of the actual density generator , that is generally unknown, so it can be applied in the considered semiparametric framework. Let us now focus our attention on the estimation of the covariance matrix . Also in this case, we have to rely on estimators that do not make use of apriori information on . As a consequence, Maximum Likelihood estimator is not an option. Here, we list five “semiparametric” estimators that we are going to take into account in the efficiency study of the MUSIC algorithm.
Iva The Sample Covariance Matrix (SCM)
The wellknown SCM estimate of is given by:
(16) 
is the ML estimator when the data are Gaussian, but its performance drastically decreases in heavytailed scenarios.
IvB The Normalized (or Sign) and the Kendall’s Tau SCM
IvC Tyler’s and Huber’s estimators
The Tyler’s and Huber’s estimates are the convergence points of the following iterative algorithm:
(20) 
where the starting point is . The weight function for Tyler’s estimator is defined as (see e.g. [7]):
(21) 
whereas the one for Huber’s estimator is given by:
(22) 
where is a tuning parameter and indicates the distribution of a chisquared random variable with degrees of freedom. The parameter is usually chosen as [7]. Note that Tyler’s estimator is the minimax robust estimator of the scatter matrix for CESdistributed data, while Huber’s one represents a compromise between the robustness of Tyler’s estimator (that can be obtained for ) and the efficiency at Gaussianity of the SCM () [7].
V The IAAAPES algorithm
The Iterative Adaptive Approach for Amplitude and Phase EStimation (IAAAPES) is a least squaresbased algorithm that, as the MUSIC algorithm, does not exploit any information on the data distribution but, unlike MUSIC, does not rely of the estimation of the snapshot covariance matrix [20]. Here, a short description of the IAAAPES method is provided. For additional details and and discussions, we refer the reader to [20]. Let be a grid of possible spatial frequencies and let be a diagonal matrix whose diagonal elements are the powers of the potential sources with spatial frequencies . Following [20], we introduce the matrix
(23) 
where and is a matrix whose columns are the steering vectors for each spatial frequency in the grid . The IAAAPES cost function is defined as:
(24) 
and, by minimizing w.r.t. the signal parameter , we get:
(25) 
Finally, the sources spatial frequencies can be identified as the elements of whose indices characterize the smallest diagonal elements of:
(26) 
Note that, since to implement (25), we need an estimate of (see the definition of the matrix in (23)), the estimation of has to be implemented in an iterative algorithm as detailed in [20].
Vi Numerical results
This section is dedicated to the numerical assessment of the efficiency of MUSIC and IAAAPES algorithms w.r.t. the SSCRB in (14). In the following simulations, we assume to have two sources at spatial frequencies and . The noise power while the source covariance matrix is:
(27) 
with . In Figures 1 and 2 where the efficiency is assessed as function of the SignaltoNoise ratio (SNR), and are chosen as and , while in Figures 3 and 4, where the efficiency is assessed as function of the nonGaussianity of the collected data, and are chosen according to dB and dB. In all our simulations, the number of snapshots is , and the number of Monte Carlo runs in . The tuning parameter for the Huber’s estimator is . The IAAAPES algorithm has been implemented with a maximum number of iterations equal to 30. As Mean Square Error indices, we use:
(28) 
where is the Frobenius norm and . As lower bound, we plot the following index:
(29) 
The efficiency study has been conducted for two CES distributions: the complex  and the Generalized Gaussian (GG) distributions. A brief description of these two distributions along with the relevant calculation needed to obtain a closed form for SSCRB in (14) is given below.
Via distributed data
The pdf related to the complex distribution is [9]:
(30) 
and then From (4), we have that:
(31) 
As discussed in Sec. II, we have to constrain the density generator of the distribution to satisfy the constraint in (5). It is immediate to verify that the constraint is satisfied by choosing . Note that, for small values of the shape parameter the distribution have tails heavier that the Normal one, while the distributed data tends to be Gaussian. Using the integral in [21, pp. 315, n. 3.194 (3)], we get:
(32) 
ViB GGdistributed data
The pdf related to the GG distribution is [7, Sec. IV.B]:
(33) 
and then From (4), we have that:
(34) 
The GG distribution could have heavier tails () and lighter tails () as compared to the Normal one (). As discussed in [7, Sec. IV.B], in order to satisfy the constraint is (5), the scale parameter as to be chosen as . Using the integral in [21, pp. 370, n. 3.478 (1)], we get:
(35) 
In Figures 1 and 2 the MSE of the MUSIC and IAAAPES and the related SSCRB are reported as function of the SNR, while in Figures 3 and 4 the efficiency w.r.t the SSCRB of the two DOA estimation algorithms is investigated as function of the nonGaussianity of the collected data. Here, some observations:

In the presence of heavy tailed data, the MUSICTyler and the MUSICHuber algorithms present the best DOA estimation performance in both low and high SNR regimes. On the other hand, as expected, the MUSIC algorithm with the SCM performs poorly in nonGaussian scenarios.

None of the considered DOA estimation algorithms is efficient w.r.t. the SSCRB. However, for a sufficiently large SNR value, MUSICTyler and MUSICHuber algorithms are almost efficient.
Vii Concluding remarks
In this paper, the statistical efficiency of the MUSIC (with different scatter matrix estimates) and the IAAAPES estimators has been assessed in the presence of CESdistributed data whose density generator is unknown and has been considered as an infinitedimensional, nuisance parameter. The SSCRB is the proper bound to be calculated to assess the efficiency of any estimator in such a scenario. Numerical results have shown that the best performance are achieved by exploiting the MUSIC algorithm together with the Tyler’s or Huber’s estimate of the scatter matrix. However, none of the considered estimators is an efficient one w.r.t. the SSCRB. This open problem, together with the experimental validation of the measurement model adopted in this paper, will be addressed in future works.
References
 [1] H. Krim and M. Viberg, “Two decades of array signal processing research: the parametric approach,” IEEE Signal Processing Magazine, vol. 13, no. 4, pp. 67–94, July 1996.
 [2] H. Van Trees, Optimum Array Processing Detection, Estimation and Modulation Theory, Part IV. Wiley, New York, 2002.
 [3] P. Stoica and A. Nehorai, “Performance study of conditional and unconditional directionofarrival estimation,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 38, no. 10, pp. 1783–1795, Oct 1990.
 [4] P. Stoica, E. G. Larsson, and A. B. Gershman, “The stochastic CRB for array processing: a textbook derivation,” IEEE Signal Processing Letters, vol. 8, no. 5, pp. 148–150, May 2001.
 [5] J. P. Delmas, “Performance bounds and statistical analysis of DOA estimation,” in Academic Press Library in Signal Processing: Volume 3, ser. Academic Press Library in Signal Processing, A. M. Zoubir, M. Viberg, R. Chellappa, and S. Theodoridis, Eds. Elsevier, 2014, vol. 3, pp. 719 – 764.
 [6] E. Ollila and V. Koivunen, “Influence function and asymptotic efficiency of scatter matrix based array processors: Case MVDR beamformer,” IEEE Transactions on Signal Processing, vol. 57, no. 1, pp. 247–259, Jan 2009.
 [7] E. Ollila, D. E. Tyler, V. Koivunen, and H. V. Poor, “Complex elliptically symmetric distributions: Survey, new results and applications,” IEEE Transactions on Signal Processing, vol. 60, no. 11, pp. 5597–5625, 2012.
 [8] K. J. Sangston, F. Gini, and M. S. Greco, “Coherent radar target detection in heavytailed compoundgaussian clutter,” IEEE Trans. on Aerospace and Electronic Systems, vol. 48, no. 1, pp. 64–77, Jan 2012.
 [9] S. Fortunati, F. Gini, M. S. Greco, A. M. Zoubir, and M. Rangaswamy, “Semiparametric CRB and SlepianBangs formulas for Complex Elliptically Symmetric distributions,” submitted to IEEE Transactions on Signal Processing, 2019. [Online]. Available: http://arxiv.org/abs/1902.09541
 [10] A. Mennad, S. Fortunati, M. N. E. Korso, A. Younsi, A. M. Zoubir, and A. Renaux, “SlepianBangstype formulas and the related Misspecified CramérRao bounds for complex elliptically symmetric distributions,” Signal Processing, vol. 142, pp. 320 – 329, 2018.
 [11] A. Hjørungnes, ComplexValued Matrix Derivatives. Cambridge University Press, 2011.
 [12] P. Bickel, C. Klaassen, Y. Ritov, and J. Wellner, Efficient and Adaptive Estimation for Semiparametric Models. Johns Hopkins University Press, 1993.
 [13] A. Tsiatis, Semiparametric Theory and Missing Data. Springer series in statistics, 2006.
 [14] S. Fortunati, F. Gini, M. Greco, A. M. Zoubir, and M. Rangaswamy, “A fresh look at the semiparametric CramérRao bound,” in EUSIPCO, Sep. 2018, pp. 261–265.
 [15] S. Fortunati, F. Gini, M. S. Greco, A. M. Zoubir, and M. Rangaswamy, “Semiparametric inference and lower bounds for real elliptically symmetric distributions,” IEEE Transactions on Signal Processing, vol. 67, no. 1, pp. 164–177, Jan 2019.
 [16] P. Stoica and A. Nehorai, “Music, maximum likelihood, and CramérRao bound,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 5, pp. 720–741, May 1989.
 [17] S. Visuri, H. Oja, and V. Koivunen, “Subspacebased directionofarrival estimation using nonparametric statistics,” IEEE Transactions on Signal Processing, vol. 49, no. 9, pp. 2060–2073, Sep. 2001.
 [18] F. Gini and M. Greco, “Covariance matrix estimation for cfar detection in correlated heavy tailed clutter,” Signal Processing, vol. 82, no. 12, pp. 1847 – 1859, 2002.
 [19] S. Bausson, F. Pascal, P. Forster, J. Ovarlez, and P. Larzabal, “First and secondorder moments of the normalized sample covariance matrix of spherically invariant random vectors,” IEEE Signal Processing Letters, vol. 14, no. 6, pp. 425–428, June 2007.
 [20] T. Yardibi, J. Li, P. Stoica, M. Xue, and A. B. Baggeroer, “Source localization and sensing: A nonparametric iterative adaptive approach based on weighted least squares,” IEEE Transactions on Aerospace and Electronic Systems, vol. 46, no. 1, pp. 425–443, Jan 2010.
 [21] I. S. Gradshteyn and M. Ryzhik, Tables of Integrals, Series, and Products(7th edition). Academic Press, Orlando, Florida, 2007.