Jean Pierre Delmas

Victor Gabillon

ASYMPTOTIC PERFORMANCE ANALYSIS OF PCA ALGORITHMS BASED ON THE WEIGHTED SUBSPACE CRITERION

Keywords: Principal component analysis, weighted subspace criterion, asymptotic performance analysis

This paper studies the asymptotic distribution of the eigenvectors estimated by some PCA algorithms based on the weighted subspace criterion. This enables us to analyse how the choice of the weighting matrix affects the algorithm's performance, an issue previously overlooked.

INTRODUCTION

Research in subspace and component-based techniques was initiated in the Statistics literature in the middle of the last century through the problem of linear feature extraction solved by the Karhunen-Loève Transform. However the interest of the signal processing community in adaptive subspace and component-based schemes remains strong as it is evident from the numerous articles and reports published in this area each year. In contrast, only a few contributions have been dedicated to the theoretical performance analysis of these adaptive algorithms.

In particular, stochastic approximation gradient-like algorithms devoted to principal and minor subspace analysis (PSA and MSA) or principal and minor component analysis (PCA and MCA) have been mainly studied in stationary environments from two points of view. The first is their convergence in a decreasing step-size situation, using the stability theory of an associated ordinary differential equation (ODE) (see e.g., [START_REF] Oja | Principal components, minor components and linear neural networks[END_REF], [START_REF] Dehaene | Continuous-time matrix algorithms, systolic algorithms and adaptive neural networks[END_REF] and the references therein). In a constant step-size situation, it has been shown [START_REF] Kuan | Convergence of learning algorithms with constant learning rates[END_REF] that the sequence of estimates can be approximated by the associated ODE, in the sense of weak convergence of random processes as the step size tends to zero. The second concerns the asymptotic Gaussian case when both the step size tends to zero and the number of samples tend to ∞ and has been derived for only a few algorithms (see e.g., [START_REF] Yang | Asymptotic distribution of recursive subspace estimators[END_REF], [START_REF] Delmas | Asymptotic performance analysis of subspace adaptive algorithms introduced in the neural network literature[END_REF], [START_REF] Delmas | Performance analysis of an adaptive algorithm for tracking dominant subspace[END_REF]).

The purpose of this paper is to use the approach developed in [START_REF] Delmas | Asymptotic performance analysis of subspace adaptive algorithms introduced in the neural network literature[END_REF], to derive the asymptotic distribution of the estimate issued from a family of adaptive PCA and MCA algorithms derived from the weighted criterion introduced in [START_REF] Oja | Principal component analysis by homogeneous neural networks, Part I: The weighted subspace criterion[END_REF]. One of these gradient-descent like algorithms has just been the subject of a new fast orthogonalization step [START_REF] Bartelmaos | Fast adaptive algorithms for minor component analysis using Householder transformation[END_REF]. From the estimation accuracy point of view, the choice of the weighting matrix crucially affects the algorithm's performance as been shown by simulations in [START_REF] Bartelmaos | Fast adaptive algorithms for minor component analysis using Householder transformation[END_REF]. Hence, it would be of interest to analyse the effect of this weighting matrix on the algorithm's performance and derive, if possible, its optimal value in some sense. This point will be an application of our analysis.

This paper is organized as follows. In Section 2, we give an overview of the family of adaptive PCA and MCA algorithms derived from the weighted criterion introduced in [START_REF] Oja | Principal component analysis by homogeneous neural networks, Part I: The weighted subspace criterion[END_REF]. Then, in Section 3, after presenting a bref review of a general Gaussian approximation result, closed-form expressions of the covariance in the limiting distributions of the eigenvector estimators in a constant step size environment are given by solving Lyapunov equations. Finally, we present briefly in Section 4 some simulation results to illustrate the optimal choice of the weighting matrix.

REVIEW OF THE ALGORITHMS UNDER STUDY

For a given n × n covariance matrix R x = E(xx T ) of a Gaussian distributed, zero mean real random vector x, denote by λ 1 ≥ . . . ≥ λ n the eigenvalues of R x and by v 1 , . . . , v n corresponding normalized eigenvectors. Consider the problem of adaptively estimating r < n normalized eigenvectors corresponding to the r largest [or smallest] distinct eigenvalues (λ 1 , . . . , λ r ) with λ r > λ r+1 [resp. λ n-r+1 , . . . , λ n with

λ n-r > λ n-r+1 ] of R x .
The PCA and MCA stochastic approximation algorithms that are studied in this paper, are derived from a small change in the subspace criterion [START_REF] Yang | Adaptive eigensubspace algorithms for direction or frequency estimation and tracking[END_REF] that reduces its ambiguous solution to a unique set of eigenvectors. If

J(W) def = Tr(ΩW T R x W) = r i=1 ω i w T i R x w i (1)
where

W def = [w 1 , ..., w r ] ∈ R n×r is an orthogonal matrix
and Ω def = Diag(ω 1 , .., ω r ) with ω 1 > ... > ω r > 0, we have the following lemma introduced for PCA in [START_REF] Oja | Principal component analysis by homogeneous neural networks, Part I: The weighted subspace criterion[END_REF], proved in [START_REF] Xu | Least mean square error reconstruction principle for self organizing neural net[END_REF], then recently in [START_REF] Bartelmaos | Fast adaptive algorithms for minor component analysis using Householder transformation[END_REF] with a unifying proof for PCA and MCA.

Lemma 1

The maximization or minimization of J(W) [START_REF] Oja | Principal components, minor components and linear neural networks[END_REF] has the unique solution {±v 1 , ..±v r } or {±v n-r+1 , ..±v n } respectively. All the stationary points of J(W) are saddle points, except the global maximum and global minimum solutions specified above.

The maximization or minimization of the weighted subspace criterion J(W) subject to the constraint W T W = I r can be adaptively solved by constrained gradient-ascent/descent techniques using the instantaneous estimate x k x T k of R x . Since the gradient function of J(.) under orthonormality constraint (see e.g., [START_REF] Dehaene | Continuous-time matrix algorithms, systolic algorithms and adaptive neural networks[END_REF]Ch.7.2]) is equal to R x WΩ -WΩW T R x W, the PCA and MCA stochastic approximation algorithms, that we denote here by weighted Oja algorithm 1 (WOja 1 ), follow clearly

W k+1 = W k ± µ k (x k x T k W k Ω -W k ΩW T k x k x T k W k ).
(2) Another approach consists in using an unconstrained gradient ascent/descent algorithm followed by an orthonormalization step. Using the gradient function R x WΩ of J(.) without constraint, this gives

W k+1 = {W k ± µ k x k x T k W k Ω}G k+1 , (3)
in which G k+1 is a matrix depending on

W k+1 def = W k ± µ k x k x T k W k Ω,
which approximately orthonormalizes the columns of W k+1 . Thus, W k has approximately orthonormal columns for all k. Depending on the form of matrix G k+1 , variants of (3) are obtained. If G k+1 is the symmetric square root inverse of W T k+1 W k+1 , we obtain [10, ch.4] after a power expansion in µ k where all but the zero and first order terms are omitted, the following stochastic approximation algorithm denoted here by weighted Oja algorithm 2 (WOja 2 )

W k+1 = W k ± µ k [x k x T k W k Ω - 1 2 W k ΩW T k x k x T k W k - 1 2 W k W T k x k x T k W k Ω].( 4
)
Instead of deriving a stochastic approximation algorithm from a specific orthonormalization matrix G k+1 , an analogy with Oja's algorithm [START_REF] Oja | Principal components, minor components and linear neural networks[END_REF]

W k+1 = W k + µ k [I n -W k W T k ]x k x T k W k
has been used in [START_REF] Oja | Principal component analysis by homogeneous neural networks, Part I: The weighted subspace criterion[END_REF] to derive the following algorithm denoted here by weighted Oja algorithm 3 (WOja 3 )

W k+1 = W k ±µ k x k x T k W k -W k ΩW T k x k x T k W k Ω -1 . Using y k def = W T k x k , y Ω,k def = Ωy k and z k def = W k y Ω,k
, we note that WOja 1 , WOja 2 and WOja 3 algorithms require only O(nr) operations at each update.

To the best of our knowledge, no complete theoretical performance analysis of these algorithms has been carried out until now. However, the eigenvectors [±v 1 , ..., ±v r ], have been proved to be locally asymptotically stable points for the associated ODE corresponding to the WOja 1 [START_REF] Xu | Least mean square error reconstruction principle for self organizing neural net[END_REF] and WOja 3 [START_REF] Oja | Principal component analysis by homogeneous neural networks, Part II: Analysis and extensions of the learning algorithms[END_REF] PCA algorithms. As for the MCA algorithms, they have been considered as unstable by simulations [START_REF] Bartelmaos | Fast adaptive algorithms for minor component analysis using Householder transformation[END_REF]. Techniques to stabilize these MCA algorithms have been proposed in [START_REF] Bartelmaos | Fast adaptive algorithms for minor component analysis using Householder transformation[END_REF] consisting to an additional step of exact orthonormalization. The denoted MCA-OFRANS algorithms has been proposed where the matrix G k+1 of (3) is the exact symmetric square root inverse of 2), the denoted MCA-OOja algorithm has been derived. However these two MCA algorithms remain sensitive to numerical rounding errors that are eliminated by using an equivalent Householder transform implementation of this last step of exact orthonormalization.

W T k+1 W k+1 . Applied to W k+1 def = x k x T k W k Ω -W k ΩW T k x k x T k W k ) of (

ASYMPTOTICAL PERFORMANCE ANALYSIS

A short review of a general Gaussian approximation result

In this section, we evaluate the asymptotic distributions of eigenvector estimators given by the previous algorithms. For this purpose, we shall use the following result [13, Th.2, p.108]. Consider a constant step size recursive stochastic algorithm

θ k+1 = θ k + µf (θ k , x k ) (5)
with

x k = g(ξ k )
, where ξ k is a Markov chain independent of θ k . Suppose that the parameter vector θ k converges almost surely to the unique asymptotically stable point θ * in the corresponding decreasing step size algorithm. Consider the continuous Lyapunov equation

DC θ + C θ D T + G = O (6)
where D and G are, respectively, the derivative of the mean field and the covariance of the field of algorithm ( 5)

D def = E[ ∂f ∂θ (θ, x k )] θ=θ * , G def = ∞ k=-∞ Cov[f (θ * , x k ), f (θ * , x 0 )].
If all the eigenvalues of the derivative D of the mean field have strictly negative real parts, then, in a stationary situation, when µ → 0 and k → ∞, we have the convergence in distribution

µ -1/2 (θ k -θ * ) L → N (0, C θ ), (7)
where C θ is the unique symmetric solution of the Lyapunov equation (6).

Asymptotic distributions of eigenvector estimators

To characterize the derivative of the mean field and the covariance of the field of the WOja 1 , WOja 2 and WOja 3 algorithms, we use the Vec operator which turns the n × r matrix W into the nr×1 vector parameter vec(W). Thus these three PCA/MCA algorithms with constant step size are written in the form

vec(W k+1 ) = vec(W k ) ± µf (vec(W k , x k x T k )) (8)

Local characterization of the field

Derivative of the field: Using a first order expansion of the field f (., .) with respect to its first argument, we obtain after some algebraic manipulations the following derivatives associated with the three PCA algorithms

D 1 = Ω ⊗ R x -ΛΩ ⊗ I n -I r ⊗ (WΛΩW T ) -[ΛW T ⊗ WΩ]K n,r D 2 = Ω ⊗ R x -ΛΩ ⊗ I n - 1 2 [I r ⊗ (WΛΩW T ) + Ω ⊗ (WΛW T )] - 1 2 [(ΛW T ) ⊗ WΩ) + (ΛΩW T ⊗ W)]K n,r D 3 = I r ⊗ R x -Λ ⊗ I n -Ω -1 ⊗ (WΛΩW T ) -[ΛΩ -1 W T ⊗ WΩ]K n,r ,
where

W def = [±v 1 , .., ±v r ], Λ def = Diag(λ 1 , .
., λ r ) and K l,m denotes the vec-permutation matrix that transforms vec(A T ) to vec(A) for any l × m matrices.

From these expressions, and from those associated with the MCA algorithms, the following theorem is proved in [START_REF] Delmas | Asymptotic performance analysis of PCA algorithms issued from the weighted subspace criterion[END_REF] Theorem 1 The eigenvalues of the derivative D of the mean field of the WOja 1 PCA algorithm are strictly negative real and those of the WOja 2 and WOja 3 PCA algorithms have strictly negative real parts. In contrast, there exist eigenvalues with positive real parts for all the three MCA algorithms.

Covariance of the field: Noting that field (8) is linear in x k x T k and using the covariance of vec(x k x T k ) given by [14, p. 57]

Cov(vec(xx T )) = R x ⊗ R x + (R x ⊗ R x )K n,
n , for independent Gaussian distributed data x k , allow us to derive, after simple algebraic manipulations, the following expressions of the covariance of the field

G 1 = ΛΩ 2 ⊗ R x -2ΛΩ ⊗ (WΛΩW T ) + Λ ⊗ (WΛΩ 2 W T ) + [2ΛΩW T ⊗ WΛΩ -ΛW T ⊗ WΛΩ 2 -ΛΩ 2 W T ⊗ WΛ]K n,r G 2 = ΛΩ 2 ⊗ R x - 1 2 ΛΩ ⊗ (WΛΩW T ) - 3 4 ΛΩ 2 ⊗ (WΛW T ) + 1 4 Λ ⊗ (WΛΩ 2 W T ) + [ 1 2 ΛΩW T ⊗ WΛΩ - 1 4 ΛW T ⊗ WΛΩ 2 - 1 4 ΛΩ 2 W T ⊗ WΛ]K n,r G 3 = Λ ⊗ R x -2ΛΩ -1 ⊗ (WΛΩW T ) + ΛΩ -2 ⊗ (WΛΩ 2 W T ) + [2ΛW T ⊗ WΛ -ΛΩW T ⊗ WΛΩ -1 -ΛΩ -1 W T ⊗ WΛΩ]K n,r .

Solution of the Lyapunov equation

For independent observations x k and for the investigated algorithms, which can be written in a form similar to [START_REF] Delmas | Asymptotic performance analysis of subspace adaptive algorithms introduced in the neural network literature[END_REF] for which the derivative of the mean field have strictly negative real parts (Theorem 1), the hypotheses of the model of Benveniste et al ([13, Th.2, p.108]) are fulfilled. However, the underlying assumption for the results by Benveniste et al is that the solution of the corresponding stochastic approximation type algorithms with decreasing step size, almost surely converges to the unique asymptotically stable solution of the associated ODE. Since the normalized eigenvectors are defined up to a sign, the global attractor W is not unique. Still, the practical use of the Benveniste results in such situation is usually justified by using formally a general approximation result ([13, Th.1, p.107]). Furthermore, the almost sure convergence of the associated decreasing step size algorithms are not strictly fulfilled for these WOja PCA algorithms. This a.s. convergence would need a boundedness condition, whose satisfaction is a challenging problem as discussed in [START_REF] Delmas | Asymptotic performance analysis of subspace adaptive algorithms introduced in the neural network literature[END_REF] If we allow ourselves the Benveniste results in our situation, the Lyapunov continuous equations can be solved exactly. Using the approach of [START_REF] Delmas | Asymptotic performance analysis of subspace adaptive algorithms introduced in the neural network literature[END_REF] that consists in representing the matrices D and G in the orthonormal basis constituted by the nr columns of the r block diagonal matrix Diag(V, .., V)

with V def = [v 1 , .., v n ],
the following closed-form expressions of the covariance matrix C W solution of the Lyapunov equation (6) are derived [START_REF] Delmas | Asymptotic performance analysis of PCA algorithms issued from the weighted subspace criterion[END_REF].

Theorem 2

The covariance matrices C W of the asymptotic distribution that appears in [START_REF] Oja | Principal component analysis by homogeneous neural networks, Part I: The weighted subspace criterion[END_REF] read

C W = 1≤i≤r 1≤k =i≤n b k,i (e r i e r T i ⊗ v k v T k ) + 1≤i =j≤r c i,j (e r i e r T j ⊗ v j v T i ) (9)
with for the WOja 1 , WOja 2 and WOja 3 PCA algorithms respectively

b WOja 1 k,i = λ i λ k (ω i -ω k ) 2(λ i -λ k ) 1 1≤k =i≤r + λ i λ k ω i 2(λ i -λ k ) 1 r<k c WOja 1 i,j = - λ i λ j (ω i -ω j ) 2(λ i -λ k ) 1 1≤i =j≤r b WOja 2 k,i = λ i λ k (ω i -ω k ) 4(λ i -λ k ) 1 1≤k =i≤r + λ i λ k ω i 2(λ i -λ k ) 1 r<k c WOja 2 i,j = - λ i λ j (ω i -ω j ) 4(λ i -λ k ) 1 1≤i =j≤r b WOja 3 k,i = f k,i 1 1≤i =k≤r + λ i λ k 2(λ i -λ k ) 1 r<k c WOja 3 i,j = g i,j 1 1≤i =j≤r
where e r i and 1 A denote the ith unit vector in R r and the indicator function of the condition A and where

f k,i = λiλ k ω k ω i ω k ω i - ω i ω k (λi-λ k ) ω k ω i -λi ω i ω k (1+ ω i ω k ) 2(λ i -λ k )( λ i ω i ω k + λ k ω k ω i ) g i,j = λ i λ j ω i ω j - ω j ω i λ i (2 ω i ω j - ω j ω i )+λ j (2 ω j ω i - ω i ω j ) 2(λ j -λ i )( λ j ω j ω i + λ i ω i ω j )

ILLUSTRATIVE EXAMPLE

Simple measures of performance of the studied PCA algorithms can be derived from Theorem 2 if we suppose that both the first and second moments of the limiting distribution of µ -1/2 (W k -W) are equal to the corresponding asymptotic moments. Among them we have in the steady state,

E(W k ) -W 2 Fro = o(µ), Cov(vec(W k )) = µC W + o(µ) which implies E W k -W 2 Fro = µTr(C W ) + o(µ) ≈ µ 1≤i≤r 1≤k =i≤n b k,i , E w k,i -v i 2 Fro = µTr([C W ] i,i )+o(µ) ≈ µ 1≤k =i≤n b k,i and [15] E W T k W k -I r 2 Fro = o(µ).
We note that from the expressions of b k,i given by Theorem 2, the MSE of each estimated eigenvector w k,1 , . . . , w k,r can be arbitrarily fixed by an appropriate choice of (µ, ω 1 , .., ω r ).

Fig. 1 Learning curves of the MSE E W k -W 2

Fro averaging 1000 independent runs for the 3 PCA WOja algorithms compared to µTr(C W ) (0) for different values 0.9 (1), 0.6 (2), 0.4 (3), 0.1 (4) of ω 1 /ω 2 for n = 3, p = 2 and (λ 1 , λ 2 , λ 3 ) = (1, 0.6, 0.2).

Concerning the minimization of the misadjustment E W k -W 2 Fro , we note that it vanishes for ω 1 = .. = ω r for which the speed of convergence tends to ∞. Consequently a misadjustment/convergence speed tradeoff crops up for the choice of the "optimal" weighting matrix. To solve this dilemma, a solution consists in fixing the theoretical value µTr(C W ) and searching for the value of (ω 1 , . . . , ω r ) that maximizes the speed of convergence. Noting that the set of parameters (µ, ω 1 , . . . , ω r ) = (µ/c, cω 1 , . . . , cω r ) [resp. (µ, ω 1 , . . . , ω r ) = (µ, cω 1 , . . . , cω r )] (for arbitrary value of c > 0) yields the same WOja 1 and WOja 2 algorithms [resp. WOja 3 algorithm], this speed of convergence depends only on the ratios ω 2 /ω 1 , .., ω r /ω r-1 . To illustrate this choice, Fig. 1 displays the similar behavior of the three studied PCA WOja algorithms in a particular scenario. A thorough analysis in relation with (λ 1 , . . . , λ r ) will be presented in [START_REF] Delmas | Asymptotic performance analysis of PCA algorithms issued from the weighted subspace criterion[END_REF].


Jean Pierre Delmas

Soumaya Sallem
email: soumaya.sallem@it-sudparis.eu

Pascal Chevalier
email: pascal.chevalier@fr.thalesgroup.com

SENSITIVIY OF SAIC AND MAIC CONCEPTS TO RESIDUAL FREQUENCY OFFSETS

This paper presents a performance analysis of the single antenna interference cancellation (SAIC) and the multiple antenna interference cancellation (MAIC) techniques in the presence of non null frequency offsets of both the signal of interest (SOI) and interference, that are neither corrected nor compensated, in terms of signal to interference plus noise ratio (SINR) and bit error rate (BER). General theoretical expressions of the SINR and BER are given for BPSK SOI and interference. To obtain engineering insight, the particular case of the SAIC with no return to zero (NRZ) pulse shape filters is considered, where simple expressions are given and analyzed. It is proved in particular that the performance more deteriorates for non null frequency offset of the interference than for non null frequency offset of the SOI. Finally illustrative examples are presented in order to specify the validity domain of our approximations and to quantify the obtained results in the context of the global system for mobile communication (GSM) standard.

INTRODUCTION

For more than a decade, there has been an increasing interest in optimal widely linear (WL) processing [START_REF] Picinbono | Widely linear estimation with complex data[END_REF] in radiocommunication contexts involving rectilinear signals, such as binary phase shift keying (BPSK) signals or quasi rectilinear signals such as continuous phase modulation (CPM) with modulation index 1/2, or such as offset quadrature amplitude modulation (OQAM).

In particular, it has been pointing out in [START_REF] Chevalier | New Insights into optimal widely linear array receivers for the demodulation of BPSK, MSK and GMSK signals corrupted by noncircular interferences -Application to SAIC[END_REF] that SAIC may be performed by WL filters in the context of rectilinear SOI and interference. Then it has been shown in [START_REF] Trigui | Performance bounds for cochannel interference cancellation within the current GSM standard[END_REF], that the Gaussian minimum shift keying (GMSK) modulation whose linearized approximation was introduced in [START_REF] Laurent | Exact and approximate construction of digital phase modulations by superposition of amplitude modulated pulses[END_REF], may be interpreted as a BPSK modulation after a simple algebraic operation of derotation, displaying the great interest of optimal WL filtering for cochannel interference mitigation in the GSM network.

Concerning the performance of these SAIC and MAIC techniques, only a few contributions have appeared in the literature. Among them, [START_REF] Trigui | Performance bounds for cochannel interference cancellation within the current GSM standard[END_REF] have given some bounds on the maximum likelihood sequence estimation (MLSE) for cochannel interference cancellation within the current GSM standard and [START_REF] Chevalier | New Insights into optimal widely linear array receivers for the demodulation of BPSK, MSK and GMSK signals corrupted by noncircular interferences -Application to SAIC[END_REF] have presented some enlightening results about the behavior, properties and performance of the SAIC and MAIC techniques for the reception of a BPSK, MSK or GMSK SOI corrupted by interference of the same kind. In this latter case, no performance analysis has been given concerning the loss in performance in terms of SINR and BER in the presence of residual frequency offsets of the SOI and interference that are practically unavoidable. Therefore, it is of paramount importance to specify how these residual frequency offsets degrade the performance of the SAIC and MAIC receivers.

The purpose of this paper is to quantify this sensitivity. The paper is organized as follows. After the introduction of the observation model and data statistics given in Section 2, Section 3 reviews the optimal BPSK SAIC and MAIC receivers. A performance analysis of these receivers in the presence of residual frequency offsets of the BPSK SOI and interference is presented in Section 4 with a particular attention paid to SAIC with NRZ pulse shape filter. Finally, illustrative examples are given in order to specify the validity domain of our approximations and to quantify the obtained results in the context of the GSM standard in Section 5. Section 6 contains the conclusion.

HYPOTHESES AND DATA STATISTICS

Let us consider an array of N narrow-band sensors. Each sensor is assumed to receive a BPSK SOI corrupted by a noncircular total noise composed of a BPSK interference and background noise. Note that in cellular radiocommunication networks, such an interference may be generated by the network itself (e.g., from signals coming from neighboring cells using the same carrier frequency). Hence, this interference have the same waveform and modulation as the SOI. The complex envelopes of the SOI and the interference are given respectively by

s(t) = ∑ n a n υ(t -nT ), and
j(t) = ∑ n b n υ(t -nT -t j ),
where a n = ±1 and b n = ±1 are independent sequences of independent equiprobable symbols, T is the symbol duration, t j ∈ [0, T ) is the time origin of the interference assuming an optimal sampling time for the SOI at t = kT, k ∈ Z and υ(t) is a square root raised cosine Nyquist filter. The vector x(t) of complex amplitudes of the signals at the output of these sensors is given by

x(t) = µ s s(t)e i2π∆ f s t h s + µ j j(t)e i2π∆ f j t h j + n(t), (1)
where h s and h j are the channel impulse response vectors of the SOI and interference respectively. Note that model (1) assumes propagation channel with no delay spread, which occurs for example, for flat fading channel or ideal propagation. h s = e iφ s s and h j = e iφ j j, where φ s and φ j , s and j correspond to the phases and the steering vectors (such that their first components is one) of the SOI and interference, respectively. ∆ f s and ∆ f j , are the frequency offsets of the SOI and interference respectively. µ s and µ j control the power of the SOI and interference. The background noise n(t) is assumed independent of the SOI and interference, zero-mean Gaussian circular, temporally and spatially white. The sampled observation vector x υ (kT ) def = x(t) ⊗ υ(-t) * / t=kT (where ⊗ is the convolution operation) ob- tained after a matched filtering operation to the pulse shape filter υ(t) and a decimation operation at the symbol rate, is given by

x υ (kT ) = µ s a k e i2π∆ f s kT I 0 h s + µ s ∑ n =k a n e i2π∆ f s kT I k-n h s + µ j ∑ n b n e i2π∆ f j kT J k-n (t j )h j + n υ (kT ), (2)
where

I n def = υ * (-τ)υ(nT -τ)e -i2π∆ f s τ dτ J n (t j ) def = υ * (-τ)υ(nT -τ -t j )e -i2π∆ f j τ dτ. ( 3
)
The extended model is given by

x υ (kT ) = [x υ (kT ) T , x υ (kT ) H ] T with x υ (kT ) = µ s a k h s (k, k) + µ s ∑ n =k a n h s (k, n) + µ j ∑ n b n h j (k, n) + n υ (kT ), (4)
where

h s (k, n) def = [e i2π∆ f s nT I k-n h T s , e -i2π∆ f s nT I * k-n h H s ] T , h j (k, n) def = [e i2π∆ f j (nT +t j ) J k-n (t j )h T j , e -i2π∆ f j (nT +t j ) J * k-n (t j )h H j ] T , n υ (kT ) def = [n υ (kT ) T , n υ (kT ) H ] T .
The second order statistics of the data considered in this paper are defined by R x (k)

def = E[ x υ (kT ) x H υ (kT )
] and we assume that E[ n υ (kT ) n υ (kT ) H ] = η 2 I where η 2 is the mean power of the background noise per sensor and I is the 2N × 2N identity matrix.

SAIC AND MAIC RECEIVERS

Optimal receiver

Under the assumption of equiprobable SOI symbol sequences, without any residual offset and stationary noncircular Gaussian distributed total noise of extended covariance matrix denoted specifically here by R ñ, the MLSE receiver which minimizes the output sequence error rate is given [START_REF] Chevalier | New Insights into optimal widely linear array receivers for the demodulation of BPSK, MSK and GMSK signals corrupted by noncircular interferences -Application to SAIC[END_REF] by a WL filtering of x υ (kT ). This one is the so-called WL spatial matched filter w MLSE defined by

w MLSE def = R -1 n h s , (5)
with

h s def = [h T s , h H s ]
T , whose output y(kT ) is real-valued, followed by a zero threshold detector. This filter is proportional to the filter w MMSE which minimizes the mean square error (MSE) between the output w H x υ (kT ) of the WL filter w = [w T , w H ] T and a k and is given by

w MMSE def = R -1 x r x,a (6)
= [µ s /(1 + µ 2 s h H s R -1 n h s ]R -1 n h s def = β w MLSE , where R x def = E[ x υ (kT ) x H υ (kT )], r x,a def = E[ x υ (kT )a k ].
Note that this WL MMSE gives the same output SINR that the MLSE receiver (5).

Implementation

Without any knowledge about h s , but if a training sequence (a k ) k=1,...,K is available after a synchronisation process, we use the following estimated WL MMSE receiver.

w(K) def = R -1 x (K) r x,a (K), (7)
with

R x (K) def = 1 K ∑ K k=1 x υ (kT ) x H υ (kT ) and r x,a (K) def = 1 K ∑ K k=1 x υ (kT )a k .
Then a sequence (a k ) k=K+1,...,K+L of information symbols is transmitted1 for which the output y(kT ) of the estimated WL MMSE filter is given from ( 4)

for k = K + 1, . . . , K + L by y(kT ) = w(K) H x υ (kT ), k = K + 1, . . . , K + L, (8) = µ s a k w(K) H h s (k, k) + µ s ∑ n =k a n w(K) H h s (k, n) + µ j ∑ n b n w(K) H h j (k, n) + w(K) H n υ (kT ). (9)

PERFORMANCE ANALYSIS

Assumptions

The performance analysis of such a scheme is challenging because y(kT ) given by ( 9) is a random variable (RV) depending on the RVs (a k , b k , n υ (kT )) k=1,...,K of the training period through the estimate w(K) and on the RVs (a k , b k , n υ (kT )) k=K+1,...,K+L of the information period. To simplify the performance analysis, we assume that for K "sufficiently large", the loss of SINR with respect to those obtained with

w(K) def = R -1 x (K)r x,a (K) (10)
is "very weak", where

R x (K) def = 1 K K ∑ k=1 R x (k) and r x,a (K) def = 1 K K ∑ k=1 r x,a (k) with r x,a (k) def = E[ x υ (kT )a k ].
Note that in the absence of frequency offset for which x υ (kT ) is stationary, this assumption is valid for relatively small values of K. More precisely, it is proved in [5, Chap. 6.1.2] that to obtain a loss of SINR lower than 3dB, K = 2N samples are needed for not too small output SINR. We will consider throughout our performance analysis that w(K) can be replaced by w(K) given by (10), without affecting the SINR at the output y(kT ).

Theoretical SINR, SNR and INR

With w(K) given by (10) in (8) and with the independence of a k , b k and n υ (kT ), for k = K + 1, . . . , K + L, the SINR, the interference to noise ratio (INR) and the signal to noise ratio (SNR) at the output of the estimated WL MMSE filter are given by

SINR(kT ) = µ 2 s [ w H (K) h s (k, k)] 2 /µ 2 s ∑ n =k [ w H (K) h s (k, n)] 2 + µ 2 j ∑ n [ w H (K) h j (k, n)] 2 + η 2 || w(K)|| 2 , (11) INR(kT ) = µ 2 s ∑ n =k [ w H (K) h s (k, n)] 2 + µ 2 j ∑ n [ w H (K) h j (k, n)] 2 /η 2 || w(K)|| 2 , (12) SNR(kT ) = µ 2 s [ w H (K) h s (k, k)] 2 /η 2 || w(K)|| 2 . ( 13
)

Theoretical BER

But the relevant criterion to evaluate the loss in performance is the output BER. Using ( 9)

y(kT ) = a k α k + ∑ n =k a n α n,k + ∑ n b n β n,k + n k , ( 14
)
with α k def = µ s w H (K) h s (k, k), α n,k def = µ s w H (K) h s (k, n), β n,k def = µ j w H (K) h j (k, n) and n k def = w H (K) n υ (kT )
which is a real valued zero-mean Gaussian RV of variance η 2 || w(K)||2 and conditioning on specific values of the sequence (a n ) n =k , (b n ), the BER at time kT is clearly given by BER(kT

) = Q α k -∑ n =k a n α n,k -∑ n b n β n,k η 1/2 2 || w(K)|| , where Q(v) def = +∞ v 1 √
2π e -u 2 /2 du. This gives from the total probability formula and from the assumption of equiprobable sequences

(a n ) n =k , (b n ) contained in (14), the following BER BER(kT ) = 1 2 I+J 2 I+J ∑ s=1 Q α k -∑ i =k a s i α i,k -∑ j b s j β j,k η 1/2 2 || w(K)|| (15) where (. . . , a s k-1 , a s k+1 , . . . , b s k-1 , b s k , b s k+1 , . . .)
denotes the 2 I+J different (I + J)-uplets 2 of binary symbols (a i ) i∈I , (b i ) j∈J with (a i ) i∈I are the SOI inter symbol interference and (b i ) j∈J are the interference symbols associated with the SOI symbol a k .

Particular case of SAIC with NRZ pulse shape filter

To give an interpretation of the loss in performance from the different theoretical expressions (11), ( 12), ( 13) and (15) that are lacking of engineering insights, we consider the particular case of SAIC (N = 1) with t j = 0 and with the following NRZ filter

υ(t) = 1/ √ T if -T 2 ≤ t ≤ T 2 0 elsewhere
for which the values of I n and J(t j ) defined in (3), are given by

I n = sin(π∆ f s T ) π∆ f s T δ (n) and J n def = J n (0) = sin(π∆ f j T ) π∆ f j T δ (n).
Consequently x υ (kT ) given by (4) reduces to

x υ (kT ) = µ s a k h s (k, k) + µ j b k h j (k, k) + n υ (kT ), (16)
with no interference inter symbol coming from the SOI and interference, where h s (k, k) = I 0 (e i2π∆ f s kT e iφ s , e -i2π∆ f s kT e -iφ s ) T and h j (k, k) = J 0 (e i2π∆ f j kT e iφ j , e -i2π∆ f j kT e -iφ j ) T . This gives

R x (K) = µ 2 s I 2 0 1 α s (K)e 2iφ s α * s (K)e -2iφ s 1 + µ 2 j J 2 0 1 α j (K)e 2iφ j α * j (K)e -2iφ j 1 + η 2 I r x,a (K) = µ s I 0 α s (K)e iφ s α * s (K)e -iφ s ,
where α s (K) kT . Using these expressions of R x (K) and r x,a (K) in (10), the inversion matrix lemma gives after straightforward but cumbersome algebraic derivations w(K) = (w(K), w(K) * ) T with

def = 1 K ∑ K k=1 e i4π∆ f s kT , α j (K) def = 1 K ∑ K k=1 e i4π∆ f j kT and α s (K) def = 1 K ∑ K k=1 e i2π∆ f s
w(K) ∝ (1 + ε s + ε j )α s (K)e iφ s -[ε s α s (K)e 2iφ s + ε j α j (K)e 2iφ j ]α * s (K)e -iφ s , ( 17
)
where ∝ denotes proportional up to a real-valued constant and ε

s def = µ 2 s I 2 0 η 2 and ε j def = µ 2 j J 2 0 η 2 .
Consequently the general expressions of (11), ( 12), ( 13) and (15) reduce to

SINR(kT ) = 2ε s [ℜ(w * e iψ s,k e iφ s )] 2 2ε j [ℜ(w * e iψ j,k e iφ j )] 2 + |w(K)| 2 , (18) INR(kT ) = 2ε j [ℜ(w * e iψ j,k e iφ j )] 2 |w(K)| 2 , ( 19
) SNR(kT ) = 2ε s [ℜ(w * e iψ s,k e iφ s )] 2 |w(K)| 2 , ( 20
BER(kT ) = 1 2 Q( SNR(kT ) + INR(kT )) + Q( SNR(kT ) -INR(kT )) . ( 21
)
We assume to get engineering insight that the residual frequency offsets are sufficiently weak3 such that

(K + L)|∆ f s |T 1 and (K + L)|∆ f j |T 1 which implies that ε s ≈ µ 2 s η 2 and ε j ≈ µ 2 j η 2 .
Furthermore, we assume that the power of the interference is large with respect to the power of the background noise, i.e., ε j 1. We note [START_REF] Chevalier | New Insights into optimal widely linear array receivers for the demodulation of BPSK, MSK and GMSK signals corrupted by noncircular interferences -Application to SAIC[END_REF] that in the absence of residual offset and for ψ def = φ j -φ s = 0, the SAIC receiver has no rejection capability, i.e., the INR at the output of the optimal SAIC receiver is INR= 2ε j = 2µ 2 j η 2 . So we assume that ψ is such that |ψ| (K + L)|∆ f j |T and |ψ| (K + L)|∆ f s |T . To proceed on, we must specialize the case where there is residual offset of interference or SOI only.

Frequency offset of interference only

In this case ∆ f s = 0, and thus α s (K) = α s (K) = 1 and ψ s,k = 0. With respect to the SINR, INR and SNR given by the optimal SAIC without residual offset, we obtain from ( 18), ( 19) and ( 20)

SINR(kT ) ≈ SINR 1 + 2ε j [π{2k -K -1}∆ f j T ] 2 , ( 22
) whereas INR(kT ) ≈ 2ε j [π{2k -K -1}∆ f j T ] 2 with INR ≈ 1 2ε j tan 2 (ψ)
1 and SNR(kT ) ≈ SNR ≈ 2ε s sin 2 (ψ) for k = K +1, . . . , K +L. We see that the SNR is preserved in contrast to the INR and SINR which strongly degrade in the presence of residual offset of interference for ε j 1. This loss in performance naturally increases with |∆ f j |T and (using 2k -K = 2(k -K) + K) with the position k -K of the information symbol, but also strongly with the input INR ε j and the size K of the training sequence.

Frequency offset of SOI only

In this case ∆ f j = 0, and thus α j (K) = 1 and ψ j,k = 0. The simplification of the expressions of SINR(kT ), INR(kT ) and SNR(kT ) is more complex because we must distinguish ψ = ±π/2 from the case where ψ is not in the neighborhood of ±π/2. For φ = ±π/2 we have

SINR(kT ) ≈ SINR 1 - (2π∆ f s T ) 2 ((K + 1) 2 + 4k 2 ) 8 , (23) whereas INR(kT ) ≈ 1 2ε j (π(K + 1)∆ f s T ) 2 with INR = 0 and SNR(kT ) ≈ SNR 1 -4(πk∆ f s T ) 2 with SNR= 2ε s for k = K + 1, . . . , K + L.
For ψ not in the neighborhood of ±π/2, we have

SINR(kT ) ≈ SINR 1 - 4πk∆ f s T tan(ψ) , ( 24
)
whereas INR(kT

) ≈ INR 1 + 4π(K+1)∆ f s T sin(2ψ) with INR ≈ 1 2ε j tan 2 (ψ) 1 and SNR(kT ) ≈ SNR 1 -4πk∆ f s T tan(ψ)
with SNR= 2ε s sin 2 (ψ) for k = K + 1, . . . , K + L. In a similar way for ψ = ±π/2, the loss in performance increases with |∆ f s |T , with the position k -K of the information symbol and the size K of the training sequence, but no longer with ε j . But for ψ not in the neighborhood of ±π/2, depending of the sign of ∆ f s T / tan(ψ), the SINR can locally decrease or increase for (K + L)|∆ f s |T 1, independently of ε j .

Comparisons

Comparing the loss in performance due to residual offsets of the interference or the SOI from the aforementioned expressions is not easy. But comparing ( 22), ( 23) and (24) using the assumption that ε j 1, we see that a loss in SINR greater than e.g. 3dB is obtained for weaker value of |∆ f j |T than of |∆ f s |T , and larger is ε j , weaker is this ratio |∆ f j |/|∆ f s |. Consequently, under the assumption that ε j 1, the SAIC is less sensitive to residual offsets of the SOI than of the interference. Note that this property is similar to the sensitivity to steering errors in spatial beamforming for which the loss in SINR is much more sensitive to interference steering error than to SOI steering error. This property will be specified and confirmed in the next section in other scenarios.

ILLUSTRATIONS

Throughout this section we use ε s = 10dB and ε j = 20dB with K = 26 and L = 58 (GSM standard).

Validation of the assumptions

Here, ψ def = φ j -φ s is fixed to π/3. Figs. 1 and2 show the theoretical (exact (21) and approximate ( 22) and ( 24)) and empirical (Monte Carlo with 10000 runs) SINR as a function of the position k -K of the information symbol for ∆ f s = 0, ∆ f j T = 5 10 -4 and ∆ f j = 0, ∆ f s T = 5 10 -4 respectively, for the SAIC receiver with NRZ pulse shape filter. The SINRs at the middle of the burst (30th symbol) are plotted in Figs. 3 and4 as a function of ∆ f j T and ∆ f s T for ∆ f s T = 0 and ∆ f j T = 0, respectively. We see from these figures that all the theoretical SINR fit the empirical ones. However, we note in Fig. 1 that the empirical SINR perfectly fits the theoretical SINR in contrast to Fig. 2. This shows that the assumptions of Subsection 4.1 is better justified in case of frequency offset of interference than of SOI. This is likely due to the estimated filter w(K) in ( 17) that is more perturbed to frequency offset of SOI than of interference. Furthermore our simplified expressions ( 22) and ( 24) are roughly valid up to ∆ f j T = 10 -3 and ∆ f s T = 5 10 -3 respectively.

Practical applications

We consider now the practical case in which the phases of the SOI and interference for the SAIC receiver, and the phases and the directions of arrival for the MAIC receiver are totally unknown. Figs. [START_REF] Monzingo | Introduction to adaptive arrays[END_REF] shows the theoretical BER given by (21) and the empirical BER (Monte Carlo with 20000 runs) averaged over the phases φ s and φ j , and the L information symbols for the SAIC receiver with an NRZ pulse shape filter, as a function of ∆ f j T for ∆ f s T = 0 and ∆ f s T for ∆ f j T = 0. This BER is compared to those obtained with a raise cosine pulse shape filter with a roll-off of 0.22 in Fig. 6. Note that it is derived by 20000 Monte Carlo runs because the derivation from (21) is too computationally demanding due to the truncation of υ(t) to 15 samples. We see from Fig. 5 that the averaged theoretical BER fits the empirical ones and that the raise cosine pulse shape filter roughly does not modify the loss in performance. Note that the SAIC is less sensitive to residual offsets of the SOI than of the interference. More precisely, we note an increase of the BER of ten per cent for ∆ f j T = 1.5 10 -4 or for ∆ f s T = 5 10 -4 , i.e., 40Hz or 135Hz for the GSM standard (T = 1/270) ms, respectively. We note that for ε j = 10dB and ε j = 30dB these values become ∆ f j T = 4 10 -4 and ∆ f j T = 1.2 10 -4 , whereas the value of ∆ f s T keeps the same value as predicted by our approximations ( 22) and (23).

Finally Figs.7 shows the averaged theoretical BER given by (21) and the empirical BER (Monte Carlo with 10000 runs) with respect to the phases φ s and φ j , the directions of arrival θ s and θ j and the L information symbols for the MAIC receiver (with N = 2 omnidirectional sensors equispaced half a wavelength apart) with an NRZ pulse shape filter as a function of ∆ f j T for ∆ f s T = 0 and ∆ f s T for ∆ f j T = 0. We see that MAIC receiver roughly presents the same sensitivity than the SAIC receiver.

CONCLUSION

In this paper, we have presented a theoretical performance analysis of the loss in performance of the SAIC and MAIC receiver in the presence of residual frequency offsets of BPSK SOI and interference. We have proved for NRZ pulse shape filters than the SAIC receiver is less sensitive to residual offset of the SOI than of the interference for strong interference w.r.t. background noise and some rules of thumb about the normalized tolerable residual offsets have been given. In particular for a burst structure similar to those of the GSM, an increasing of BER of ten per cent is obtained for ∆ f j = 40Hz or ∆ f s = 135Hz. These properties of sensitivity have been extended by Monte Carlo experiments to raise cosine pulse shape filters, MAIC receivers and extensive scenarios of phases and directions of arrival.

Extension of this work to MSK and GMSK signals will be considered in a future contribution to analyze the behavior of the SAIC and MAIC receivers for cochannel mitigation with nonnull residual offsets in the GSM network.

) with ψ s,k def =

 def 2π∆ f s kT and ψ j,k def = 2π∆ f j kT and[START_REF] Chevalier | New Insights into optimal widely linear array receivers for the demodulation of BPSK, MSK and GMSK signals corrupted by noncircular interferences -Application to SAIC[END_REF]

Fig. 1 Fig. 2

 12 Fig.1SINR for ∆ f j T = 5.10 -4 and ∆ f s T = 0.

Fig. 3

 3 Fig.3 SINR at 30-th symbol for ∆ f s T = 0.

Fig. 4

 4 Fig.4SINR at 30-th symbol for ∆ f j T = 0.

Fig. 6

 6 Fig.5 Averaged BER, NRZ pulse

Fig. 7

 7 Fig.7 Averaged BER, NRZ pulse for MAIC (N = 2)

In practice in each burst, the training sequence is preceded and followed by this information sequence with e.g., K =

and L = 58 in the GSM standard.

where I = J -1 with J is such that JT represents the length of the "significant" part of υ(τ)υ * (τ -t)dτ.

Note that when (K + L)|∆ f s |T ∼ 1 or (K + L)|∆ f j |T ∼ 1, the SOI or the interference are seen by the SAIC receiver as second order circular and no rejection capability is possible.


Jean-Pierre Delmas
email: jean-pierre.delmas@it-sudparis.eu

Houcem Gazzah
email: hgazzah@sharjah.ac.ae

CRB ANALYSIS OF NEAR-FIELD SOURCE LOCALIZATION USING UNIFORM CIRCULAR ARRAYS

Keywords: Cramer Rao bound, near field source localization, uniform circular array, isotropy

This paper is devoted to the Cramer Rao bound (CRB) on the azimuth, elevation and range of a narrow-band near-field source localized by means of a uniform circular array (UCA), using the exact expression of the time delay parameter. After proving that the conditional and unconditional CRB are generally proportional for constant modulus steering vectors, we specify conditions of isotropy w.r.t. the distance and the number of sensors. Then we derive very simple, yet very accurate non-matrix closed-form expressions of different approximations of the CRBs.

INTRODUCTION

A considerable literature has been dedicated to the CRB on the direction of arrival (DOA) of narrow-band sources. However, most of it is limited the far-field case, meaning that a planar wavefront is impinging on the sensors array. The nearfield assumption implies that the curvature of the waves has to be taken into account, resulting into more complicated models parameterized by the source DOAs and range. Hence, we can identify only very few studies of this kind. Inspired by the subspace-based DOA algorithms, early ones were based on an approximate propagation model based on second-order Taylor expansion of the time delay parameter [START_REF] Grosicki | A weighted linear prediction method for near field source localization[END_REF][START_REF] El Korso | Conditional and unconditional Cramer Rao bounds for near-field Source localization[END_REF]. Only lately has the exact time delay formula been used [START_REF] Begriche | Exact Cramer Rao bound for near field source localization[END_REF] for the only uniform linear arrays (ULA).

The UCA has been popular with many applications [START_REF] Wang | Signal processing algorithms for fire localization using temperature sensor arrays[END_REF][START_REF] Giridhara | Rapid Localization of damage using a circular sensor array and Lamb wave based triangulation[END_REF][START_REF] He | Near-field beamforming analysis for acoustic emission source localization[END_REF] and dedicated estimation algorithms have been proposed for far-field source localization [START_REF] Mathews | Eigenstructure techniques for 2-D angle estimation with uniform circular array[END_REF][START_REF] Chan | MUSIC and maximum likelihood techniques on two-dimensional DOA estimation with uniform circular array[END_REF]. Isotropy is a compelling feature of UCA. Not only does it mean that the CRB on the azimuth and elevation of a single source is uniform for all azimuths [START_REF] Gazzah | Optimum ambiguity free directional and omnidirectional planar antenna arrays for DOA estimation[END_REF], but also that the two estimates are uncorrelated and that the CRB is not affected by the indeterminacy about the wave speed of propagation [START_REF] Oktel | Source localization with isotropic arrays[END_REF]. We are interested in analyzing the behavior of UCA for near-field sources, and in particular to what extend is isotropy maintained. Furthermore, very simple non-matrix closed-form expressions of ap-proximations of the CRB on the azimuth, elevation and range are given, where the sixth-order approximation of the CRB on the azimuth turns out to be very accurate.

The paper is organized as follows. Section 2 formulates the problem and specifies the data model. Section 3 is devoted to the CRB derivations. After proving that the so-called conditional and unconditional CRB are generally proportional for constant modulus steering vectors, we specify conditions of isotropy w.r.t. the range and the number of sensors. We derive very simple non-matrix closed-form expressions of different approximations of the CRB. Finally some numerical illustrations are given in Section 4. Corrupted by an additive noise with complex envelope n(k), the complex envelope of the signals collected by the array of sensors is modelled as

DATA MODEL

x k = s k a(α) + n k , k = 1, ..., K,
where a(α) is the so-called steering vector parameterized by α = [θ, φ, r] T . Referenced to the origin, its is given by a(α) = [e iτ1(α) , ..., e iτp(α) , ..., e iτP (α) ] T where τ p (α) = 2π dp λ0 with d p = SO -SC p and λ 0 is the propagation wavelength. In fact, τ p is given by

τ p (α) = 2π r λ 0 1 -1 -2 r 0 r cos θ p sin φ + r 2 0 r 2 .
The usual statistical properties about n k and s k are the following: (i) s k and n k are independent, (ii) (n k ) k=1,...,K are independent, zero-mean circular Gaussian distributed with covariance σ 2 n I P , (iii) (s k ) k=1,...,K are assumed to be either deterministic unknown parameters (the so-called conditional or deterministic model), or independent zero-mean circular Gaussian distributed with variance σ 2 s (the so-called unconditional or stochastic model).

CRB DERIVATION

General expression of the CRB

General compact expressions of the CRB, concentrated on the DOA parameter alone have been derived for these two models of sources (see e.g. [START_REF] Stoica | Performances study of conditional and unconditional direction of arrival estimation[END_REF] for one parameter per source and [10, Appendix D] for several parameters per source). Specialized to a single source with one or several parameters per source, these expressions CRB det (α) and CRB sto (α) have all been studied independently in the literature (see e.g. the recent papers [START_REF] Vu | A Cramer Rao bounds based analysis of 3D antenna array geometries made from ULA branches[END_REF] and [START_REF] Begriche | Exact Cramer Rao bound for near field source localization[END_REF] which even concludes by "extension of this work for stochastic sources is under consideration").

In fact, it can be proved (proof is not given here by lack of space) that the expressions CRB det (α) and CRB sto (α) for arbitrary parametrization α = [α 1 , ...., α L ] T of the steering vector a(α) related to array geometry or polarization such that a(α) 2 = P , are proportional1 , an issue previously overlooked. More precisely:

CRB sto (α) = 1 + σ 2 n P σ 2 s CRB det (α) = [FIM(α)] -1 (1)
where

σ 2 s def = 1 K K k=1 |s 2 k | for the deterministic model of sources, with FIM(α) = c σ P D H (α)D(α) -D H (α)a(α)a H (α)D(α) (2) where c σ def = 2Kσ 4 s σ 2 n (σ 2 n +P σ 2 s ) and D(α) def = ∂a(α) ∂α1 , ..., ∂a(α) ∂αL . In the addressed problem, D(α) = ∂a(α) ∂θ , ∂a(α) ∂φ , ∂a(α) ∂r = [a ′ θ , a ′ φ , a ′ r ].
Thanks to (1), we only consider the stochastic source model.

Exact Fisher information matrix

The elements of the FIM (2) are, for i, j = 1, 2, 3, given by

FIM i,j (α) = c σ [P (a ′ H αi a ′ αj ) -(a ′ H αi a)(a ′ H αj a) * ],
where a ∂αi . Consequently, we get Taylor series expansions of 1/β p and 1/ β p w.r.t. ǫ = r 0 /r, followed by elementary trigonometric relations, show that all elements of the FIM depend on the azimuth θ only through the sums P p=1 cos kθ p and P p=1 sin kθ p for k = 1, 2, ..., which can be easily simplified thanks to

FIM 1,1 (α) = c ′ σ sin 2 φ   P P p=1 u 2 p - P p=1 u p 2   FIM 2,2 (α) = c ′ σ cos 2 φ   P P p=1 v 2 p - P p=1 v p 2   FIM 1,2 (α) = c ′ σ sin(2φ) 2 P p=1 u p P p=1 v p -P P p=1 v p u p FIM 3,3 (α) = c ′′ σ   P P p=1 (1 + w p ) 2 - P p=1 [1 + w p ] 2   FIM 1,3 (α) = c
P p=1 e ikθp = P e ikθ if k/P ∈ N 0 otherwise (3)
It turns to be that while the FIM is periodic in θ of period 2π/P , as one may expect, it is not constant anymore, i.e. the UCA is no longer strictly isotropic in the near-field.

Isotropy of the UCA

We analyze the way isotropy, guaranteed at the antenna farfield [START_REF] Gazzah | Optimum ambiguity free directional and omnidirectional planar antenna arrays for DOA estimation[END_REF], is deteriorated as a result of a decreasing source range r or a decreasing number of sensors P . Careful examination of the Taylor expansion of the elements of the FIM w.r.t. ǫ, where only the θ dependence is retained, yields to

FIM i,j (α) = ∞ k=0 P p=1 g (i,j) k (cos θ p , sin θ p ) ǫ k , (4)
where g (i,j) k is a polynomial expression of cos θ p and sin θ p of degree k + 2, k + 1 or k for (i, j = 1, 2), (i = 1, 2, j = 3) or (i = j = 3) respectively. By linearizing this polynomial, we have for example for i, j = 1, 2:

g (i,j) k (cos θ p , sin θ p ) = k+2 ℓ=0 c (i,j) ℓ,k cos(ℓθ p )+ k+2 ℓ=1 s (i,j) ℓ,k sin(ℓθ p )
where c (i,j) 0,k = 0 for odd degrees of g (i,j) k . Then, using (3), focusing on θ and carefully studying the first terms of the Taylor expansion in ǫ (4), we obtain

FIM i,j (α) = ⌊(P -3)/2⌋ k=0 b i,j 2k ǫ 2k + ∞ k=P -2 b i,j k (θ)ǫ k (5) FIM i,j (α) = ∞ k=P -2 b i,j k (θ)ǫ k (6) FIM 2,3 (α) = P -1 k=3 b i,j k ǫ k + ∞ k=P b i,j k (θ)ǫ k (7) FIM 3,3 (α) = ⌊(P -1)/2⌋ k=2 b i,j 2k ǫ 2k + ∞ k=P b i,j k (θ)ǫ k , (8)
for P ≥ 4 and i = j = 1 or i = j = 2 (5), i = 1, j = 2, 3

, where the terms b i,j k and b i,j 2k do not depend on θ. Note that zero-order terms of the block [FIM(α)] (1;2,1:2) derived from ( 5) and [START_REF] He | Near-field beamforming analysis for acoustic emission source localization[END_REF] give the FIM(θ, φ) of the far-field case. Using the decoupling of θ and φ of the far-field case that we uncover by ( 6), we have CRB FF (θ) = 1/b 1,1 0 and CRB FF (φ) = 1/b2,2 0 . To illustrate the behavior of the FIM w.r.t. P , FIM1,1(α) cσcr 0 sin 2 φ is given by the following expressions for P = 3, 4, 5, 6:

P 3 1 -ǫ 2 sin φ cos 3θ + o(ǫ 2 ) 4 1 -ǫ 2 (cos 2 φ -sin 2 φ sin 4θ) + o(ǫ 2 ) 5 1 -ǫ 2 cos 2 φ -ǫ 3 sin 3 φ cos 5θ + o(ǫ 3 ) 6 1 -ǫ 2 cos 2 φ -ǫ 4 (1 -3 sin 2 φ + 2 sin 4 φ -sin 4 φ cos 6θ) + o(ǫ 4 ) with c r0 def = 2π r0 λ0 2 P 2
2 and where lim ǫ→0 o(ǫ)/ǫ = 0. Consequently, from ( 5)-( 8), the following can be concluded: For a number P of sensors, the FIM of the UCA does not depend on the azimuth up to the order P -3 in r 0 /r. More precisely, we see from [START_REF] He | Near-field beamforming analysis for acoustic emission source localization[END_REF], that the parameters θ and (φ, r) are decoupled up to the order P -1 in r 0 /r. Consequently, the azimuth CRB does not depend on the azimuth up to the order P -1 in r 0 /r, in contrast to the elevation and range CRBs that require an order smaller or equal to P -3. So, for r or P fixed, isotropy increases when P or r increases, respectively, and the azimuth CRB is much less sensitive to the azimuth than those of the elevation and range.

Closed-form expression of the CRB

To further improve our understanding of the near-field CRBs, we investigate different closed-form expressions of the approximations of the CRB on α. Note that to obtain these approximations, we have to elaborate a little bit the necessary order of the Taylor expansions of each term of FIM(α) because the order in ǫ of these terms may be different.

For P > 8, we have proved that 2 This expression of FIM(α) allows us to derive after some manipulations of FIM -1 (α) the following closed-form expressions of the CRBs:

FIM 1,1 (α) c σ c r0 = sin 2 φ 1 -ǫ 2 cos 2 φ + ǫ 4 g 1 (sin 2 φ) +ǫ 6 g 2 (sin 2 φ) + o(ǫ 7 ) FIM 1,2 (α) = o(ǫ 7 ) FIM 1,3 (α) = o(ǫ 7 ) FIM 2,2 (α) c σ c r0 = cos 2 φ 1 -ǫ 2 1 - 5 2 sin 2 φ + o(ǫ 3 ) FIM 3,3 (α) c σ c r0 = sin 2 φ 16r 2 0 sin 2 φ + ǫ 2 g 3 (sin 2 φ) ǫ 4 + o(ǫ 7 ) FIM 2,3 (α) c σ c r0 = - sin φ cos φ r 0 1- 7 8 sin 2 φ ǫ 3 + o(ǫ 3 ),
CRB(θ) = CRB FF (θ) 1 + cos 2 φ r 2 0 r 2 + sin 2 φ cos 2 φ r 4 0 r 4 + h 1 (sin 2 φ) r 6 0 r 6 + o r 7 0 r 7 (9)
CRB(φ) = CRB FF (φ) 1+h 2 (sin 2 φ) r 2 0 r 2 + o r 2 0 r 2 (10)
CRB(r)= 16 c σ c r0 sin 4 φ r 4 r 2 0 1+h 3 (sin 2 φ) r 2 0 r 2 +o r 2 0 r 2 , (11)
where

CRB FF (θ) = 1/(c σ c r0 sin 2 φ) CRB FF (φ) = 1/(c σ c r0 cos 2 φ)
denote the CRB on the azimuth and elevation in the far field case [START_REF] Gazzah | Optimum ambiguity free directional and omnidirectional planar antenna arrays for DOA estimation[END_REF], and where we have introduced Using the decoupling between θ and (φ, r), note that the second-order (resp. fourth-order) expansion in r 0 /r of CRB(θ) in ( 9) is still valid for only P > 4 (resp. P > 6).

Interestingly, when the source is known to be in the (x, y) plane, we deduce from the FIM(α) for φ = π/2 that

CRB(θ) = CRB FF (θ) 1 + o r 7 0 r 7 CRB(r) = 16 c σ c r0 r 4 r 2 0 1 - r 2 0 2r 2 + o r 2 0 r 2 ,
where here CRB FF (θ) = 1/c σ c r0 . Note that the first order term r 0 /r vanishes in all the expressions of the CRB for P > 6. Furthermore, for a source located in the plane (x, y), the CRB on the azimuth is very insensitive to the range. This contrasts with the near field expression of the CRB for the ULA [START_REF] Begriche | Exact Cramer Rao bound for near field source localization[END_REF] for which, CRB(θ) includes a first order term in r 0 /r and thus is much more sensitive to the range than for the UCA. Finally, note the simplicity of our closed-form expressions w.r.t. the complicated non interpretable closed-form expressions obtained for the ULA [START_REF] El Korso | Conditional and unconditional Cramer Rao bounds for near-field Source localization[END_REF][START_REF] Begriche | Exact Cramer Rao bound for near field source localization[END_REF].

ILLUSTRATIONS

In order to characterize the level of isotropy, we introduce a non-isotropy criterion

ρ = sup θ |CRB(θ) -CRB(θ)| CRB(θ)
(where CRB(θ) denotes the mean of CRB(θ) w.r.t. θ) that we illustrate in Fig. 2 for a UCA with half-wavelength intersensors spacing, and a source emitting at 1 MHz 3 . This figure shows that the isotropy is much more sensitive to P than to r. The validity of some approximate closed-form expressions of the CRB is illustrated for a source located with an azimuth θ = 70 • and elevation φ = 70 • . Fig. 3 and4 compare the approximate ratios CRB(θ)/CRB FF (θ) and CRB(φ)/CRB FF (φ) given by ( 9) and [START_REF] Abeida | Efficiency of subspace-based DOA estimators[END_REF] to the exact ones and Fig. 5 compares the approximate CRB(r) [START_REF] Vu | A Cramer Rao bounds based analysis of 3D antenna array geometries made from ULA branches[END_REF] to the exact one. We see from these three figures that the proposed approximations are very accurate from the ratio r/r 0 = 2 for P = 10 and from r/r 0 = 4 for P = 7.

3 Note that these characteristics impact only r 0 , but not the different ratios of CRB.

AFig. 1

 1 Fig.1 Uniform circular array and source DOAs.

v p w p where u p def =- 1 /2 p , β p def = 1 -

 def11 sin θ p / β p , v p def = cos θ p / β p , w p def = ǫv p sin φ -β 2ǫ cos θ p sin φ + ǫ 2 and ǫ def = r0r . There, we also define constants c

with g 1 (sin 2 φ) def = 1 -3 sin 2 φ + 2 sin 4 φ g 2 (

 112 sin 2 φ) def = -1 + 6 sin 2 φ -10 sin 4 φ + 5 sin 6 φ g 3 (sin 2 φ) def = 16 -33 sin 2 φ + 35 2 sin 4 φ.

h 1 (

 1 sin 2 φ) def =sin 2 φ + 3 sin 4 φ -2 sin 6 φ h 2 (sin 2 φ)

Fig. 2

 2 Fig.2Non-isotropy criterion ρ w.r.t. r and P .

Fig. 3

 3 Fig.3 Approximate and exact ratios CRB(θ)/CRB FF (θ).

Fig. 4

 4 Fig.4Approximate and exact ratios CRB(φ)/CRB FF (φ).

Fig. 5

 5 Fig.5Ratio of the approximate CRB(r) to the exact one.

Note that (1) is no longer valid for parametrizations for which a(α) depends on α.

For P > 6, FIM 1,1 (α) = cσcr 0 sin 2 φ[1-ǫ 2 cos 2 φ+ǫ 4 g 1 (sin 2 φ)+ o(ǫ 5 ) and FIM 1,2 (α) = FIM 1,3 (α) = o(ǫ 5) and the other FIM expressions remain valid.


Houcem Gazzah
email: hgazzah@sharjah.ac.ae

Jean-Pierre Delmas
email: jean-pierre.delmas@it-sudparis.eu

Optimizing V Antenna Arrays Using a Bayesian DOA Estimation Criterion

Experimentations have shown V-shaped uniform antenna arrays to be near-optimum for estimating the direction of arrival of a far-field source, whether the source position is fixed or random. We consider, as performance measure, the expected Cramer-Rao bound, normalized (for comparison purposes) to that of the commonly used uniform circular arrays. We study in details the behavior of V arrays in this context. For large-sized V arrays, the performance measure shows a simple expression, enabling analytical solution of the subsequent array-geometry optimization problem. We obtain closed-form expressions of the orientation, shape and performance of optimal V arrays and learn about their ability to benefit from the available prior about the source direction.

I. INTRODUCTION

Direction finding is a major research field of statistical signal processing that has connections to a wide range of electronic systems and applications, civil and military. The direction of arrival (DOA) of an emitting source is estimated from the signals it induces at a set of sensors that form the antenna array. The most studied scenario is that of sensors displayed in a plane and a source anywhere in the array farfield.

The array geometry, i.e. sensors positions in the plane, has a notorious impact on DOA estimation, at least from two pointsof-view. The first one is ambiguity, i.e. the ability to resolve sources with distinct DOAs. The second one is accuracy, i.e. the ability to distinguish two array output response vectors associated with different sources.

Array ambiguity is a difficult issue, already because there are many definitions of ambiguity: first-order vs higher order (two or more steering vectors are taken into consideration) and strict vs wide sense ambiguity (strict-singularity or nearsingularity of the matrix of steering vectors). By ambiguityfree, it is meant that two steering vectors cannot be colinear if they correspond to two different DOAs. Still, they can be close to be colinear, a situation referred to as wide-sense ambiguity [START_REF] Athley | On radar detection and direction finding using sparse arrays[END_REF], [START_REF] Gavish | Array geometry for ambiguity resolution in direction finding[END_REF]. There are no conditions to strictly avoid wide-sense ambiguity. However, many studies tend to confirm that uniform arrays suffer much less ambiguities [START_REF] Athley | On radar detection and direction finding using sparse arrays[END_REF], from this point-of-view. We adopt a conservative approach of array ambiguity where we focus on uniform antenna arrays, because "For most of the applications, uniform arrays are the first choice . . . " [START_REF] Birinci | Optimization of nonuniform array geometry for DOA estimation with the constraint on gross error probability[END_REF]Sec. I]. For uniform arrays, adjacent sensors are separated by a constant (half-the-wavelength) spacing [START_REF] Birinci | Optimization of nonuniform array geometry for DOA estimation with the constraint on gross error probability[END_REF]. Most notably, they include the Uniform Linear Array used for sources coplanar with the array, and the Uniform Circular Array (UCA) used when sources may appear anywhere in the 3D space. Since we study the latter more general scenario, we will use the UCA as our reference antenna, to which we compare the proposed uniform V arrays.

The second aspect, however, has remained largely unquantified, resulting into a very few theoretical results [START_REF] Baysal | On the geometry of isotropic arrays[END_REF], [START_REF] Hawkes | Effects of sensor placement on acoustic vector-sensor array performance[END_REF] and heuristic-only geometry optimization techniques [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF], [START_REF] Bevelacqua | Optimizing antenna array Geometry for interference suppression[END_REF]. Only recently has a simplification of the Cramer-Rao Bound (CRB) been obtained for the single source case [START_REF] Gazzah | Cramer-Rao bounds for antenna array design[END_REF], telling a lot about the interaction between the array geometry and the estimation accuracy. It has been profitably used to determine the best antenna arrays, for instance, those that minimize the CRB [START_REF] Gazzah | Optimum Ambiguity-Free Directional and Omni-Directional Planar Antenna Arrays for DOA Estimation[END_REF], [START_REF] Gazzah | Optimum Antenna Arrays for Isotropic Direction Finding[END_REF], [START_REF] Gazzah | Direction Finding Antenna Arrays for the Randomly Located Source[END_REF].

The most noticeable outcome of previous optimizations of uniform antenna arrays is that optimum uniform arrays show to have a shape close to V. This fact is true in both deterministic and probabilistic scenarios, i.e. if the source DOA is fixed (and unknown) or random (with known distribution). The V array is also convenient in practice because it is easy to manoeuvre (for instance to adapt the V orientation and shape to the available source prior, as we explain later).

The objective of this paper is to study in details the performance of V arrays in a scenario where the source DOA is not deterministic, but randomly located, following a known probabilistic distribution. It is reasonable to expect the available prior to help improving DOA estimation accuracy. In cellular communication networks, for example, we collect DOA statistics of base station users. The information is used to improve the quality of service and/or offer location-based services [START_REF] Ahonen | Mobile terminal location for UMTS[END_REF]. A variety of modeling techniques [START_REF] Pedersen | A stochastic model of the temporal and azimuthal dispersion seen at the base station in outdoor propagation environments[END_REF], [START_REF] Wong | Landmobile radiowave multipaths DOA-distribution: Assessing geometric models by the open literatures empirical datasets[END_REF] are used for this purpose.

For such a scenario, the CRB is no more a valid performance criterion. W.r.t. the theory of Bayesian estimation, Expected CRB (ECRB) is more relevant, and opportunistically, inherits the simple structure of the CRB. Consequently, analysis conducted in [START_REF] Gazzah | Optimum Ambiguity-Free Directional and Omni-Directional Planar Antenna Arrays for DOA Estimation[END_REF], [START_REF] Gazzah | Optimum Antenna Arrays for Isotropic Direction Finding[END_REF] for fixed DOA, can be extended leading to analytical results which allow rich interpretation and are of great applicability. The first contribution of this paper is to justify the ad-hoc ECRB criterion used in [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF], [START_REF] Gazzah | Direction Finding Antenna Arrays for the Randomly Located Source[END_REF] in the Bayesian framework. The second one is to give a detailed proof of an important result (presented in [START_REF] Gazzah | Direction Finding Antenna Arrays for the Randomly Located Source[END_REF] without justification) about the performance of V arrays. The analytical expressions are followed by a detailed discussion and original illustrations.

The proposed paper is only the second paper (after [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF]) to address array geometry adaptation based on a known distribution of the random source DOA. In comparison to [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF], this paper (i) does not ignore the array ambiguity issue; (ii) does not restrict the source to be in the array plane; (iii) does not restrict the sensors to be in some closed region (iv) includes many analytical results compared to [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF] that contains only the result that, if sensors are restricted to be in a closed connected region, then the optimal sensors placement is along the boundary of the region, as one may expect; (v) shows, for instance, that the array performance depends on the characteristic function at -2, rather than on the complete distribution function; (vi) shows that a major claim [6, Beginning of Sec. V] (that isotropic arrays are optimal for uniformly distributed DOA) is misleading. This is true because the sensors are also required to be inside a disk. As we do not impose such a constraint, we are able to find antenna arrays with better performance.

The paper is organized as follows: First, in Sec. II, we introduce the observation model and recall previous results. We introduce our Bayesian performance criterion in Sec. III, then, apply it to optimize V-structured antenna arrays in Sec. IV. Finally, a conclusion is given in Sec. V.

II. DATA MODEL AND PREVIOUS RESULTS

①

❤ ❤ ❤ ❤ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ z x Source y Θ ρ m Φ φ m
Sensor m Fig. 1. Planar array and source DOAs.

As illustrated by Fig. 1, a planar antenna array is composed of M identical and omni-directional sensors disposed, each, in the (x, y) plane. The position of the m-th sensor is given by its cartesian coordinates x m and y m , or, equivalently, by the complex number γ m =x m + iy m . A source is radiating towards the antenna a narrow-band signal characterized by a wavelength λ. For a source located in the antenna far-field, its DOA is given by a set of two angular parameters: the azimuth angle Φ and the elevation angle Θ, as shown in Fig. 1. Snapshots outputted by the M sensors are collected and used to infer about the azimuth and elevation angles using a variety of algorithms. Some of them are capable of achieving the famous CRB.

Our analysis is conducted under the following assumptions: (i) circular, Gaussian, zero-mean and mutually independent source and noise signals with respective power σ 2 s and σ 2 n ; (ii) independent and identically distributed source snapshots and (iii) noise samples collected from different sensors and/or at different time indexes are mutually independent. Originally available in a very sophisticated expression [START_REF] Porat | Analysis of the asymptotic relative efficiency of the MUSIC algorithm[END_REF], the CRB on (Φ, Θ) has been derived from the inverse of the Fisher Information Matrix (FIM) J F associated with the unknown deterministic parameters (Φ, Θ, σ 2 s , σ 2 n ) in a compact expression that has been greatly simplified (under no additional assumptions) in [START_REF] Gazzah | Cramer-Rao bounds for antenna array design[END_REF] under the following convenient sinusoidal expression:

C Φ,Θ = C ΦΦ C ΦΘ C ΘΦ C ΘΘ , with C ΦΦ = A sin 2 (Θ) B (Φ) C ΘΘ = A cos 2 (Θ) B Φ + π 2 C ΦΘ = A sin(2Θ) ℑ [S 1 exp(-2jΦ)] |S 1 | 2 -S 2 0 where constant A = 1 4π 2 σ 2 n σ 2 s 1 + σ 2 n Mσ 2 s
depends on the observation SNR, the antenna size M and the number of snapshots N ; and

B (Φ) = S 0 + ℜ [S 1 exp(-2jΦ)] S 2 0 -|S 1 | 2
is a function of the source azimuth angle and the arraygeometry dependent constants

S 0 = M m=1 γ m λ 2 - 1 M M m=1 γ m λ 2 , S 1 = M m=1 γ m λ 2 - 1 M M m=1 γ m λ 2 .

III. PERFORMANCE CRITERIA

In practice, there are situations where an a priori knowledge is available about the DOA parameters with (Φ, Θ) following some known joint distribution. The performance of any DOA estimator is measured by the associated MSE matrix [16, (18)

p. 4] R( Φ, Θ) = E( Φ -Φ) 2 E( Φ -Φ)( Θ -Θ) E( Θ -Θ)( Φ -Φ) E( Θ -Θ) 2 .
(1) As the powers σ 2 s and σ 2 n remain in practice unknown deterministic parameters, the theoretical framework to derive lower bound on (1) is the hybrid CRB firstly introduced in [START_REF] Rockah | Array shape calibration using sources in unknown locations. Part I; far filed sources[END_REF], which states that the MSE matrix associated with (Φ, Θ, σ 2 s , σ 2 n ) is lower bounded by the so-called Bayesian CRB given by the inverse of the hybrid information matrix J H = J D + J P [16, (74-75) p.13] with J D = E Φ,Θ (J F ) and J P = J P (Φ, Θ) 0 0 0 , where J P (Φ, Θ) is the prior information matrix w.r.t. the DOA parameters. Consequently we have

R( Φ, Θ) ≥ [(J D + J P ) -1 ] (1:2,1:2)
with

[(J D + J P ) -1 ] (1:2,1:2) ≤ [J -1 D ] (1:2,1:2) ≤ ECRB, (2) where ECRB = E Φ,Θ [J -1 F ] (1:2,1:2) = E Φ,Θ (C Φ,Θ
) is the expectation of the conditional CRB (ECRB) [16, p. 6]. We note that the first inequality in ( 2) is tight for N ≫ 1 if the prior source DOA are not too accurate because in our data model introduced in Sec. II, J F and thus J D is proportional to the number N of snapshots in contrast to J P (Φ, Θ) which depends only on the prior distribution of (Φ, Θ). For example, for azimuth and elevation angles independently Gaussian distributed with prior variance σ 2 Φ and σ 2 Θ , J P (Φ, Θ) =

σ -2 Φ 0 0 σ -2 Θ .
Consequently in practice, the equivalence between J H and J D occurs each time the variance of the DOA Φ and Θ estimated by an efficient algorithm are much smaller than the respective prior variance σ 2 Φ and σ 2 Θ . The second inequality in (2) is derived form the Jensen's inequality applied to the usual FIM J F , which is assumed definite positive, considered here as a random matrix

[E Φ,Θ (J F ] -1 ≤ E Φ,Θ (J -1 F ).
Noting that the MSE given by the maximum a posteriori estimator asymptotically (w.r.t. N ) attains E Φ,Θ (C Φ,Θ ), we propose to base the Bayesian optimization on this bound.

E Φ,Θ (C Φ,Θ ) = CΦΦ CΦΘ CΘΦ CΘΘ .
Furthermore this bound inherits the convenient structure of the CRB displayed in Sec. II. in contrast to the 2 × 2 submatrix [E Φ,Θ (J F ] -1 (1:2,1:2) which has not a simple compact form. Noting that the MSE given by the maximum a posteriori estimator asymptotically (w.r.t. N ) attains E Φ,Θ (C Φ,Θ ), we propose to base the Bayesian optimization on this bound.

E Φ,Θ (C Φ,Θ ) = CΦΦ CΦΘ CΘΦ CΘΘ .
Furthermore, this bound inherits the convenient structure of the CRB displayed in Sec. II, in contrast to matrix

[E Φ,Θ (J F )] -1
which does not have a simple compact form. Finally, we note that the ECRB was used in [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF] as a cost function to be minimized, but only the azimuth DOA angle was considered. We assume, as plausible in practice, that azimuth and elevation are independent. Then, we have

CΦΦ = A E 1 sin 2 (Θ) E [B (Φ)] CΘ,Θ = A E 1 cos 2 (Θ) E B Φ + π 2 ,
where Φ is taken in [-π, π), in order to handle symmetrical distributions.

In the above expressions, terms that do not depend on the array geometry can be disposed of by normalizing to any reference antenna. For convenience, we choose the widelyused UCA, one with M sensors spaced by a distance d.

It verifies B UCA =B (Φ) = 4λ 2 sin 2 (π/M )/(M d 2 ). If we denote ϕ Φ (x) =E [exp(jΦx)]
as the characteristic function relative to the random parameter Φ, we straightforwardly obtain

CΦΦ CΦΦ | UCA = 1 B UCA S 0 + ℜ [S 1 ϕ Φ (-2)] S 2 0 -|S 1 | 2 CΘΘ CΘΘ | UCA = 1 B UCA S 0 -ℜ [S 1 ϕ Φ (-2)] S 2 0 -|S 1 | 2 .
Our criterion is the ECRB normalized to that of the UCA, as expressed above. The latter is a constant, so that minimizing the normalized ECRB is optimal in the general case. A planar array that minimizes this criterion will outperform any antenna array from the ECRB point-of-view.

IV. OPTIMIZATION OF LARGE SIZED V ARRAYS

Array optimization has been conducted in [START_REF] Gazzah | Optimum Ambiguity-Free Directional and Omni-Directional Planar Antenna Arrays for DOA Estimation[END_REF], [START_REF] Gazzah | Optimum Antenna Arrays for Isotropic Direction Finding[END_REF], [START_REF] Gazzah | Direction Finding Antenna Arrays for the Randomly Located Source[END_REF] for curve-like uniform arrays formed of a sequence of regularly spaced sensors. Thanks to an adequate parametrization of the sensors positions using angular (so bounded) parameters, array optimization can be achieved by exhaustive search. For a number of M sensors, the total number of search points is extremely huge. Likely enough, the results of the exhaustive search has shown optimal array to have a shape close to V. Conveniently enough, V array geometry is function of the angles characterizing the orientations of the two branches of the V shape and of the length of the branches, i.e. M . The normalized ECRBs will depend on these three parameters only and optimization will be much easier. For instance, it can be achieved analytically. Furthermore, such optimization will be found to be independent from M .

As shown in Fig. 2, the V array is made of an odd number M of sensors. Sensor 1 is placed at the origin while, for k = 1, • • • , (M -1)/2, sensors 2k and 2k + 1 are positioned such that γ 2k = kd exp(j∆ 1 ) and γ 2k+1 = kd exp(j∆ 2 ). When M tends to infinity, it has been shown that [START_REF] Gazzah | Direction Finding Antenna Arrays for the Randomly Located Source[END_REF] the normalized ECRBs become independent from M . For instance, they are equal to

r ǫ (∆ 1 , ∆ 2 ) = 3 2π 2 5 -3 cos(∆ 1 -∆ 2 ) +ǫ [5 cos(∆ 1 -∆ 2 ) -3] ℜ [exp[j(∆ 1 + ∆ 2 )]ϕ Φ (-2)] sin 2 (∆ 2 -∆ 1 )
.

where ǫ = 1 (resp. -1) for the azimuth angle Φ (resp. elevation angle Θ). Our objective is to study in details functions r 1 (∆ 1 , ∆ 2 ) and r -1 (∆ 1 , ∆ 2 ), and specifically, angles ∆ 1 and ∆ 2 for which the V array achieves a minimal r 1 (∆ 1 , ∆ 2 ) or r -1 (∆ 1 , ∆ 2 ).

A. ECRB-Minimizing V arrays

We prove in Appendix that min ∆1,∆2 r ǫ (∆

1 , ∆ 2 ) = 3 4π 2 (3 -5α) 2 5 -3α -4 √ 1 -α 2 , if α = 3 5 (3) = 24 5π 2 , if α = 3 5 . ( 4
)
where the key-parameter α =|ϕ Φ (2)| expresses the amount of prior available about the source azimuth angle. The proof of ( 3) is detailed in Appendix, while ( 4) is obtained by continuity of (3).

In Appendix, we also prove that

r ǫ (∆ † 1 , ∆ † 2 ) = min ∆1,∆2 r ǫ (∆ 1 , ∆ 2 ) iff cos(∆ † 1 -∆ † 2 ) = 4 √ 1 -α 2 + 3α -5 5α -3 (5) cos(∆ † 1 + ∆ † 2 ) = ǫ E [cos (2Φ)] α (6) sin(∆ † 1 + ∆ † 2 ) = ǫ E [sin (2Φ)] α . (7)
W.r.t. the other cost function r -ǫ (∆ 1 , ∆ 2 ), the V array characterized by the so-defined ∆ † 1 and ∆ † 2 achieves

r -ǫ (∆ † 1 , ∆ † 2 ) = 3 16π 2 5α -3 √ 1 -α 2 × 3 + 5α + (5 + 3α)(3 -5α) 3α -5 + 4 √ 1 -α 2 .
The proof of the above follows immediately. Proof Using the same notations as Appendix, and bearing in mind that α ǫ,u = α, and following steps similar to those that led to [START_REF] Gazzah | Optimum Ambiguity-Free Directional and Omni-Directional Planar Antenna Arrays for DOA Estimation[END_REF], we can prove that

r -ǫ (δ 3 , δ 4 ) = 3 2π 2 5 + 3α -(5α + 3) cos(δ 3 ) 1 -cos 2 (δ 3 ) .
Implementing ( 13) and ( 14), where we denote β =β ǫ,u , we get

2π 2 3 r -ǫ (δ 3 , δ 4 ) = 5+3α 2 √ β 2 -1 - √ β 2 -1+sign(5α-3)β + 3+5α 2sign(5α-3) √ β 2 -1 . As a consequence, sign (5α -3) β 2 -1 4π 2 3 r -ǫ (δ 3 , δ 4 ) = 5 + 3α β -sign (5α -3) β 2 -1 + 3 + 5α, i.e. sign (5α -3) √ 1 -α 2 |5α -3| 16π 2 3 r -ǫ (δ 3 , δ 4 ) = 5 + 3α β + sign (3 -5α) 4 √ 1-α 2 |3-5α|
+ 3 + 5α, and finally

√ 1 -α 2 5α -3 16π 2 3 r -ǫ (δ 3 , δ 4 ) = 5 + 3α 3α-5 3-5α + 4 √ 1-α 2 3-5α + 3 + 5α = (5 + 3α)(3 -5α) 3α -5 + 4 √ 1 -α 2 + 3 + 5α

B. Interpretation

First, notice that the angle ∆ † 1 -∆ † 2 between the two branches of the optimal V array does not depend on the angle of interest (i.e. not on ǫ). It is shown in Fig. 3 as function of the available azimuth prior, as expressed by α. The angle (parameter) of interest (or also ǫ) affects the orientation of the (optimal) V array, which is not unique by the way. To identify V arrays that satisfy (5)-( 7), we let ∆ 3 = arccos (4 √ 1 -α 2 + 3α -5)/(5α -3) and ∆ 4 be the angle in [0, 2π[ whose cosine and sine are given by the right-hand sides in ( 6)- [START_REF] Bevelacqua | Optimizing antenna array Geometry for interference suppression[END_REF], respectively. Let's also define

∆ 5 =(∆ 4 -∆ 3 )/2 and ∆ 6 =(∆ 3 + ∆ 4 )/2.
On one hand, ∆ 4 is, by definition, in [0, π). Given that ∆ 1 + ∆ 2 is in [0, 4π), the system with ( 6) and ( 7) has two possible solutions for

∆ 1 + ∆ 2 . If ǫE [sin (2Φ)] is positive, then ∆ 1 + ∆ 2 = ∆ 4 . Otherwise, ∆ 1 + ∆ 2 = ∆ 4 + 2π. On the other hand, ∆ 3 is, by definition, in [0, π). Given that ∆ 1 -∆ 2 is in [-2π, 2π
), based on (6), there are 4 possible values for

∆ 1 -∆ 2 . It is equal to either ∆ 3 , -∆ 3 , ∆ 3 -2π or 2π -∆ 3 .
Hence, we have 8 possible cases for the tuple (∆ 1 , ∆ 2 ) which are: (∆ 6 , ∆ 5 ), (∆ 5 , ∆ 6 ), (∆ 6 -π, ∆ 5 + π), (∆ 5 + π, ∆ 6 -π), (∆ 6 + π, ∆ 5 + π), (∆ 5 + π, ∆ 6 + π), (∆ 6 , ∆ 5 + 2π) and (∆ 5 + 2π, ∆ 6 ). By definition, ∆ 5 and ∆ 6 are in [-π, π) and [0, 2π), respectively. Hence, some of these solutions correspond to identical V arrays. Actually, among these V arrays, only two are different. They are (∆ 5 , ∆ 6 ) and (∆ 5 + π, ∆ 6 + π). These two V arrays are obtained one from the other by symmetry w.r.t. to the origin.

C. Distributions of

Φ such that |ϕ Φ (2)| = 1
The largest reduction of the ECRB is obtained for |ϕ Φ (2)| = 1. This is the case of a deterministic prior, but not only. More generally, we prove that α is equal to 1 iff the random azimuth has two possible values spaced by π, with arbitrary probabilities (hence, including the deterministic case). This is not surprising since the CRB is identical for two such directions. To establish the proof, let's, first, prove the following lemma. Lemma Let I be a bounded interval and f (x) be a summable complex-valued function over I. Then, Hence, the argument of f (Φ) exp(2jΦ), i.e., 2jΦ is a constant. Let's denote it Φ 0 , chosen in [-π, π). Then, Φ is equal to any number of the form Φ 0 /2 + kπ, where k is any integer. There are two such numbers in [-π, π)

V. CONCLUSION

We study in details the performance of uniform V arrays w.r.t. the estimation of the azimuth and elevation angles of a randomly located source. The expected CRB, normalized to that of the UCA, is our performance measure. It depends on the direction of the source (in fact, of the only azimuth angle) and the orientation of the branches of the V shape. Analytical expressions are obtained and closed-form solution of the optimal V array: its shape, orientation and accuracy.

APPENDIX

We denote δ 3 =∆ 2 -∆ 1 and δ 4 =∆ 1 + ∆ 2 and rewrite

r ǫ (δ 3 , δ 4 ) = 3 2π 2 5 -3 cos(δ 3 ) + ǫ [5 cos(δ 3 ) -3] × {cos (δ 4 ) E [cos (2Φ)] + sin (δ 4 ) E [sin (2Φ)]} sin 2 (δ 3 ) .
Minimization of r ǫ w.r.t. δ 4 straightforwardly leads to

tan(δ 4 ) = E [sin (2Φ)] /E [cos (2Φ)] . (8)
Thus, cos(δ 4 ) = uE [cos (2Φ)] /α, where u = ±1 and

α = {E [sin (2Φ)]} 2 + {E [cos (2Φ)]} 2 , i.e., |ϕ Φ (-2)| = |ϕ Φ (2)
|, which, by definition, is less than one. Replacing δ 4 by its expression, we rewrite r ǫ as follows

r ǫ (δ 3 , δ 4 ) = 3 2π 2 5 -3 cos(δ 3 ) + α ǫ,u [5 cos(δ 3 ) -3] sin 2 (δ 3 ) ,
where α ǫ,u =uǫα. This is a function of cos(δ 3 ),

r ǫ (δ 3 , δ 4 ) = 3 2π 2 5 -3α ǫ,u + (5α ǫ,u -3) cos(δ 3 ) 1 -cos 2 (δ 3 ) . (9)
The minimization of r ǫ w.r.t. cos(δ 3 ) leads to the second order equation cos 2 (δ 3 ) + 2β ǫ,u cos(δ 3 ) + 1 = 0, where

β ǫ,u =(5 -3α ǫ,u )/(5α ǫ,u -3). (10)
Notice that

β 2 ǫ,u -1 = 16(1 -α 2 )/((3 -5α ǫ,u ) 2 ) (11)
is always positive because α is in [0, 1], by definition. Hence, cos(δ 3 ) is equal to the one among -β ǫ,u + β 2 ǫ,u -1 and -β ǫ,u -β 2 ǫ,u -1 that has an absolute value lower than 1. Notice that the product of these two values is 1, so one and only one is a valid cosine. Now, one can verify from (10) that β ǫ,u -1 = 8 (1 -α ǫ,u ) / (5α ǫ,u -3) and β ǫ,u + 1 = 8 (1 + α ǫ,u ) / (5α ǫ,u -3). Both have the same sign, since their respective numerators are positive.

On one hand, if α ǫ,u ≥ 3/5, then β ǫ,u -1 and β ǫ,u + 1 are both positive. Hence, we have

-β ǫ,u ≤ -1 ≤ β ǫ,u , implying that (β ǫ,u -1) 2 ≤ β 2 ǫ,u -1 ≤ (β ǫ,u +1) 2 , or also that β ǫ,u -1 ≤ β 2 ǫ,u -1 ≤ β ǫ,u + 1. Hence, cos(δ 3 ) = -β ǫ,u + β 2 ǫ,u -1.
On the other hand, if α ǫ,u ≤ 3/5, then β ǫ,u -1 and β ǫ,u + 1 are both negative. Hence, we have

β ǫ,u ≤ -1 ≤ -β ǫ,u , implying that (β ǫ,u + 1) 2 ≤ β 2 ǫ,u -1 ≤ (β ǫ,u -1) 2 , or also that -β ǫ,u + 1 ≤ β 2 ǫ,u -1 ≤ -β ǫ,u -1. Hence, cos(δ 3 ) = -β ǫ,u -β 2 ǫ,u -1. As a conclusion, cos(δ 3 ) = sign (5α ǫ,u -3) β 2 ǫ,u -1 -β ǫ,u . (12)
In order to update (9), we write that

cos 2 (δ 3 ) = 2β 2 ǫ,u -1 -2sign (5α ǫ,u -3) β ǫ,u β 2 ǫ,u -1 sin 2 (δ 3 ) 2 = 1 -β 2 ǫ,u +sign (5α ǫ,u -3) β ǫ,u β 2 ǫ,u -1 (13)
= sign (5α ǫ,u -3)

β 2 ǫ,u -1 × -sign (5α ǫ,u -3) β 2 ǫ,u -1 + β ǫ,u
= -sign (5α ǫ,u -3) β 2 ǫ,u -1 cos(δ 3 ),

so that we rewrite (9) as 2π .

To determine whether the minimum of r ǫ is met with u equal to 1 or -1, we study the sign of r ǫ (δ 3 , δ 4 )| u=1r ǫ (δ 3 , δ 4 )| u=-1 , which, w.r.t. the above, is also that of (3 -5α ǫ,1 ) As a conclusion, r 1 (resp. r -1 ) is minimized by δ 3 and δ 4 for which u = 1 (resp. u = -1). Hence u = ǫ and α ǫ,u = α.

For such δ 3 and δ 4 minimizing r ǫ we have, on one hand, by application of [START_REF] Ahonen | Mobile terminal location for UMTS[END_REF], cos(δ 3 ) = sign (5α -3) β 2 ǫ,ǫ -1 -β ǫ,ǫ . By means of ( 10) and [START_REF] Gazzah | Direction Finding Antenna Arrays for the Randomly Located Source[END_REF], this can be proved to be equal to [START_REF] Hawkes | Effects of sensor placement on acoustic vector-sensor array performance[END_REF]. On the other hand, cos(δ 4 ) = ǫE [cos (2Φ)] /α, in addition to (8) lead to [START_REF] Oktel | A Bayesian approach to array geometry design[END_REF] and [START_REF] Bevelacqua | Optimizing antenna array Geometry for interference suppression[END_REF].

The minimum achieved by r ǫ is, independently from ǫ,

2 Fig. 2 .

 22 Fig. 2. The V-shaped antenna array.

Fig. 3 .

 3 Fig.3. The angle between the two branches of the optimal V array as function of parameter α.

I

  f (x)dx = I |f (x)| dx iff arg [f (x)] is constant over I almost everywhere. Proof Consider the constant z = I f (x)dx. Then I f (x)dx = e -j arg(z) z = ℜ e -j arg(z) z = I ℜ e -j arg(z) f (x) dx that is meant, here, to be equal to I |f (x)| dx. However, ℜ e -j arg(z) f (x) is lower than its magnitude, i.e., lower than |f (x)|. Both being real-valued functions, their respective integrals over I can be equal iff the two functions are identical almost everywhere over I. By writing that ℜ e -j arg(z) f (x) = |f (x)|ℜ e j arg[f (x)]-j arg(z) = |f (x)| cos {arg[f (x)] -arg(z)} must be equal to |f (x)|, we conclude that arg[f (x)] -arg(z) must be a constant (0 modulo 2π), where arg(z) is itself a constant. This terminates the proof. In order to apply the above lemma, we rewrite |ϕ Φ (2)| = 1 as π -π f (Φ) exp(2jΦ)dΦ = π -π |f (Φ) exp(2jΦ)|dΦ = 1, where f (Φ) is the PDF associated with the random variable Φ.

+ 3 - 1 . 5 - 4 √ 1 -α 2 + 5 - 2 3 2 3

 315412522 2 3 r ǫ (δ 3 , δ 4 ) = 5 -3α ǫ,u 2 β 2 ǫ,u -1 -β 2 ǫ,u -1 + sign (5α ǫ,u -3) β ǫ,u 5α ǫ,u 2sign (5α ǫ,u -3) β 2 ǫ,u -Using (10) and (11), we simplify the above as follows3α ǫ,u ) |5α ǫ,u -3| -3α ǫ,u -|5α ǫ,u -3| , r ǫ (δ 3 , δ 4 ) = 5-3αǫ,u -4 √ 1-α 2 +5-3αǫ,u -1 or also √ β 2 ǫ,u -1 |5αǫ,u-3| π r ǫ (δ 3 , δ 4 ) =

2

 2 which can be found to be equal toǫα 150α 2 + 240 √ 1 -α 2 -246 25 -9α 2 -40 √ 1 -α 2 + 16(1 -α 2 ) = -6ǫα.

1 ,

 1 which can be verified to be equal to (3).