THEORETICAL ANALYSIS OF AN IMPROVED COVARIANCE MATRIX ESTIMATOR IN NON-GAUSSIAN NOISE F. Pascal123 , P. Forster2, J-P. Ovarlez1 , P. Larzabal3. 1
ONERA, Chemin de la Huni`ere, F-91761 Palaiseau Cedex, France 2 GEA, 1 Chemin Desvalli`eres, F-92410 Ville d’Avray, France 3 ENS Cachan, 61 Avenue du Pr´esident Wilson, F-94235 Cachan Cedex, France Email :
[email protected],
[email protected],
[email protected],
[email protected] where
ABSTRACT This paper presents a detailed theoretical analysis of a recently introduced covariance matrix estimator, called the Fixed Point Estimate (FPE). It plays a significant role in radar detection applications. This estimate is provided by the Maximum Likelihood Estimation (MLE) theory when the non-Gaussian noise is modelled as a Spherically Invariant Random Process (SIRP). We study in details its properties: existence, uniqueness, unbiasedness, consistency and asymptotic distribution. We propose also an algorithm for its computation and prove the convergence of this numerical procedure. These results will allow to study the performance analysis of the adaptive CFAR radar detectors (GLRT-LQ, BORD,...).
Non-Gaussian noise characterization has gained many interests since experimental radar clutter measurements showed that these data are correctly described by non-Gaussian statistical models. One of the most tractable and elegant non-Gaussian model comes from the so-called Spherically Invariant Random Process (SIRP) theory. A SIRP is the product of a Gaussian random process - called speckle - with a non-negative random variable - called texture. This model leads to many results [1, 2, 3, 4]. The basic problem of detecting a complex signal corrupted by an additive SIRP noise c in a m-dimensional complex vector y allowed to build several Generalized Likelihood Ratio Tests like the GLRT-Linear Quadratic (GLRT-LQ) in [1, 2] or the Bayesian Optimum Radar Detector (BORD) in [3, 4]. Let us recall some SIRP theory results. A noise modelled as a SIRP is a non-homogeneous Gaussian process with random power. More precisely, a SIRP [5] is the product of a positive random variable τ (texture) and a m-dimensional independent complex Gaussian vector x (speckle) with zero mean covariance matrix M = E(xx† ) with normalization tr (M) = m, where † denotes the conjugate transpose operator : √
τx.
The SIRP PDF expression is: Z pm (c) =
+∞ 0
gm (c, τ ) p(τ ) dτ ,
(1)
„ † −1 « c M c 1 exp − . (π τ )m |M| τ
(2)
In many problems, non-Gaussian noise can be characterized by SIRPs but the covariance matrix M is generally not known b is required. Obviously, it has to satisfy the and an estimate M b = m. M-normalization: tr (M) In the literature [6, 7], the Normalized Sample Covariance Matrix Estimate (NSCME) defined as follows is usually used:
b NSCM E M
1. PROBLEM STATEMENT AND BACKGROUND
c=
gm (c, τ ) =
mX = N i=1 N
ci c†i
c†i ci
!
mX = N i=1 N
xi x†i
!
x†i xi
,
(3)
√ where 1 ≤ i ≤ N , ci = τi xi are m-dimensional independent complex SIRP vectors and xi are m-dimensional independent complex Gaussian vectors with zero mean and normalized covariance matrix M. This estimate is remarkably independent of the texture statistics. However, despite of this interesting property, the NSCME (3) suffers the following drawbacks: • it is a biased estimate; • the resulting adaptive Generalized Likelihood Ratio (GLR) is not independent of matrix M characteristics. In the next sections, we propose to introduce and analyze an improved estimator of M: the FPE estimator. b FP 2. THE FIXED POINT ESTIMATOR M Conte and Gini in [8, 9] have shown that the Maximum Likelihood b of M is a solution of the following equation: estimator M ! N X ci c†i m b = . (4) M N i=1 c† M b −1 ci i Existence and uniqueness of the above equation solution have been investigated in a previous paper [10] and we briefly recall here the main results in Theorem 1.
Let function f be defined as: b = f (M)
mX N i=1
!
ci c†i
N
−1
b c†i M
,
(5)
to fill these gaps by establishing the bias, the consistency and the asymptotic distribution of this random matrix: this is the aim of the next section.
ci
b and notice that f can be rewritten as follows, which shows that M is also texture statistics independent: ! N xi x†i mX b . (6) f (M) = N i=1 x† M b −1 xi i Theorem 1
b F P PROPERTIES 3. M In order to simplify the notations in the proofs, let us define the following quantities: b =M b f p; • M b M−1/2 − Im where Im is the m × m iden• ∆ = M−1/2 M tity matrix;
b f p , which 1. the function f admits a single fixed point, called M verifies b f p) = M b f p. f (M (7)
• X = vec(∆) where X is the vector containing all the elements of ∆ and vec denotes the operator which reshapes the m × n matrix elements into a mn column vector;
2. Let us consider the recurrence relation
• denotes the transpose operator.
b t+1 = M
m b t) “ ” f (M b t) tr f (M
(8)
b f p is unbiased and it is a consistent estimator. Proposition 1 M
b f p. b t −−−→ M Then M t→∞
The point one of the theorem shows that the Maximum Likeb f p exists and is unique. lihood estimate M Relative Error convergence
0
−2
10
Proposition 2 We have the following original results with the above notations: » – √ ` ´ Re(X) law law −−−−−→ N 02m2 , C , where −−→ rep1. N Im(X) N→+∞ resents the convergence in distribution; «2 „ ` ´ m+1 2. N E XX = C1 where m « „ m 1 C1 = ; (9) m + 1 P − m vec(Im )vec(Im )
−4
10
−8
||M − M || / ||M || for the matrix 2−norm t+1 t t ||Mt+1 − Mt|| / ||Mt|| for the matrix 1−norm ||M − M || / ||M || for the Frobenius norm
||M
t+1
t
t
− M || / ||M ||
−6
10
−10
10
t+1
t
Proof 1 This proof will appear in a forthcoming journal paper. It is too long to fit here and we decided to detail the proof of the next proposition. b f p covariance matrix 3.2. M
10
10
3.1. Bias and consistency
t
−12
10
−14
10
−16
10
1
10 Number of Iterations
40
80
Fig. 1. Illustration of the convergence to the fixed point. Figure 1 illustrates the second point of the above theorem. On b t /M b t has been plotb t) − M this figure, the relative error f (M b0 = M b NSCM E (3). ted versus t with initial value M Other simulations have been performed with different initial b 0 (ex : sample Gaussian covariance estimate, matrix with values M uniform PDF, deterministic Toeplitz matrix), each of them conducting to the same value with an extremely fast convergence : b t −−−→ 0 b t) − M f(M t→∞
b f p is only given in (7) as an implicit function of the data: But M b f p . So, it is a significant there is no closed form expression for M b f p by its properties. The statistical propfeature to characterize M erties have never been investigated. The purpose of this paper is
´ ` 3. N E XX† =
„
m+1 m „
«2 C2 where
« m 1 . (10) I m + 1 m2 − m vec(Im )vec(Im ) » – Re(X) Notice that C, which is the covariance matrix of , Im(X) ` †´ ` ´ is fully characterized by the two quantities E XX and E XX . C2 =
b = M + ∆1 . Notice that Proof 2 First, let us define M ∆1 = M1/2 ∆M1/2 .
(11)
b consistency. For large N , ∆1 0m×m because of the M We suppose N to be large enough to ensure the validity of the first order expressions.
Notice that the vec reorganizes the matrix elements ` operator ´ as follows: if H = hij 1≤i,j≤m and vec(H) = (vk )1≤k≤m2 , then hij = vk for k = (j − 1)m + i .
We can write b M
−1
` ` ´´−1 M Im + M−1 ∆1 ` ´ −1 −1 M−1 `Im + M−1 ∆1 ´ −1 − M ∆1 M ´ `Im−1 M − M−1 ∆1 M−1 .
= =
Let us set in equation (18): A = vec
For N large enough, this implies that: N X xi x∗i b m ` −1 ´ . M ∗ N i=1 xi M − M−1 ∆1 M−1 xi
(12)
Hence,
∆1
« N „ xi x∗i m X ` ´ − M. N i=1 x∗i M−1 − M−1 ∆1 M−1 xi
(13)
Let us define yi = M−1/2 xi where yi are gaussian iid complex vectors with identity covariance matrix. Then, N mX yi y∗i ` ´ −Im , N i=1 y∗i Im − M−1/2 ∆M−1/2 yi (14) or equivalently using expression (11),
M−1/2 ∆1 M−1/2
∆
N m X N i=1
y y∗ „ i i∗ « − Im . y ∆y y∗i yi 1 − i ∗ i yi yi
(15)
! .
(19)
Using the Central Limit Theorem (CLT), the right hand side of equation (18) satisfies: » – √ ` ´ Re(A) law N −−−−−→ N 02m2 , G , (20) Im(A) N→+∞ » – Re(A) where G is the covariance matrix of , Im(A) while B at the left hand side, by the Strong Law of Large Numbers (SLLN), has the following property: „ « Di a.s. . (21) B −−−−−→ C2 = Im2 − m E N→+∞ (y∗i yi )2 Thus, from standard probability convergence considerations, we obtain the first point of proposition 2. Moreover, for large N , we have the following equations: ` ´ −1 ` ´ E `XX ´ = C−1 2 E `AA ´ C2 . (22) † E XX† = C−1 C−1 2 E AA 2 Let us now turn to the closed form expression of C2 . « „ D , C2 = Im2 − m E (y∗ y)2
We obtain at the first order, for large N, „ «« N „ m X yi y∗i y∗i ∆yi − Im . ∆ 1+ ∗ N i=1 y∗i yi yi yi
« N „ m X yi y∗i − Im N i=1 y∗i yi
(23)
(16)
where D is the m2 ×m2 matrix defined by D = (dkl )1≤k,l≤m2 with dkl = yp yq y p yq .
To find the explicit expression of ∆ in terms of data, the above expression can be reorganized as:
` ´ Now, as we have y = y1 , . . . , ym ∼ N (0m , Im ) , we can p rewrite yj = rj /2 exp(iθj ) where ∀j ∈ {1 . . . m} , rj and θj are independent variables, with rj ∼ χ2 (2) and θj ∼ U([0, 2π]) .
∆−
« « N „ N „ m X m X yi y∗i y∗i ∆yi yi y∗i − Im . (17) N i=1 y∗i yi y∗i yi N i=1 y∗i yi
To solve this m2 -system, the above equation has to be rewritten as: ! « N „ m X yi y∗i B X vec (18) − Im , N i=1 y∗i yi where • X = vec(∆) , • B = Im2
´ ` We can notice that E exp(i(θp − θq + θq − θp )) = 0 if and only if
N m X Di − , N i=1 (yi y∗i )2 2
2
Then, „ ∀ 1 ≤« p, q, p , q ≤ m , let us set for the elements of D in (23): E=E (y∗ y)2 „ « y p y q y p y q Ekl = E . (y∗ y)2 Thus, with the above notations, we have: ! √ r p r q r p r q ´ ` Ekl = E E exp(i(θp − θq + θq − θp )) `Pm ´2 j=1 rj
(i) (dkl )1≤k,l≤m2
• Di is the m × m matrix defined by Di = with dkl = yp y q y p yq k = p + m(q − 1) with 1 ≤ p, q ≤ m . and l = p + m(q − 1) with 1 ≤ p , q ≤ m
1. p = q = p = q , 2. p = q , p = q et p = p , 3. p = p , q = q et p = q , which is equivalent to: 1. k = l = p + m(p − 1) ,
2. k = p + m(p − 1) , l = p + m(p − 1) and p = p , i.e. k = l 3. k = p + m(q − 1) , l = p + m(q − 1) and p = q .
Curves "PFA−threshold" − CFAR property
N = 20; m = 10;
−1
10
In summary, we obtain the closed form expression of C2 given by (10). In a similar way, we have derived the following results: „ « m 1 where P is the • C1 = m + 1 P − m vec(Im )vec(Im ) m2 × m2 nilpotent matrix defined as follows:
−2
10 PFA
Finally, the non zero elements of the matrix E are: ! rp2 2 1. Ep+m(p−1),p+m(p−1) = E `Pm ´2 = m(m + 1) r j=1 j ! r p r p 1 2. Ep+m(p−1),p +m(p −1) = E `Pm ´2 = m(m + 1) r j j=1 ! rp rq 1 3. Ep+m(q−1),p+m(q−1) = E `Pm ´2 = m(m + 1) j=1 rj
−3
10
−4
Gaussian K−distribution Student−t Cauchy Laplace M estimated with N=20 M Theoretic
10
1
10
2
10
3
10
4
10 threshold λ
5
10
6
10
7
10
8
10
Fig. 2. Relation ”PFA-threshold” for different SIRP noises
6. REFERENCES
Pij = 0 and for 1 ≤ p, p ≤ m : - Pp+m(p−1), p+m(p−1) = 1 , - Pp+m(p −1), p +m(p−1) = 1 , ` ´ ´ ` • N E AA = C1 and N E AA† = C2 .
[1] E. C ONTE , M. L OPS AND G. R ICCI, ”Asymptotically Optimum Radar Detection in Compound-Gaussian Clutter”, IEEE Trans.-AES, 31(2) (April 1995), 617-625.
4. APPLICATIONS
[2] F. G INI, ”Sub-Optimum Coherent Radar Detection in a Mixture of K-Distributed and Gaussian Clutter”, IEE Proc.Radar, Sonar Navig, 144(1) (February 1997), 39-48.
In radar detection, the following test is used to detect, in an mvector observation z, a complex known signal s whose Doppler characteristics are represented here by its steering vector p:
[3] E. JAY, J.P. OVARLEZ , D. D ECLERCQ AND P. D UVAUT, ”BORD : Bayesian Optimum Radar Detector”, Signal Processing, 83(6) (June 2003), 1151-1162
ˆ M) b = Λ(
−1
H1 b z|2 |p† M ≷ λ. −1 −1 † † b b (p M p)(z M z) H0
(24)
In a previous paper [10], we have derived the closed form expression of the relation between the Probability of False Alarm (PFA) and the detection threshold λ, under Gaussian noise assumption b equal to the sample covariance matrix (Wishart disand for M tributed). As illustrated on the Figure 2, results obtained in this paper allowed us to show that the previous relation is surprisingly still b equal to the FPE. Indeed, valid under any SIRP noise and for M the covariance matrix of the Wishart matrix (not done in this paper) is the same near to a multiplicative scalar factor as the FPE covariance matrix established in this paper. Then, the asymptotic behavior of the previous estimates is the same near to a multiplicative scalar factor. 5. CONCLUSIONS AND OUTLOOK We have established in this paper the theoretical statistical properties (unbiasedness, consistency, asymptotic distribution) of the MLE SIRP kernel covariance matrix estimate. Simulation results have confirmed the validity of this work. These properties will be of practical use in radar detection for establishing the statistics of the GLRT-LQ or BORD radar detectors. These adaptive detectors built with the FPE will have the following significant feature: they will be SIRP-CFAR (independent of the texture statistics as well as the structure of the covariance matrix M of the SIRP gaussian kernel).
[4] E. JAY, ”D´etection en Environnement non-Gaussien”, Ph.D. Thesis, University of Cergy-Pontoise / ONERA, France, June 2002. [5] K. YAO ”A Representation Theorem and its Applications to Spherically Invariant Random Processes”, IEEE Trans.-IT, 19(5)(September 1973), 600-608. [6] E. C ONTE , M. L OPS , G. R ICCI, ”Adaptive Radar Detection in Compound-Gaussian Clutter”, Proc. of the European Signal processing Conf., September 1994, Edinburgh, Scotland. [7] F. G INI , M.V. G RECO , AND L. V ERRAZZANI, ”Detection Problem in Mixed Clutter Environment as a Gaussian Problem by Adaptive Pre-Processing, Electronics Letters, 31(14)(July 1995), 1189-1190. [8] E. C ONTE , A. D E M AIO AND G. R ICCI, ”Recursive Estimation of the Covariance Matrix of a Compound-Gaussian Process and Its Application to Adaptive CFAR Detection”, IEEE Trans.-SP, 50(8)(August 2002), 1908-1915. [9] F. G INI , M. G RECO, ”Covariance Matrix Estimation for CFAR Detection in Correlated Heavy Tailed Clutter”, Signal Processing, special section on Signal Processing with Heavy Tailed Distributions, 82(12)(December 2002), 1847-1859. [10] F. PASCAL , J.P. OVARLEZ , P. F ORSTER AND P. L ARZA BAL , ”Constant False Alarm Rate Detection in Spherically Invariant Random Processes”, Proc. of the European Signal processing Conf., September 2004, 2143-2146, Vienna, Austria.