Statistical Analysis of the Performance of MDL Enumeration for Multiple-Missed Detection in Array Processing

Abstract

An accurate performance analysis on the MDL criterion for source enumeration in array processing is presented in this paper. The enumeration results of MDL can be predicted precisely by the proposed procedure via the statistical analysis of the sample eigenvalues, whose distributive properties are investigated with the consideration of their interactions. A novel approach is also developed for the performance evaluation when the source number is underestimated by a number greater than one, which is denoted as “multiple-missed detection”, and the probability of a specific underestimated source number can be estimated by ratio distribution analysis. Simulation results are included to demonstrate the superiority of the presented method over available results and confirm the ability of the proposed approach to perform multiple-missed detection analysis.

1. Introduction

Source enumeration is a critical step in array signal processing and widely used in many scenarios [1]. The accuracy and the tendency of enumeration will significantly affect the performance of succeeding algorithms, such as direction-of-arrival (DOA) estimation [2] or blind source separation [3]. Minimum description length (MDL) criterion derived by Rissanen [4], or its equivalent criterion derived by Schwarz under the name of Bayesian information criterion (BIC) [5], is one of the most commonly used enumeration methods for its low complexity and asymptotic consistency which ensures a correct estimation as the sample size tends to infinity [6]. A lot of algorithms have been proposed to improve the MDL criterion for a performance promotion, low computation complexity or robustness in various environments [7,8,9,10,11,12,13,14,15]. Dayan and Rausley have presented a norm-based improved MDL (iMDL) algorithm in [7] by nonlinear rescaling of the sample eigenvalues and the corresponding normalized indexes. By using the training sequence of the desired signal, a minimum mean square error (MMSE) based MDL method has been developed by Huang in [8] to get a more accurate estimation of the source number. Huang and So have also employed the linear shrinkage estimation of noise subspace covariance matrix instead of sample covariance matrix in MDL criterion in [9] to achieve a more reliable detection in severe environments where the number of snapshots is comparable or even smaller than the number of sensors. For the cases at a small sample size, the probability density function of the sample eigenvalues has been taken into consideration in MDL as an essential supplement in [10,11]. To handle the coherent signals contaminated by colored noise, Zhen and Si [12] have whitened the sample eigenvalues to eliminate the inequality of eigenvalues caused by colored noise. Fishler and Poor [13] have proposed a robust-MDL (RMDL) method with proven consistency for source enumeration under non-uniform noise situations while Huang et al. [14,15] have improved MDL by introducing a multi-stage Wiener filter by using the filtered component variances or MMSE rather than the sample eigenvalues, which can offer computational simplicity and robustness to non-uniform noise.

Most of the above methods will give an accurate estimation of source number under the assumption of an infinite sample size. However, only a limited number of observations is available in practice particularly in the applications with latency requirements. Thus, the performance analysis of MDL at a finite sample size is of great practical value.

The statistical performance of MDL criterion for source enumeration has been extensively analyzed in [16,17,18,19,20,21,22,23,24,25,26,27]. Since the method is eigendecomposition-based, the statistical property of the sample eigenvalues has been investigated considerably. The distributions of the sample eigenvalues which were derived from the multivariate statistical theory [28], have been used by the authors in [16,17,18]. However, the performance estimations are found biased when the sample size is not sufficiently large. Recently, the random matrix theory [29] approach has been proposed to solve the enumeration problem in array processing in [30,31,32,33]. Asymptotic distributions of the sample eigenvalues have been given to rectify the bias by taking in the influence on signal eigenvalues from the noise subspace under the assumption of large dimension. For non-circular or non-Gaussian cases, statistical analysis has been made in [22,23] by taking the fourth-order statistics of the signals into consideration. In [34], the interactions between signal eigenvalues are considered by Lawley on the distribution analysis of the sample eigenvalues. The authors in [24,25,26,27] are able to predict the probability of underestimation in close accordance at a moderate sample size by a combination of different theories.

Although many analyses on the performance of MDL are available, few discuss the cases that the source number is underestimated by a number greater than one, which is denoted as “multiple-missed detection” here. Since the enumeration performance of MDL is signal-to-noise ratio (SNR) dependent, the source numbers estimated under varying noise levels may be different. This inconsistency in enumeration can be attributed to either the variation of actual source number or the disturbance of noise. By the statistical analysis on multiple-missed detection, the probability of a specific enumeration result can be estimated as an important reference for the attribution of the enumeration discrepancy.

In this paper, we propose a new procedure for the performance evaluation of MDL, which can predict the estimation results of MDL precisely at a finite sample size by considering the interactions between signal eigenvalues. A novel approach is also developed for the multiple-missed detection analysis. Thus the deterioration of enumeration performance with the degradation of SNR can be estimated.

The remainder of this paper is organized as follows: the problem formulations are given in Section 2, including the array signal model and the theoretical analysis of underestimation. Section 3 introduces the statistical analysis of underestimation by discussing the distributive property of the sample eigenvalues. Simulation results that illustrate the superior performance of the proposed method and the performance for multiple-missed detection are presented in Section 4. Finally, conclusions are drawn in Section 5.

2. Problem Formulation

2.1. Array Signal Model

Consider q narrowband far-field and incoherent sources impinging on a sensor array of p elements (p > q). The observed signals can be modeled as a superposition of source signals corrupted by additive circular Gaussian noise, which can be written as:

$$x=As+n$$

(1)

where A is the p × q array steering matrix composed of q linearly independent column vectors, s is the q-dimensional source signal vector with nonsingular covariance matrix R_S = E [ss^H] where (.)^H stands for conjugate transpose, n is the source-independent noise vector with zero mean and covariance σ²I where I is the p × p identity matrix. Signals and noises are assumed to be i.i.d. and complex circular Gaussian distributed. The p × p population covariance matrix R is calculated as:

$$R=E[x{x}^H]=A{R}_S{A}^H+{\mathrm{\sigma }}^{2}I$$

(2)

whose population eigenvalues in descending order are given by:

$${\mathrm{\lambda }}_1\ge {\mathrm{\lambda }}_2\ge \mathrm{...}\ge {\mathrm{\lambda }}_q>{\mathrm{\lambda }}_{q+1}=\mathrm{...}={\mathrm{\lambda }}_p={\mathrm{\sigma }}^2$$

(3)

The first q eigenvalues of R are contributed by both the source signals and the noise, which are called the signal eigenvalues. The last p – q eigenvalues are contributed by noise only, which are called the noise eigenvalues. The population covariance matrix can be estimated using the sample covariance matrix $\widehat{R}$:

$$\widehat{R}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{x}_i{x}_i^H}$$

(4)

where x₁, ... , x_N are the independent and identically distributed snapshots of x. The corresponding sample eigenvalues of $\widehat{R}$ in descending order are given by:

$$l_1>l_2>\mathrm{...}>l_q>l_{q+1}>\mathrm{...}>l_p$$

(5)

2.2. Source Enumeration of MDL

Assuming that all the observations are i.i.d. complex circular Gaussian random vectors with zero mean, the MDL estimator formulation is given by:

$$MDL(k)=-L(k)+y(k)=N(p-k)\ln \frac{A_k}{G_k}+\frac{1}{2}\left[k(2p-k)+1\right]\ln N$$

(6)

where L(k) is the log-likelihood term, y(k) is the penalty term:

$$A_k\triangleq \frac{1}{p-k}{\displaystyle \sum _{i=k+1}^p{l}_i}$$

(7)

and:

$$G_k\triangleq {\displaystyle \prod _{i=k+1}^p{l}_i^{1/(p-k)}}$$

(8)

represent the arithmetic and geometric means of the last p – k eigenvalues respectively. The estimated source number is denoted by $\widehat{q}$ which can be derived as follows:

$$\widehat{q}=\underset{k}{\arg \min }MDL(k),k=0,\mathrm{...},p-1$$

(9)

Let H_q denote the hypothesis that the true number of sources is q. The probability of incorrect estimation P_e is defined as:

$$P_e=P\left(\widehat{q}\ne q|H_q\right)=P_m+P_f$$

(10)

where the probability of missed detection P_m is defined as:

$$P_m=P\left(\widehat{q}<q|H_q\right)$$

(11)

and the probability of false alarm P_f is defined as:

$$P_f=P\left(\widehat{q}>q|H_q\right)$$

(12)

Ding and Kay [35] have proven that MDL is inconsistent at high SNR with a finite sample size. However, for the relatively large penalty term, MDL has a trend of underestimation of the source number at low SNR. For example, in the setting of p = 5, q = 2, N = 50 and SNR = 3 dB, the probability of false alarm P_f is 0.0013 while the probability of underestimation P_m is 0.5029 based on a 10,000 trial Monte Carlo simulation. Zhang et al., showed in [18] that for a moderate number of snapshots, the probability of false alarm using the MDL criterion is approaching zero. So the probability of incorrect estimation P_e is dominated by the probability of underestimation P_m, which can be expressed as:

$$P_e\approx P_m=P(\widehat{q}<q|H_q)\approx P\left(MDL(q-1)<MDL(q)|H_q\right)$$

(13)

Using Equation (6) in Equation (13) we obtain:

$$P_m\approx P\left(\frac{A_{q-1}^{p-q+1}}{l_q{A}_q^{p-q}}<\exp \left(\frac{(2p-2q+1)\ln N}{2N}\right)\right)$$

(14)

According to the definition of A_k in Equation (7), we can rewrite A_q−1 as:

$$A_{q-1}=\frac{{\displaystyle \sum _{i=q}^p{l}_i}}{p-q+1}=\frac{l_q+{\displaystyle \sum _{i=q+1}^p{l}_i}}{p-q+1}\text{=}\frac{l_q+(p-q)A_q}{p-q+1}$$

(15)

In order to simplify Equation (14), we define:

$${\mathrm{\rho }}_i=l_i/A_i,{\mathrm{\rho }}_i>1$$

(16)

when i = q, we will have:

$${\mathrm{\rho }}_q=l_q/A_q,{\mathrm{\rho }}_q>1$$

(17)

Substituting Equations (15) and (17) into Equation (14), we can get:

$$P_m\approx P\left(f_1({\mathrm{\rho }}_q)<C_1\right)$$

(18)

where:

$$\mathrm{w}{f}_1(x)\triangleq \left(x+p-q\right)x^{-\frac{1}{p-q+1}}$$

(19)

and:

$$C_1=(p-q+1)\exp \left(\frac{(2p-2q+1)\ln N}{2N(p-q+1)}\right)$$

(20)

The function f₁(x) is a monotonically increasing function in the region of x > 1, therefore we can transform Equation (17) into a simpler form as:

$$P_m\approx P\left({\mathrm{\rho }}_q<{\mathrm{\rho }}_{C1}\right)$$

(21)

in which:

$$f_1({\mathrm{\rho }}_{C1})=C_1$$

(22)

Since Equation (21) cannot be solved analytically, we can use the Newton-Raphson method to find a very accurate solution of ${\mathrm{\rho }}_{C1}$ numerically from the initial value derived by binomial expansion:

$${\mathrm{\rho }}_{C1}^0=D_1+\sqrt{D_1^2-1}$$

(23)

where:

$$D_1={\left(\frac{C_1}{p-q+1}\right)}^{p-q+1}$$

(24)

For a more complicated situation, we will discuss the cases of multiple-missed detection as follows. Consider that the true source number q is underestimated by a number greater than or equal to d, we define the corresponding probability as:

$$P_{md}=P\left(\widehat{q}\le q-d|H_q\right),d<q$$

(25)

in order to distinguish the probabilities for different underestimated source numbers. So:

$$P_{md}\approx P\left(MDL(q-d)<MDL(q-d+1)|H_q\right)$$

(26)

By using Equation (6) in Equation (26), we can derive:

$$P_{md}\approx P\left(\frac{A_{q-d}^{p-q+d}}{l_{q-d+1}{A}_{q-d+1}^{p-q+d-1}}<\exp \left(\frac{(2p-2q+2d-1)\ln N}{2N}\right)\right)$$

(27)

Let i = q – d + 1 in Equation (16), we will have:

$${\mathrm{\rho }}_{q-d+1}=l_{q-d+1}/A_{q-d+1},{\mathrm{\rho }}_{q-d+1}>1$$

(28)

and obtain a simpler form of Equation (26):

$$P_{md}\approx P\left(f_d({\mathrm{\rho }}_{q-d+1})<C_d\right)$$

(29)

where:

$$f_d(x)\triangleq \left(x+p-q+d-1\right)x^{-\frac{1}{p-q+d}}$$

(30)

$$C_d=(p-q+d)\exp \left(\frac{(2p-2q+2d-1)\ln N}{2N(p-q+d)}\right)$$

(31)

Since f_d(x) is also a monotonically increasing function in the region of x > 1, Equation (29) can be rewritten as:

$$P_{md}\approx P\left({\mathrm{\rho }}_{q-d+1}<{\mathrm{\rho }}_{Cd}\right)$$

(32)

in which:

$$f_1({\mathrm{\rho }}_{Cd})=C_d$$

(33)

The threshold ${\mathrm{\rho }}_{Cd}$ can be calculated numerically from the initial value:

$${\mathrm{\rho }}_{Cd}^0=D_d+\sqrt{D_d^2-1}$$

(34)

where:

$$D_d={\left(\frac{C_d}{p-q+d}\right)}^{p-q+d}$$

(35)

We can find that P_m is a particular form of P_md when d = 1. So the problem of underestimation probability turns into the statistical performance analysis of ρ_i which will be discussed in the next section. The expectation of the estimated source number can be calculated by:

$$E(\widehat{q})=q(1-P_{m1})+P_{m(q-1)}+{\displaystyle \sum _{d=1}^{q-2}(q-d)P_{md}(1-P_{m(d+1)})},\widehat{q}\in [1,q]$$

(36)

which would be an effective indicator of the extent of underestimation.

3. Performance Analysis of Multiple-Missed Detection

According to Equation (16), the statistics of ρ_i are determined by the distributions of l_i and A_i. Many researches have been done on the statistical properties of the signal sample eigenvalue l_i and the arithmetic mean of the noise sample eigenvalues A_q, and can be mainly divided into multivariate statistical theory [28], random matrix theory [29] and Lawley’s theory [34]. The multivariate statistical theory is derived based on large sample asymptotics and requires a large value of sample size N. Random matrix theory has been proposed to investigate the spectral properties of random matrices with the assumption of high-dimension and large sample asymptotic regime. Nadakuditi and Edelman have concluded in [32] that for a signal-free sample covariance matrix formed from a p × N matrix of observations with i.i.d. Gaussian samples of zero mean and identical variance ${\mathrm{\sigma }}^2$, the sample eigenvalues will follow the Marchenko-Pastur distribution and their arithmetic mean will converge to Gaussian distribution asymptotically as p, N →∞ with $p/N\to c\in (0,\infty )$, i.e.,

$$p\left(\frac{1}{p}{\displaystyle \sum _{i=1}^p{l}_i}-{\mathrm{\sigma }}^2\right)\stackrel{D}{\to }N\left(0,\frac{2{\mathrm{\sigma }}^{4}c}{\mathrm{\beta }}\right)$$

(37)

where $\stackrel{D}{\to }$ denotes convergence in distribution and β = 1 or 2 for real or complex values respectively. For a q-signal-bearing case defined in Section 2.1, as p, N → ∞ with $p/N\to c\in (0,\infty )$, if all the signal eigenvalues of the population covariance matrix are larger than the critical value ${\mathrm{\sigma }}^2(1+\sqrt{c})$ where ${\mathrm{\sigma }}^2$ stands for the value of noise eigenvalue, which means no phase transition phenomenon, the distribution of signal sample eigenvalues can be described as following. If the signal eigenvalue ${\mathrm{\lambda }}_i>{\mathrm{\sigma }}^2(1+\sqrt{c})$ has multiplicity of one for i ≤ q and $\sqrt{N}\left|c-p/N\right|\to 0$, the distribution of l_i converges to Gaussian distribution as:

$$\sqrt{N}\left(l_i-{\mathrm{\lambda }}_i\left(1+\frac{{\mathrm{\sigma }}^{2}c}{{\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2}\right)\right)\stackrel{D}{\to }N\left(0,\frac{2}{\mathrm{\beta }}{\mathrm{\lambda }}_i^2\left(1-\frac{c{\mathrm{\sigma }}^4}{{\left({\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2\right)}^2}\right)\right)$$

(38)

The distribution analysis of A_q can be performed as the arithmetic mean of the sample eigenvalues of a signal-free sample covariance matrix formed from a (p – q) × N matrix of observations. Although the random matrix theory is derived under the assumption of large dimension, some simulation results have shown that it may also work well in some low-dimension cases. Lawley has constructed a matrix with the same eigenvalues as the sample covariance matrix by using the sampling errors. By comparing the diagonal elements, the statistics of sample eigenvalues are derived under the assumption of Gaussian distribution. All the three theories assume that l_i and A_q follow Gaussian distribution asymptotically and the expectations and variances are listed in Table 1 for comparison.

Note that the expectations and variances of l_i in random matrix theory and Lawley’s theory have an additional term compared with those in multivariate statistical theory. The random matrix theory has taken the disturbance from noise subspace into account while Lawley’s theory has considered the interactions between signal eigenvalues in addition. When the sample size N is sufficiently large, the additional terms will diminish to zero and the distributions of all the three theories will equal to each other. Similarly, the expectation of A_q in Lawley’s theory is different from the others for including the bias induced by the signal eigenvalues. Thus, we employ the expectation and variance of l_i in Lawley’s theory for an accurate analysis, and additional terms with the order higher than O(N⁻¹) in the expectations or O(N⁻²) for variances are omitted as they are relatively small and decay rapidly with increasing sample size. The phenomenon of phase transition is not considered for simplicity in this paper.

Table 1. Comparison of three theories on the expectations and variances of l_i and A_q.

	Multivariate Statistics	Random Matrix	Lawley
E(li) *1	${\mathrm{\lambda }}_i$	${\mathrm{\lambda }}_i+\frac{c{\mathrm{\lambda }}_i{\mathrm{\sigma }}^2}{{\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2}$ *3	${\mathrm{\lambda }}_i+{\displaystyle \sum _{j\ne i}^p\frac{{\mathrm{\lambda }}_i{\mathrm{\lambda }}_j}{N({\mathrm{\lambda }}_i-{\mathrm{\lambda }}_j)}}$
Var(li) *1	$\frac{{\mathrm{\lambda }}_i^2}{N}$	$\frac{2{\mathrm{\lambda }}_i^2}{\mathrm{\beta }N}\left[1-\frac{c{\mathrm{\sigma }}^4}{{({\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2)}^2}\right]$ *2,*3	$\frac{2{\mathrm{\lambda }}_i^2}{\mathrm{\beta }N}\left[1-\frac{1}{N}{\displaystyle \sum _{j\ne i}^p{\left(\frac{{\mathrm{\lambda }}_j}{{\mathrm{\lambda }}_i-{\mathrm{\lambda }}_j}\right)}^2}\right]$ *2
E(Aq)	${\mathrm{\sigma }}^2$	${\mathrm{\sigma }}^2$	${\mathrm{\sigma }}^2-{\displaystyle \sum _{i=1}^q\frac{{\mathrm{\lambda }}_i{\mathrm{\sigma }}^2}{N({\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2)}}$
Var(Aq)	$\frac{{\mathrm{\sigma }}^4}{N(p-q)}$	$\frac{2{\mathrm{\sigma }}^4}{\mathrm{\beta }N(p-q)}$ *2	/

*¹ i ≤ q for the expectations and variances of l_i. E(.) and Var(.) are the mathematical expectation and variance respectively; *² β = 1 for real-valued signals and β = 2 for complex-valued signals. The original formula of Lawley’s Var(l_i) without β is revised here; *³ c is a positive finite value when p, N →∞, p/N → c.

Table 1. Comparison of three theories on the expectations and variances of l_i and A_q.

	Multivariate Statistics	Random Matrix	Lawley
E(li) *1	${\mathrm{\lambda }}_i$	${\mathrm{\lambda }}_i+\frac{c{\mathrm{\lambda }}_i{\mathrm{\sigma }}^2}{{\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2}$ *3	${\mathrm{\lambda }}_i+{\displaystyle \sum _{j\ne i}^p\frac{{\mathrm{\lambda }}_i{\mathrm{\lambda }}_j}{N({\mathrm{\lambda }}_i-{\mathrm{\lambda }}_j)}}$
Var(li) *1	$\frac{{\mathrm{\lambda }}_i^2}{N}$	$\frac{2{\mathrm{\lambda }}_i^2}{\mathrm{\beta }N}\left[1-\frac{c{\mathrm{\sigma }}^4}{{({\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2)}^2}\right]$ *2,*3	$\frac{2{\mathrm{\lambda }}_i^2}{\mathrm{\beta }N}\left[1-\frac{1}{N}{\displaystyle \sum _{j\ne i}^p{\left(\frac{{\mathrm{\lambda }}_j}{{\mathrm{\lambda }}_i-{\mathrm{\lambda }}_j}\right)}^2}\right]$ *2
E(Aq)	${\mathrm{\sigma }}^2$	${\mathrm{\sigma }}^2$	${\mathrm{\sigma }}^2-{\displaystyle \sum _{i=1}^q\frac{{\mathrm{\lambda }}_i{\mathrm{\sigma }}^2}{N({\mathrm{\lambda }}_i-{\mathrm{\sigma }}^2)}}$
Var(Aq)	$\frac{{\mathrm{\sigma }}^4}{N(p-q)}$	$\frac{2{\mathrm{\sigma }}^4}{\mathrm{\beta }N(p-q)}$ *2	/

*¹ i ≤ q for the expectations and variances of l_i. E(.) and Var(.) are the mathematical expectation and variance respectively; *² β = 1 for real-valued signals and β = 2 for complex-valued signals. The original formula of Lawley’s Var(l_i) without β is revised here; *³ c is a positive finite value when p, N →∞, p/N → c.

The expectation and variance of A_i are given only under the condition of i = q in above theories. Now we will discuss the distributive property of A_i to get general expressions when i < q. According to the definition of A_i in Equation (7), we can rewrite A_i as:

$$A_i=\frac{1}{p-i}\left((p-q)A_q+{\displaystyle \sum _{j=i+1}^q{l}_j}\right)$$

(39)

Noting that:

$$E\left((p-i)A_i+{\displaystyle \sum _{j=1}^i{l}_i}\right)=E(tr(\widehat{R}))=tr(R)={\displaystyle \sum _{j=1}^p{l}_j}$$

(40)

we can obtain:

$$E(A_i)=\frac{1}{p-i}\left\{{\displaystyle \sum _{j=i+1}^q\left[{\mathrm{\lambda }}_j+{\displaystyle \sum _{k=1}^i\frac{{\mathrm{\lambda }}_j{\mathrm{\lambda }}_k}{N({\mathrm{\lambda }}_j-{\mathrm{\lambda }}_k)}}\right]}-{\displaystyle \sum _{j=1}^i\frac{{\mathrm{\lambda }}_j{\mathrm{\sigma }}^2(p-q)}{N({\mathrm{\lambda }}_j-{\mathrm{\sigma }}^2)}+(p-q){\mathrm{\sigma }}^2}\right\},i\le q$$

(41)

which equals to the expectation of A_q in Lawley’s theory as i = q. When i < q, A_i will contain signal sample eigenvalues. We assume that the covariance between signal and noise sample eigenvalues can be ignored. The covariance between signal sample eigenvalues is given in [34] as:

$$Cov(l_i,l_j)=\frac{2{\mathrm{\lambda }}_i^2{\mathrm{\lambda }}_j^2}{\mathrm{\beta }{N}^2{({\mathrm{\lambda }}_i-{\mathrm{\lambda }}_j)}^2},i\ne j，i,j\le q$$

(42)

According to the properties of variance and complex circular Gaussian assumption, A_i is also asymptotically Gaussian distributed and the variance can be derived as:

$$Var(A_i)=Var\left(\frac{p-q}{p-i}{A}_q+\frac{{\displaystyle \sum _{j=i+1}^q{l}_j}}{p-i}\right)=\frac{{(p-q)}^{2}Var(A_q)}{{(p-i)}^2}+{\displaystyle \sum _{j=i+1}^q\frac{Var(l_j)}{{(p-i)}^2}}+2{\displaystyle \sum _{j=i+1}^q{\displaystyle \sum _{j<k\le q}\frac{Cov(l_j,l_k)}{{(p-i)}^2}}}$$

(43)

Note that Equation (43) equals to the variance of A_q when i = q, so Equation (43) is the general expression of the variance of A_i for i ≤ q. Since the distributions of l_i and A_i are asymptotically Gaussian, the distribution of ρ_i defined in Equation (16) would follow a ratio distribution of two Gaussian random variables. The probability density function of the ratio of two correlated Gaussian random variables has been derived by Hinkley in [36]. The correlation coefficient r_i between l_i and A_i can be derived as:

$$r_i=\frac{Cov(l_i,A_i)}{\sqrt{Var(l_i)Var(A_i)}}=\frac{1}{(p-i)\sqrt{Var(l_i)Var(A_i)}}{\displaystyle \sum _{j=i+1}^{q}Cov(l_i,l_j)}，i<q$$

(44)

So the probability density function f(ρ_i) of ρ_i can be obtained by using Equation (1) of [36]. Furthermore, the probability of multiple-missed detection P_md can be calculated as:

$$P_{md}=F({\mathrm{\rho }}_{Cd})={\displaystyle {\int }_{-\infty }^{{\mathrm{\rho }}_{Cd}}f({\mathrm{\rho }}_i)}d{\mathrm{\rho }}_i$$

(45)

where:

$$F(x)={\displaystyle {\int }_{-\infty }^{x}f(t)}dt$$

(46)

is the cumulative distribution function of the Gaussian ratio distribution.

4. Simulation Setup and Numerical Results

In the numerical simulations, a uniform linear array with an inter-sensor spacing of half-wavelength is employed. The observed signals are assumed to be uncorrelated complex circular Gaussian source signals contaminated by additive complex circular white Gaussian noise and some results may be invalid for arbitrary complex signals. The numbers of samples, sensors, true sources and underestimated sources are denoted by N, p, q and d respectively. The DOAs are denoted by the vector θ and SNR is short for the signal-to-noise ratio. The probability of an underestimated source number greater than or equal to d is denoted by P_md and the expectation of the estimated source number is denoted by $E(\widehat{q})$. All the simulation results are obtained based on 10,000 Monte Carlo trials.

4.1. Evaluation of the Proposed Method for Underestimation Analysis

The methods presented in [25,26,27] are used for comparison which based on the statistical analysis of sample eigenvalues as well. Haddadi et al. [25] use the expectations in Lawley’s theory and the variance in multivariate statistical theory with the neglect of the variances of noise sample eigenvalues while Huang et al. [26] employ the variance in random matrix theory. Lu and Zoubir [27] have incorporated the expectations in Lawley’s theory and the variances in random matrix theory to predict the estimation results of MDL precisely. All the four procedures use the same expectations of l_i and A_q whereas different variances are selected. A comparison among different procedures is shown in Table 2.

Four experimental settings are listed as following by varying the SNR or the sample size N:

• Setting 1 (see Figure 1): N = 50, p = 20, q = 3, θ = {−5°, 0°, 5°}, SNR = [−10,−4] dB.
• Setting 2 (see Figure 2): N = 1000, p = 20, q = 3, θ = {−5°, 0°, 5°}, SNR = [−15.5,−12.5] dB.
• Setting 3 (see Figure 3): N = [100,500], p = 30, q = 3, θ = {−5°, 5°, 10°}, SNR = −12 dB.
• Setting 4 (see Figure 4): N = [1000,4000], p = 30, q = 3, θ = {−5°, 5°, 10°}, SNR = −17 dB.

Table 2. Theoretical comparison of the four procedures.

	Haddadi et al.	Huang et al.	Lu &Zoubir	Ours
E(lq)	Lawley	Lawley	Lawley	Lawley
Var(lq)	Multivariate Statistics	Random Matrix	Random Matrix	Lawley
E(Aq)	Lawley	Lawley	Lawley	Lawley
Var(Aq)	- *	-	Random Matrix	Random Matrix

* Means this term has been neglected.

Table 2. Theoretical comparison of the four procedures.

	Haddadi et al.	Huang et al.	Lu &Zoubir	Ours
E(lq)	Lawley	Lawley	Lawley	Lawley
Var(lq)	Multivariate Statistics	Random Matrix	Random Matrix	Lawley
E(Aq)	Lawley	Lawley	Lawley	Lawley
Var(Aq)	- *	-	Random Matrix	Random Matrix

* Means this term has been neglected.

As shown in Figure 1, the proposed method shows the best agreement with simulation results and outperforms the others. It is worth noting that the curves nearly intersect at the same point of P_m = 0.5 with different shapes which may correspond to the expectations and variances, respectively. The methods in [25,27] have very similar performance since the former uses a larger variance of signal sample eigenvalue while the latter considered the variance of noise sample eigenvalues as a counteract. An inconspicuous difference is found for the method in [26] due to the consideration of only the interactions between signal and noise subspaces.

The simulation results at a large sample size are presented in Figure 2. All the four methods match the simulation results pretty well and the superiority of the proposed method and the method in [26] can be confirmed in the details of Figure 2. The accurate prediction by the method in [26] may be attributed to the reason that the ignorance of the variance of noise sample eigenvalues may compensate the interactions between signal eigenvalues. Similar results are presented in Figure 3 and Figure 4 by varying the sample size N instead of SNR and the proposed method outperforms the other methods.

The outstanding performance of the proposed method can be attributed to the fact that the interactions between signal eigenvalues have been taken into consideration in estimating the variance of signal sample eigenvalue. The performance of the methods in [25,26,27] is dependent of the sample size in the simulation settings, since they ignore the interactions between signal eigenvalues which are sensitive to the sample size. In the cases when the sample size is sufficiently large, all the methods are capable to yield satisfactory results. However, MDL is widely used in practical applications where only a limited number of samples is available. Thus, the proposed method is of more practical value for its accuracy in such cases.

Figure 1. The probability of underestimation P_m1 versus SNR with N = 50, p = 20, q = 3, θ = {−5°, 0°, 5°}.

Figure 1. The probability of underestimation P_m1 versus SNR with N = 50, p = 20, q = 3, θ = {−5°, 0°, 5°}.

Figure 2. The probability of underestimation P_m1 versus SNR and details with N = 1000, p = 20, q = 3, θ = {−5°, 0°, 5°}.

Figure 2. The probability of underestimation P_m1 versus SNR and details with N = 1000, p = 20, q = 3, θ = {−5°, 0°, 5°}.

Figure 3. The probability of underestimation P_m1 versus N with p = 30, q = 3, θ = {−5°, 5°, 10°}, SNR = −12 dB.

Figure 3. The probability of underestimation P_m1 versus N with p = 30, q = 3, θ = {−5°, 5°, 10°}, SNR = −12 dB.

Figure 4. The probability of underestimation P_m1 versus N and details with p = 30, q = 3, θ = {−5°, 5°, 10°}, SNR = −17 dB.

Figure 4. The probability of underestimation P_m1 versus N and details with p = 30, q = 3, θ = {−5°, 5°, 10°}, SNR = −17 dB.

4.2. Evaluation of the Analysis on Multiple-Missed Detection

To evaluate the performance of the proposed method for multiple-missed detection, the algorithms in [25,26,27] have been adapted by the proposed approach as reference methods. Experimental settings are listed as following:

• Setting 5 (see Figure 5): N = 150, p = 30, q = 5, θ = {−7°, −5°, 0°, 3°, 11°}, SNR = [−14,−10] dB, d = 4.
• Setting 6 (see Figure 6): N = 400, p = 30, q = 5, θ = {−7°, −5°, 0°, 3°, 11°}, SNR = [−16.5,−13] dB, d = 4.

Figure 5. The probability of underestimation P_m4 versus SNR with N = 150, p = 30, q = 5, θ = {−7°, −5°, 0°, 3°, 11°}, d = 4.

Figure 5. The probability of underestimation P_m4 versus SNR with N = 150, p = 30, q = 5, θ = {−7°, −5°, 0°, 3°, 11°}, d = 4.

Figure 6. The probability of underestimation P_m4 versus SNR with N = 400, p = 30, q = 5, θ = {−7°, −5°, 0°, 3°, 11°}, d = 4.

Figure 6. The probability of underestimation P_m4 versus SNR with N = 400, p = 30, q = 5, θ = {−7°, −5°, 0°, 3°, 11°}, d = 4.

The simulation results are shown in Figure 5 and Figure 6 at different sample sizes. As indicated in Figure 5, the proposed method is more accurate in the region of P_m < 0.5 while the other three methods perform better in the region of P_m > 0.5. When the sample size is increased to 400, the discrepancy among the methods become smaller and the superiority of the proposed method is reconfirmed in the entire region. We attribute this phenomenon to the omission of the higher-order terms, since a bias can be found between the simulation results and all the predictions in Figure 5, while the bias is perfectly rectified at a slightly larger sample size in Figure 6. The accumulated error by the neglect of higher-orders terms would affect the performance analysis for multiple-missed detection and a moderate sample size may be required for an accurate estimation.

To assess the performance deterioration of enumeration with the degradation of SNR, we use the expectation of the estimated source number as the indicator which is calculated by Equation (36). Experimental settings are listed as follows:

• Setting 7 (see Figure 7): N = 150, p = 20, q = 5, θ = {−10°, −6°, 0°, 6°, 8°}, SNR = [−15,0] dB.
• Setting 8 (see Figure 8): N = 300, p = 20, q = 5, θ = {−11°, −7°, 0°, 2°, 10°}, SNR = [−16,−2] dB.

Figure 7 and Figure 8 show that all the four methods are able to perform the analysis for multiple-missed detection precisely and the capability of the proposed approach for multiple-missed detection is verified. The mean absolute error (MAE) between the prediction and simulation results is selected to assess the performance of different methods quantitatively. The MAEs are listed in Table 3 and the superiority of our method can be confirmed.

Table 3. MAE of the methods in predicting the enumeration results.

MAE	Haddadi et al.	Huang et al.	Lu & Zoubir	Our Method
Setting 7	0.0143	0.0123	0.0140	0.0110
Setting 8	0.0054	0.0047	0.0055	0.0043

Table 3. MAE of the methods in predicting the enumeration results.

MAE	Haddadi et al.	Huang et al.	Lu & Zoubir	Our Method
Setting 7	0.0143	0.0123	0.0140	0.0110
Setting 8	0.0054	0.0047	0.0055	0.0043

Figure 7. The expectation of the estimated source number $E(\widehat{q})$ versus SNR and details with N = 150, p = 20, q = 5, θ = {−10°, −6°, 0°, 6°, 8°}.

Figure 7. The expectation of the estimated source number $E(\widehat{q})$ versus SNR and details with N = 150, p = 20, q = 5, θ = {−10°, −6°, 0°, 6°, 8°}.

Figure 8. The expectation of the estimated source number $E(\widehat{q})$ versus SNR and details with N = 300, p = 20, q = 5, θ = {−11°, −7°, 0°, 2°, 10°}.

Figure 8. The expectation of the estimated source number $E(\widehat{q})$ versus SNR and details with N = 300, p = 20, q = 5, θ = {−11°, −7°, 0°, 2°, 10°}.

5. Conclusions

This paper presents an accurate performance analysis for the underestimation performance of the MDL source enumeration method at a finite sample size in array processing. Theoretical derivations and statistical analyses have been performed with the consideration of the interactions between signal eigenvalues to obtain an improved estimation of the probability of underestimation. A new approach is also proposed to evaluate the performance of multiple-missed detection cases by ratio distribution analysis and can be employed by the eigenvalue-analyzing methods. Simulation results show the superiority of the proposed analysis, and verify the ability of the proposed approach in evaluating the deterioration of enumeration performance with the degradation of SNR, which may be a valuable reference for practical applications.