Chinese Journal of Aeronautics 25 (2012) 83-93

Contents lists available at ScienceDirect

Chinese Journal of Aeronautics

journal homepage: www.elsevier.com/locate/cja

JOURNAL OF

AERONAUTICS

A Weighted Inner Product Estimator in the Geometric Algebra of Euclidean 3-Space for Source Localization Using an EM Vector-sensor

JIANG Jingfei, ZHANG Jianqiu*

Department of Electronic Engineering, Fudan University, Shanghai 200433, China Received 21 January 2011; revissed 3 March 2011; accepted 28 June 2011

Abstract

In this paper, the source localization by utilizing the measurements of a single electromagnetic (EM) vector-sensor is investigated in the framework of the geometric algebra of Euclidean 3-space. In order to describe the orthogonality among the electric and magnetic measurements, two multivectors of the geometric algebra of Euclidean 3-space (G3) are used to model the outputs of a spatially collocated EM vector-sensor. Two estimators for the wave propagation vector estimation are then formulated by the inner product between a vector and a bivector in the G3. Since the information used by the two estimators is different, a weighted inner product estimator is then proposed to fuse the two estimators together in the sense of the minimum mean square error (MMSE). Analytical results show that the statistical performances of the weighted inner product estimator are always better than its traditional cross product counterpart. The efficacy of the weighted inner product estimator and the correctness of the analytical predictions are demonstrated by simulation results.

Keywords: cross product; electromagnetic; geometric algebra; geometric algebra of Euclidean 3-space; minimum mean square rule; direction finding; vector-sensor

1. Introduction

The cross product algorithm, based on the fact that the instantaneous electric and magnetic vectors of an electromagnetic (EM) plane wave and the direction vector of wave propagation are mutually orthogonal, has been widely applied to such areas as EM source tracking [1-2] and parameter estimations [3-10] in recent years. The reasons for its wide applications are at least threefold. Firstly, since the orthogonal relationships do not depend on the frequency of the EM waves, the cross product algorithm can be equally applied to both

♦Corresponding author. Tel.: +86-21-55664226. E-mail address: jqzhang@fudan.ac.cn

Foundation items: National Natural Science Foundation of China

(61171127); National Basic Research Program of China (2011CB302903)

1000-9361/$ - see front matter © 2012 Elsevier Ltd. All rights reserved. doi: 10.1016/S1000-9361(11)60365-8

the narrowband and wideband sources with very low computational complexity [1]. Secondly, the cross product algorithm is based on the measurements of a single EM vector-sensor, which consists of three orthogonally oriented electric antennas and three orthogonally oriented magnetic antennas. On the one hand, when the six antennas are spatially collocated in a point-like geometry [3-14], there is no need for sensor position calibration and time synchronization since only the instantaneous measurements of an uni-vector-sensor are used. On the other hand, if the six antennas are not spatially collocated but elaborately placed, for example using the ingenious and simple scheme presented in Ref. [15], the mutual coupling among the antennas can be effectively minimized in that the spacing between the antennas can be quite large. Thirdly, the cross product algorithm can be incorporated into the traditional algorithm of direction of arrival (DO A) estimation, which is based on the spatial

phase delay across a vector-sensor array. Creative synergy between these two algorithms has produced several advantages for high resolution vector-sensor array processing, including the improvement of DOA estimation accuracy [4-6], the DOA estimation without knowing the array geometry [7-8], self-initiates the multiple signal classification (MUSIC) iteration [9], array aperture extension [10], etc. As a result, if the performances of the cross product estimator are improved, better results of the aforementioned applications and algorithms can be achieved.

In this paper, we concentrate on improving the performances of the traditional cross product estimator based on the measurements of a single EM vector-sensor. A weighted inner product estimator, aimed at estimating the propagation direction vector of the incident EM waves, is proposed based on the minimum mean square error (MMSE) fusion rule. To evaluate the performance of the weighted inner product estimator, an analytical comparison to its cross product counterpart is drawn. From the comparison results, we proved that our weighted inner product estimator is always superior to the traditional cross product estimator. However, the price for the performance improvement is that extra computational efforts are required for the MMSE fusion process. In addition, theoretical analyses also reveal that both the cross product and weighted inner product estimators are biased in the presence of mutual coupling. Hence, to make our work meaningful, a simple strategy is proposed to calibrate the undesired coupling before applying our estimators to the sensor outputs.

2. Geometric Algebra of Euclidean 3-Space (G3)

Before starting our discussion, we would like to add a comment on the notation that is used. To facilitate the distinction between scalars, vectors, multivectors, and the ^-grade-vector, the type of given quantity will be reflected by its representation: scalars (including real numbers and complex numbers) are denoted by lower-case letters (a, b, •••), vectors are written as a lower-case letter with bold face (a, b, •••), multivec-tors are expressed by capitals (A, B, ••• ) (italic shaped), the ^-grade-vector part of a multivector A is represented as <A>k. A special multivector in G3 is the one with only the scalar and trivector parts. Such a multivector is mathematically isomorphic to a complex number [16-18], so we use a lower-case letter to represent it for notational simplicity in this paper.

To facilitate the subsequent discussions, we list here some basic properties of the G3 necessary for our purpose in this paper. For the interested readers, a thorough review about the geometric algebra can be found in Refs. [16]-[18], while some properties of the geometric algebra and its applications to physics were discussed in Ref. [19]. In the signal processing field, the application of the geometric algebra to image processing was introduced in Ref. [20], the techniques

of utilizing the geometric algebra to represent power under non-sinusoidal conditions were reported in Ref. [21], and the Fourier transform in the geometric algebra domain was investigated in Ref. [22].

2.1. Basics

The geometric algebra of G3, noted as £3, is an 8D algebra system and consists of scalars (0-grade-vecors), vectors (1-grade-vectors), bivectors (2-grade-vectors), and trivectors (3-grade-vectors) [17]. A generic element of the G3 is nominated as a multivector, which is a mathematical object of "mixed" dimensions and can be expressed by

A — a0 + a1e12 + a2 e31 + a3e23 + a4 ei 23 l a5 ei I a6 e2 I a7 e3 —

<A>o +<A>2 +<A>3 +<A)i e £3 (1)

where a0, a1, •••, a7 are real numbers; (A)k (k=0, 1, 2, 3) is the k-grade-vector part of A; {1}, {e1, e2, e3}, {e12, e31, e23}, and {e123} are respectively the basis elements of the scalar, vector, bivector, and trivector parts of the G3. The multiplication rules of the eight bases are summarized in Table 1. Since the multiplication obeys the left and right distributive rules with respect to addition [18], the product of any two multivectors in the G3, defined as the geometric product, can be obtained by using the multiplication rules depicted in Table 1.

Table 1 Multiplication rules of geometric algebra of G3

1 e1 e2 e3 e12 e31 e23 e123

1 1 e1 e2 e3 e12 e31 e23 e123

e1 ei 1 e12 -e31 e2 -e3 e123 e23

e2 e2 -e12 1 e23 -e1 e123 e3 e31

e3 e3 e31 -e23 1 e123 e1 -e2 e12

e\2 e12 -e2 e1 e123 -1 e23 -e31 -e3

e31 e31 e3 e123 -e1 -e1 -1 e12 -e2

e23 e23 e123 -e3 e2 e31 -e12 -1 -e1

e123 e123 e23 e31 e12 -e3 -e2 -e1 -1

Definition 1 The magnitude of a multivector A e £3 is a unique scalar calculated by [22]

1142 =< AA > 0 =£« A) k < At > k > 0

A — ^0 2 a 2 ^31 a é*23 ^4^123 + tf5e I a6^2 I a7^3 — <A)0 -<A)2 -<A)3 +<A>1 e G3

is the reverse of A and

means the magnitude of a multivector .

From Table 1 and Definition 1, it can be seen that e1223 — -1 and e/23 — -e123, which are similar to the properties usually attributed to the complex imaginary

unit j because j =-1 and the conjugation of j is -j. In fact, it can be proved that e123 is mathematically isomorphic to j, since {1, e123} constructs an algebra system isomorphic to the complex algebra [18].

Definition 2 The inner product of two vectors can be extended to two multivectors in the G3. The inner product of a vector a and a ^-grade-vector <B>k is defined by [18]

a-<B)k = 2(a<B)k - (-1)k<B)ka)

where a (B)k and (B)k a are the left and right geometric product of a and (B)k (the geometric product is non-commutative); the geometric interpretation of a• {B)k is that it is a (k-l)-grade-vector representing the complement (within the subspace defined by (B)k ) of the orthogonal projection of a onto (B)k [19]. In this way, the inner product by a vector lowers the grade of any simple multivector by one.

Definition 3 For a vector a and any multivector (B)k of grade k, their outer product is defined by [17]

a a <B)k = 2 (a<B)k +(-1)k <B)k a)

The geometric interpretation of a a (B)k is that it is a (k+1)-grade-vector spanned by the k-grade-vector (B)k and the vector a. As a result, the outer product by a vector raises the grade of any simple multivector by one.

2.2. Relations with vector algebra

The vector algebra developed by J. Willard Gibbs in 1884 fits naturally into the G3 [18]. For this purpose, we just give the relations between the outer product and the cross/inner product introduced by Gibbs.

Let a and b be any two vectors in the G3, their cross product a*b can be expressed by

a a b = e123a x b

or equivalently

a x b = -a-(be123) (7)

The inner product a • b can also be represented by

a-b = -e123(a a (e123b)) (8)

3. G3 Description of an Uni-vector-sensor

Suppose that a single EM vector-sensor locates at the origin of the coordinated system and there is a narrow-band source, parameterized by &=[$ 6 y rj\, radiating EM waves to the vector-sensor. When the medium is nonconductive, homogenous and isotropic, the noise free signals that would be received at the

" s Ex (». t ) " - sin^ cosjcos

SEy (ß, t) cos^ sin j cos

sEz (ß t) 0 - sin0

Shx (ß t ) - cos tj> cos 0 - sin j

SHy (ß t ) - sin j cos 0 cos j

_Sh (ß t)_ sin 0 0

ee,23n sin y

s(t ) = a(0) s(t)

ux sin 0 cos j

uy = sin 0 sin j

A _ cos0

sensor are

[3-10]

where O<0<n represents the source's elevation angle measured from the vertical z axis, 0<$<2n the source's azimuth angle, 0<y<n/2 the auxiliary polarization angle, -n<rj<n the polarization phase difference angle, a(&) the well-known steering vector, and p(y, n) the polarization parameter of the source [23]. In Eq. (9), the electric-field vector [£&(/) sEy(t) sEz(t)]T and the magnetic-field vector [sh*(0 sHy(t) sHz(t)] T are respectively measured by the electric and magnetic antennas of the vector-sensor. Based on the definitions of elevation angle 0 and azimuth angle the direction vector of the source, which is opposite to the wave propagation vector, is analytically expressed as [3-10]

where ux, uy, and uz are the well-known direction cosines of the source along x axis, y axis, and z axis.

For notational simplicity and without loss of generality, the explicit dependence on & will be omitted in the following discussions. s(t)= ||s(0|| esv is the

complex envelope of the source's signal with (p being the phase angle. Since e123 is isomorphic to the complex imaginary unit j [18], we use e123 to replace j in Eq. (9) to facilitate our discussions in the framework of the G3.

From the standpoint of G3, the original received data of each vector-sensor are six multivectors (i.e., SEx, SEy, SEz, SHx, SHy, and SHz) with only scalar and trivector parts. Taking the orthogonality among the electric or magnetic measurements into consideration, we can respectively express them as

Se(0 = Sex (t)ei + SEy (t )e2 + SEz (t^ (11)

Sh(0 = Shx (Oe, + SHy (/>2 + ShZ (/>3 (12)

Based on the rules in Table 1, it is explicit that SE(t) and SH(t) are multivectors consisting the vector and the bivector parts. Compared with Eq. (9), Eqs. (11)-(12) are more coincident with the physical nature of the EM waves, for they can vividly illustrate the orthogonality among the electric and magnetic signals, while the long vector type modeling in Eq. (9) cannot. Similar to Eqs. (11)-(12), the direction vector of the source can be expressed in the following G3 form:

u — e1 sin0cos^ + e2 sin0sin^ + e3 cos0 (13)

In applications, the sensor measurements are always corrupted by noise. Then the real data of the array outputs can be modeled as

YB(t ) = SE(t ) NB(t ) YH(t ) = SH(t ) NH(t )

where Ne(0 = n^t) e1 + nEy (t) e2 + «ez (t) e3 and Nh (t) = nHx (t)e1+nHy(t)e2+nHz(t)e3 are the G3 formulations of the electric and magnetic measurement noise respectively. In general, the noise components are modeled as complex vectors [3]. Correspondingly, they can be represented in the G3 framework as

NE(t ) = ((t ) + (t ) )ei +

(nEyr (t) + el23nEyi (t ) )e2 + (nEzr (t) + ei23nEzl (t))e3 =

«Er(t ) +

ei23nEi (t) NH(t ) = (( ) + ei23nHxi (t ) )ei + (( ) + ei23nHyi (t ) ) + (nHzr (t) + el23nHzl (t))e3 =

nHr(t ) +

ei23WHi(t )

where nEr, nEi, nHr, and nHi are (real) vectors.

Similar to the complex domain formulation [3], some necessary assumptions on the signals and noise components in the G3 framework are listed as follows.

Assumption 1 Equivalent to the complex formulation [3], the source signal sequence {s(1), s(2), • • •} is a sample from a temporally uncorrelated stationary (complex) random process with zero mean and

E {s (p) (s(q) )t| — o28 p,q, E {s (p)s(q)} — 0

where s(t)=sr(t)+e123si(t) with sr(t) and si(t) being two real random variables, a2 is an unknown positive real number, 8pq the Kronecker delta, and E{ •} the expectation operator.

Assumption 2 In Eq. (16), the measurement noise of the electric signals is a three-dimensional complex random process with zero mean and

E {{(t) (NE(t) )t} — 3^8 p, q, E {NE(t) (NE(t))} — 0 where

E {(() + e123 nHmi (t) )( Hmr (t) + e123nHmi (t) )} —°H

with m=x or y or z.

Assumption 4 The electric noise is independent of the magnetic noise. The source signals and the measurement noise are uncorrelated.

4. Weighted Inner Product Estimator

This section presents a simple algorithm to estimate the direction vector u of the EM source via the measurements of a uni-vector-sensor. The idea behind this estimator is the fact that the mutually orthogonality among the electric vector, magnetic vector, and the wave propagation vector of an EM wave can be formulated by the inner product between a vector and bivector in the G3, as described below.

(16) 4.1. Two independent inner product estimators

Before rendering our discussions, we would like to make the following definitions:

v1 —-e1sin^ + e2cos^ (18)

v2 — -e1 cos^cosû-e2 sin^cos0 + e3 sin0 (19)

According to Eq. (13), Eq. (18) and Eq. (19), it can be checked that u±v1±v2 and (u, v1, v2) form a right orthogonal triad. By further defining that />1=cos y, p2r=sin y cos n, p2i = sin y sin n, sr(t) =||s(t)||cos <p, and

s(t) = |\s(t)|| sin <p, we then have

E {(((t) + e123 + nE mi (t) )(() + e123 nE mi (t) )t}

with m=x or y or z.

Assumption 3 The assumption on the magnetic noise components is similar to that of the electric ones, i.e.,

E { Nh (t) (NH(t) )t} — 3oH 8 p ,q, E {t) NH(t)} — 0

P(Y, n) — [P1 P2r + e123P2i ] (20)

s(t ) — Sr(t) + el23si(t) (21)

Substituting Eqs. (20)-(21) into Eq. (9), Eqs. (18)-(19) imply that Eqs. (11)-(12) can be expressed as

SE (t) — P1sr (t)v1 - P2r sr (t)v2 + Pi,si (t)v2 +

e123 (( (t)v1 - Pirsi (t)v2 - Pi,sr (t)v2 ) — x(t) + en3 y(t )

e123 SH (t ) = uSE (t) = ux (t) + e123uy(t)

x(t) = P1Sr(t)V1 - P2t Sr(t)V2 + P2iSi(t)V2

y (t) = ASi(t)V1 - P2rSi(t)V2 + P2iSr(t)V2

are two vectors. Since u_v1_v2, Eqs. (24)-(25) mean that x(t) _L u and j(t) _L u. Consequently, Eq. (23) equals

e123SH (t) — uSE (t) — u a x(t) + e123u a j(t) (26)

Eqs. (24)-(25) also confirm that SE(t) and e123SH(t) are multivectors with only the vector and bivector parts. Based on such an observation and the inner product definition in Eq. (4), we obtain the following two relations from Eq. (22) and Eq. (26), i.e.,

< Se(/)>1-<e123Sn (t )>2 = x(t)-(u a x(t)) = -(x(t) • x(t))u (27)

<e123Sh(/)>•< Se(/)> 2 = em (u a y (t)) • (ej23y(t)) = - (y(t) ■ y(t)) u (28)

Since the inner products x(t) • x(t) and y(t) • y(t) are scalars, the results in Eqs. (27)-(28) are two vectors parallel to u, which means, after a normalization process, u can be obtain from Eqs. (27)-(28) by using the noise free sensor data. However, the real sensor data are always corrupted by measurement noise as represented in Eqs. (16)-(17). Using the noise representations in Eqs. (14)-(15), Eqs. (27)-(28) will change to

<7e(/)>! •<e123^H(t)>2 = [<Se(/)>! + %r(t)H<e123 Sh(/)>2 + e123«Hr(t)j =

- (x(t) • x(t)) u + ej23 (x(t) a «Hr (t)) +

e123 («Er (t) a «Hr (t)) + «Er (t) • (u a x(t)) (29)

<e!23^(t )>1 -<^e(/)> 2 =

[<e123SH(t)>1 - «Hi(t)][^^E(t)> 2 + e123«Ei(t)] =

- (y (t) • y (t)) u + em (y (t) a «Hi (t)) +

e123 ((t) a «Hi (t)) + «e (t) • (u a y(t)) (30)

where the last three items on the right sides are the noise perturbation components for the inner product results. If an average process is used to reduce the perturbation, we can get two estimators of the direction vector u as

Wx = KK i [ ))i-<e,23^n(t )) 2 ] K t=,

Wy =1 it [MOV^t)) 2 ] k 1

ù y =■

, is the l2-norm of a vector.

Theorem 1 Under Assumption 1 to Assumption 4,

it is almost sure that ù x

> u and ù,,

Proof See Appendix A.

4.2. An optimal weighted inner product estimator

As shown in Eqs. (29)-(30), the noise components in the two estimators ux and uy are independent of

each other. Motivated by such an observation, the idea of fusing ux and uy into a better estimator comes natu-

really. In this subsection, we introduce a weighted inner product estimator to fulfill this idea.

Rewrite ùx and ùy as ùx = i üxiet and ùy =

^ m and «v a s ) u, it is easy

i uy,e,. Since ù i=1

to check that

= i Uoiei =i[®,ux, + (1 -®i)u y, ] <

i=1 i=1 where 0<®,<1. Let the variances of U„, ü.

ù (33)

, and Uoi

beo-2, cr2yi, and ct2 respectively. From Eq. (33), one has

^o2, + (1 -®r fei

The method we applied to obtain coi is to minimize Eq. (34), which is known as the MMSE fusion rule for a multi-sensor system [24]. Thus, one has

' 1 1 ^

and the minimum variance of Uoi is

al =■

a2" +"a2

On ay,

According to Eq. (36), it can be easily checked that < << and < <<2. As a result, if < and <2

can be obtained in advance, the optimal weight factor c can be calculated to yield the MMSE estimation of u.

In applications, the on-line optimum estimation of u requires the on-line estimation of < and <2yi respectively. Let the mth estimation of u using Eqs. (29)-3 3

(3°) be ux(m) = £ Uxi(m)e,and uy(m) = £ ayt (m)e .

i=1 i=1

When the data length for the on-line estimation of < and <r2 is M, their variances can be estimated as

= — I

uxi(m) - m i uxi (q)

M q=m-M+1

u,i(m) - mm i u,i (q)

q=m-M+1

Equations (37)-(38) can be considered as two moving window estimators with a window length M. These windows keep moving forward with the new data entering. In each time, what are seen from the windows are the most recent M estimates given by ux and uy,

and what are produced by the moving window estimators are the current variances for all the components of u x and u,,.

ùx =-

In summary, our weighted inner product estimator can be implemented by the steps as follows.

Step 1 Select the number of snapshots K.

Step 2 Choose the window length M for the variance estimations.

Step 3 Use the most recent K snapshots of sensor data to get the newest estimations of ux and uy by Eqs. (31)-(32).

Step 4 Estimate ct2 and ct2 by Eqs. (37)-(38).

Step 5 Obtain (i=1,2,3) by substituting the estimated a2 and a2 in Eq. (35).

Step 6 Calculate u0 by Eq. (33) and normalize it to a unit vector.

Step 7 Repeat Step 3-Step 6 to get the next estimate of uo.

5. A Comparison with Cross Product Estimator

Nehorai and Paldi proposed a cross product estimator [3] by averaging the Poynting vectors. Since the cross product estimator and our weighted inner product estimator have the same objective—estimating the propagation vector of the EM waves, we will draw a comparison of them in the aspects of performance and computational issues.

5.1. Performance comparison

Based on Eq. (16), Eq. (17), Eq. (24), and Eq. (25), using complex imaginary part j to replace e123 and axe1+aye1+aze3 to represent a three-dimensional vector

, the complex vector of the sensor data can

be written as

jEl )(t) — x(t) + jy(t) + %r(t) + jnEi(t) (39)

yHl) (t) — u x (x(t) + jy(t)) + nHr (t) + jnHi (t) (40)

where the superscript "l" means the tth antenna of an EM vector-sensor. The cross product estimator [3] can be equivalently represented in the G3 framework as

= K i Re {)(t ) x(yH )(t ) )*}

where Re{ • } is the operator to extract the real part of a complex vector and (•)* denotes the conjugation of a complex vector.

In Eq. (41), Re {yE°(t)x(yHl)(t))*} can be represented as

Re { )(t) x(t) )} —

u (x(t )-x(t) + y(t )-y(t)) + (x(t) X nHr(t) + y(t) X nHi(t)) +

[nEr (t) X(u X x(t))) + [nEi (t) X (u X y (t)) +

(nEr (t) X nHr (t) + nEi (t) X «Hi (t)) (42)

Using the relationships between the cross product and outer product of two vectors in Eqs. (7)-(8), it is easy to verify that

{yE}(t) x(yH')(i))} = -<Ye(Í)>1 -^(i)>2 -

<emYH(t )>1-<YE(t )> 2 (43)

According to Eqs. (31)-(32), Eq. (43) means the cross product estimator uc satisfies

w x + w

w„ + w,,

Rewrite uo and uc as uouoiei and uc — Z uciei .

i—1 i—1

Let the variance of Uoi and uci be <72oi and ct2 respectively. The following theorem holds.

Theorem 2 When the number of snapshot K is sufficiently large, the variances ct2 and crc2 satisfy

¿o2 (45)

and the equality holds if and only if

Proof See Appendix B.

Theorem 2 implies that if Eq. (46) is satisfied, the cross product estimator has a statistical performance equivalent to our weighted inner product one. If the requirement in Eq. (46) is not met, our weighted inner product estimator is always superior over the cross product counterpart.

In applications, a2 and o2yi are affected by several factors. To begin with, the requirements of the noise statistic in Assumption 2 and Assumption 3 are difficult to satisfy due to the time-varying surrounding environment [25]. Besides, the analysis in the next section will show that, in the presence of the mutual coupling, a2 and a2 are not only dependent on the complex noise, but also rely on the coupling matrix of the antennas. As the coupling matrix is unknown, its influences on a2xi and a2 are also unknown. Finally, since the cross product estimator is the final step of such algorithms as estimate signal parameters via rotation invariance techniques (ESPRIT) [6,8] and self-initiate MUSIC [9], the six-dimensional EM vector for ux and uy is the results of some manipulations of the original sensor data. ct2- and a2 are therefore related to the accuracy of the aforementioned algorithms. Consequently, it is difficult to estimate whether

a2xi and cr2yi satisfy Eq. (46) or not. Seen in this light, our weighted inner product estimator is meaningful, for it always yields the optimal estimation of the wave propagation vector.

As a side note, Eqs. (39)-(40) also indicate that the weighted inner product estimator can be equivalently implemented in the complex domain. The main manipulations are concluded as: 1) separate the array outputs into real and imaginary parts; 2) compute ux and uy from the cross product of the real and imaginary parts respectively; and 3) fuse ux and uy together using Eq. (33).

5.2. Computational issues

According to Eqs. (39)-(40), the computational complexity of ux (or uy) is almost 1/4 of the traditional cross product estimator uc. The reason is that ux (or uy) is equivalent to the cross product of real (or imaginary) part of the sensor outputs, while uc is cross product of the complex sensor outputs before extracting the real part from the result. However, the variance estimation and fusion process of the weighted inner product estimator require extra computational efforts, which are clearly the costs for the performance improvement.

6. Further Discussions

As is well-known, deviations from the true array manifold, typically resulting from mutual coupling effects, can seriously degrade the performance of many high resolution direction finding algorithms. In this section, we will investigate the performances of our estimators in the presence of mutual coupling, and provide a simple approach to compensate the unde-sired effects.

Considering the coupling among the six antennas of the vector-sensor, the received array data in Eqs. (39)-(40) are rewritten as

yEL<ß, t ) =

'y£(t )

yHC(t ). yEl )(t ) yHl)(t )

= Ca(&)s(t ) + nEH(t ) =

AyE(t )

AjH(t ),

where the 6*6 matrix C describes the coupling effects among the antennas, and

4?E(t )

AjH(t ).

= (C -1)a(0)s(t)

symbolizes the model error, which is dependent on the matrix C and the source signal s(t). As array perturbations considered herein are due to unknown mutual

coupling, the matrix C can be modeled without direction of arrival (DOA) dependence [26].

Similar to Eqs. (11)-(15), the array output in the G3 formulation is now

(t) = Ye (t) + AjEr (t) + ej23 Aje (t) (49)

e123YHc (t) = e123YH (t)

+ emAjHr(t) - Ay hi (t) (50)

where AJe = AJEr + jAJEi and AJh = AJHr + JAJhi .

Following Eqs. (49)-(50), the inner product estimator in Eq. (29) can be redefined as

<YEc(t )>1-<e123YHc(t )>2 = Y (t )>• • <et23YH (t )> 2 + Ay& (t) ■ ^ Yh (t )> + AjEr" < e123AyHr (t )> + <Ye (t)\ ■ <e!23AyHr (t )> (51)

Since AyEr (t) and AyHr (t) are dependent on the source signal s(t), it can be checked that the mathematical expectations of last three items in Eq. (51) are not zeros. Following the procedure in the proof of Theorem 1, it is explicit that ux is biased in the presence of mutual coupling. In fact, one can correspondingly verify that all the four estimators, i.e., the traditional cross product estimator (uc), the estima-tor-x (ux), the estimator-y (uy), and the weighted inner product estimator (uo), are biased. As a result, an array calibration process is therefore very necessary. For the special problem investigated in this paper, a simple method is designed as follows.

According to Eq. (47), one can find that

E{ yEHc (&, t)s' (t)} = E{Ca (0)s(t) s * (t)} +

£{«EH(t K (t)} (52)

Based on Assumption 4, one can rewrite Eq. (52) as

E{ yElHc(@, t )s* (t)} = Ca(@K2 (53)

Owing to the finite number of snapshots, Eq. (53) means that

Ca(0) * -iy f [yEHc («, q)s' (?)] (54)

Q&s q=1

where Q is the number of available snapshots.

Given the data collected from n calibration sources of known steering vector {a (@i)}(i=1, 2, —, n) and transmitting signal sequence {s,(t)}(i=1, 2, —, n), we have

CA « D

where A=[a(&l) a(62) — a(&„)] and D=[d(&l) d(&2) — d(&„)] with

d (0,) =-1T ^ [yEHc(®,, q)s*(q)] (56)

Eq. (55) implies that the least square estimate of C is given by

C = AH( AAH)-1 D

Hence, to estimate the mutual coupling matrix C, we need to know the direction and signal sequence of the n calibration sources. For more elaborate algorithms to estimate C without knowing the signal sequence, one can refer to See's method [27]. Using the estimated C, the array output can be calibrated as

yElH (&, t) = CyEHc(0, t) (58)

before implementing the weighted inner product estimator.

7. Simulation Results

In this subsection, simulation examples are presented to verify the performances of the weighted inner product estimator and the analytical predictions in Section 3. A single-source single-vector-sensor scenario with the source parameters &=[22.92° 11.46° 28.65° 23.57°] is considered. In the first simulation, three cases of Theorem 2, i.e.,

1< <1/ <2, + |K|2 / iKik

" " 1» / <2 t \\wv

/ <2 =

Wy|I2 /

wll are taken into

II x 112 > ^ xi ' yi '|| y 112

consideration. Each case is implemented by adjusting the noise and source signal powers. The estimation mean square errors (MSEs) of the ux (estimator-x) and uy (estimator-y) in Eq. (31) and Eq. (32) are respectively set as shown in Fig. 1(a). Namely, the first 150 MSEs of ux in Fig. 1(a) are designed much larger than those of uy and reverse them for the last 200 estimated MSEs, while the MSEs of the two estimators between the 151st and 350th data are almost set the same. To validate the effectiveness of the weighted inner product estimator, a window of length M=50 is chosen. Following the steps in Section 4.2, the MSEs of the weighted inner product and cross product estimators are shown in Fig. 1(b). In the middle of Fig. 1(b), it can be seen that the MSEs of the weighted inner product and cross product estimators are almost the same, which verifies the declaration of Theorem 2 when the requirement in Eq. (46) is met. Moreover, when the requirement in Eq. (46) is not met, the estimate MSEs of our weighted inner product estimator for the first 150 and last 200 data shown in Fig. 1(b) are smaller than those of the cross product one. Such results validate that the estimate errors of the weighted inner product estimator are reduced by the given weighted process, as predicted by Theorem 2.

In the second simulation, the statistical performances of the weighted inner product estimator are tested by Monte Carlo runs. The estimate root mean square (RMS) errors are calculated by RMS(ue) =

ue - u

where ueis the estimate results respect-

tively obtained by the estimator-x, estimator-y, weighted inner product or cross product estimators while u is the true one. The same source parameters as the first

Fig. 1 MSEs of different estimators.

simulation are used, the snapshots for ux and uy are

K=1 000 and the window length is M=50. For each of the nine signal-to-noise ratios (SNRs) ranging from -20-20 dB, 300 estimates of the estimator-x, estimator-y, weighted inner product and cross product estimators are used to obtain their respective average RMS estimate errors. Figure 2(a) depicts their results when the requirement in Eq. (46) is satisfied, i.e., 1 =

<1/ =| K| |2/| Kl |2. It can be seen from Fig. 2(a)

that the RMS errors of the weighted inner product and the cross product estimators are almost the same, which confirms the statement of Theorem 2 again. Slight difference among the RMS errors is due to the fact that the finite data length of K is used. Figure 2(b) illustrates the RMS error analyses of a scenario where t w l| /\\wJ . From it, one can find that

1»<x2i/<y2i

the RMS errors of the cross product estimator are larger than those of the estimator-x and smaller than those of the estimator-y, while the RMS errors of the weighted inner product estimator are the smallest. These results confirm that the weighted inner product estimator is the best.

% V N \ . \ k \ V s Weighted inner product estimator , Cross product estimator - Estimator-* •■ Estimator-j _

-j BSS! ......-.

SNR/dB

(a) Requirement in Eq. (46) is met

-»- Weighted inner product estimator -*■ Cross product estimator "S. -•- Estimator-* • Estimator-v

Fig. 2

SNR/dii

(h) Requirement in Eq, (46) is not met

RMS estimate errors for single-source single-vector-sensor case with no mutual coupling perturbation.

The performance degradation due to the mutual coupling is studied in the third experiment. The configuration of the vector-sensor and the source are identical to that of Fig. 2(b), but the array data are perturbed by a coupling matrix given by

where C1 and C2 are symmetrical Toeplitz matrices with the first rows given by ci=[1 0.2+0.1j 0.2+0.1j] and c2=[0.15+0.15j 0.1+0.05j 0.1+0.05j] respectively. The simulation results with un-calibrated sensor data are displayed in Fig. 3(a), whereas those with calibrated sensor data using Eq. (58) are depicted in Fig. 3(b). Notice that, as expected in Section 6, the mutual coupling perturbation deteriorates the performances of all estimators. The RMS errors for the un-calibrated sensor data are almost a constant with increasing SNRs. By contrast, after the calibration process, the RMS errors of the all the estimators decrease to zero with the increasing SNRs. The above results validate the correctness of the analytical predictions and the compensation approach of Section 6.

y \ Weighted inner product estimator j ■*■ Cross product estimator Estimator-* ■ Estimator-j'

—•—

-20 0 20 SNR/dB (a) Sensor data are not calibrated

i 0.6 5J

.= 0,4

2 0,2 CA

\_ ■ \ V Weighted inner product estimator ♦ Cross product estimator -*■ Estimator-* Estimator-^

-20 0 20 SN R/dl) (b) Sensor data are calibrated

Fig. 3 RMS estimate errors for single-source single-vector-sensor case in the presence of mutual coupling perturbation.

8. Conclusions

The estimation of the propagation direction vector of an EM wave has been investigated in the framework of the G3 by utilizing the measurement data from a single EM vector-sensor. In particular, two estimators for the direction vector are first reported by the inner product between a vector and a bivector in the G3. Since the noise components for the two estimators are uncorrelated, a weighted inner product estimator, in the sense of MMSE, is proposed by fusing the two estimators together. We have analytically proved that the weighted inner product estimator is always statistically optimal. The performances of the weighted inner product estimator are always superior over its cross product counterpart.

References

[1] Nehorai A, Tichavsky P. Cross-product algorithms for source tracking using an EM vector sensor. IEEE Transactions on Signal Processing 1999; 47(10): 2863-2867.

[2] Ko C, Zhang J, Nehorai A. Separation and tracking of multiple broadband sources with one electromagnetic vector sensor. IEEE Transactions on Aerospace and Electronic Systems 2002; 38(7): 1109-1116.

[3] Nehorai A, Paldi E. Vector-sensor array processing for electromagnetic source localization. IEEE Transactions on Signal Processing 1994; 42(2): 376-398.

[4] Wong K T, Zoltowski M D. Uni-vector-sensor ESPRIT for multisource azimuth, elevation, and polarization estimation. IEEE Transactions on Antennas and Propagation 1997; 45(10): 1467-1474.

[5] Wong K T. Blind geolocation/beamforming for multiple wideband-FFH with unknown hop-sequences. IEEE Transactions on Aerospace and Electronic Systems 2001; 37(1): 65-76.

[6] Zoltowski M D, Wong K T. ESPRIT-based 2D direction finding with a sparse array of electromagnetic vector-sensors. IEEE Transactions on Signal Processing 2000; 48(8): 2195-2204.

[7] Li J. Direction and polarization estimation using arrays with small loops and short dipoles. IEEE Transactions on Antennas and Propagation 1993; 41(3): 379-387.

[8] Wong K T, Zoltowski M D. Closed-form direction-finding with arbitrarily spaced electromagnetic

vector-sensors at unknown locations. IEEE Transactions on Antennas and Propagation 2000; 48(5): 671681.

[9] Wong K T, Zoltowski M D. Self-initiating MUSIC-based direction finding and polarization estimation in spatio-polarizational beamspace. IEEE Transactions on Antennas and Propagation 2000; 48(8): 1235-1245.

[10] Zoltowski M D, Wong K T. Closed-form eigenstruc-ture-based direction finding using arbitrary but identical subarrays on a sparse uniform rectangular array grid. IEEE Transactions on Signal Processing 2000; 48(8): 2205-2210.

[11] Monte L, Elnour B, Erricolo D. Distributed 6D vector antennas design for direction of arrival application. IEEE International Conference on Electromagnetic Advanced Applications. 2007; 431-434.

[12] Bull J F. Field probe for measuring vector components of an electromagnetic field. U.S. Patent No. 5300885, 1994.

[13] Au-Yeung C K, Wong K T. CRB: Sinusoid-sources' estimation using collocated dipoles/loops. IEEE Transactions on Aerospace and Electronic Systems 2009; 45(1): 94-109.

[14] Xu Y G, Liu Z W, Wong K T, et al. Virtual-manifold ambiguity in HOS-based direction-finding with electromagnetic vector-sensors. IEEE Transactions on Aerospace and Electronic Systems 2008; 44(10): 12911308.

[15] Wong K T, Yuan X. "Vector cross-product direction-finding" with an electromagnetic vector-sensor of six orthogonally oriented but spatially noncollocating di-poles/loops. IEEE Transactions on Signal Processing 2011; 59(1): 160-171.

[16] Hestenes D. New foundations for classical mechanics. 2nd ed. Boston: Kluwer Academic, 1986.

[17] Dorst L, Mann S. Geometric algebra: a computational framework for geometrical applications Part 1. IEEE Computer Graphics and Applications 2002; 22(3): 24-31.

[18] Mann S, Dorst L. Geometric algebra: a computational framework for geometrical applications Part 2. IEEE Computer Graphics and Applications 2002; 22(4): 58-67.

[19] Sabbata V D, Datta B K. Geometric algebra and applications to physics. New York: CRC Press/Taylor & Francis Group, 2007.

[20] Ell T A. Multi-vector color-image filters. IEEE International Conference on Image Processing. 2007; 245-248.

[21] Menti A, Zacharias T, Milias-Grgitis J. Geometric algebra: a powerful tool for representing power under nonsinusoidal conditions. IEEE Transactions on Circuits and Systems I: Regular Papers 2007; 54(3): 601-609.

[22] Ebling J, Schenermann G. Clifford Fourier transform on vector fields. IEEE Transactions on Visualization and Computer Graphics 2005; 11(4): 469-478.

[23] Deschamps A. Part II-geometrical representation of the polarization of a plane electromagnetic wave. IRE Proceeding. 1951; 540-544.

[24] Xu L J, Zhang J Q, Yan Y. A wavelet-based multisensor data fusion algorithm. IEEE Transactions on Instrument and Measurement 2004; 53(6): 1539-1545.

[25] Swindlehurst A L, Kailath T. A performance analysis of subspace-based methods in the presence of model er-

rors, Part I: the MUSIC algorithm. IEEE Transactions on Signal Processing 1992; 40(7): 1758-1774.

[26] See C M S. Method for array calibration in high-resolution sensor array processing. IEEE Proceeding of Radar, Sonar, and Navigation 1995; 142(3): 90-96.

[27] Rao C R. Linear statistical inference and its application. 2nd ed. Hoboken: John Wiley & Sons, 1972.

Biographies:

JIANG Jingfei received B.S. and M.S. degrees in the Electronic Engineering Department from Fudan University, in 2008 and 2011 respectively. His research interests include array signal processing and geometric algebra with applications in signal processing and image processing. E-mail: 082021030@fudan.edu.cn

ZHANG Jianqiu received B.S. degree from the Electronic Engineering Department from East of China Institute of Engineering, Nanjing, China, in 1982, M.S. and Ph.D. degrees from the Department of Electrical Engineering, Harbin Institute of Technology (HIT), Harbin, China, in 1992 and 1996, respectively. He is currently a professor with the Department of Electronic Engineering, Fudan University, Shanghai, China. From 1999 to 2002, he was a Senior Research Fellow at the School of Engineering, University of Greenwich, Medway Campus, Chatham Maritime, U.K. In 1998, he was a visiting research scientist at the Institute of Intelligent Power Electronics, Helsinki University of Technology, Espoo, Finland. He was an associate professor from 1995 to 1997 and a lecturer from 1989 to 1994 with the Department of Electrical Engineering, HIT. From 1982 to 1987, he was an assistant electronic engineer at 544th factory, Hunan, China. His main research interests are signal processing and its application for advanced sensors, intelligent instrumentation systems and control, and communications. E-mail: jqzhang@fudan.ac.cn

Appendix A: Proof for Theorem 1

Let z(t) = 5<Ye(/)>1 • <e123YH(t)> 2, we have E{z (t)} = E ( aE )u + em E (BE ) +

emE(CE) + E (d e ) (A1)

where aE (t) = - x(t) • x(t), BE (t) = x(t) a «Hr (t), CE (t) = «Er(t)a«Hr(t), and dE(t) = «Er(t) • (u ax(t)).Let x(t) =

xjej + x2e2 + x3e3 and «Hr(t) = nHr1ej + nHr2e2 + nHr3e3, then

BE = x(t) A «Hr (t) = (x1nHr2 - x2nHr1)e12 +

Hr1 x1nHr3 )e31 + (x2n Hr3 x3nHr2 )eB (A2)

Since the sources and the noise are uncorrelated (see Assumption 4) and x(t) is a linear combination of the real and imaginary parts of the source signal, we have E(BE(t))=0. Similarly, under Assumption 4, we also have E(Ce(/))=0. Using the rule that a • (b a c) = (a • b) c -(a • c) b and the well-known (a • c) b- (a • b) c = ax(b x c), dE(t) can be rewritten as

dE(t )=-«Er (t) x (u x x(t)) (A3)

Similar to x(t), uxx(t) is also a linear combination of the real and imaginary parts of the source signal. Like Eq. (A2), we also have E(dE(t))=0. With the above results, it is clear that

E( z(t )) = E(aE)u = -[ pi2 E( sy(t)) + p2r E (s2(t)) + p2i E(sl(t))]u

Let < =E(sr2(t)) and a2 =E(s 2(t)), then pi = (p2 + P2i) <sr + p22r < is a finite constant. As a result, the expectation of z(t) is a constant vector with the same direction as u. Since z(t) is an independent identical distribution (IID) random vector with finite constant expectation, by the Kolmogorov strong law of large

numbers [27], we have

be o^ri . Since wx =—^ z(t) and z(t) is a IID ranK t=1

dom vector (see Appendix A) with a constant expectation -p^u and a constant variance matrix generated by the random vector ei23BE+ei23CE+dE, by the Lind-berg-Levy central limited theorem [25], we have

■N (-p\u,CTWa

where N(a, a2) denotes a Gauss random process with a mean equivalent to a and a variance given by a2. As a result, we have

N ( ,ct2 )

=K 5z( >■

-> = - px u

where cti=ctL

Since u is a unit vector, we have

With similar process, one can also obtain that

where p2 = ( pi2 + p2r ) CTs2 + p22iCTs2 .

Appendix B: Proof for Theorem 2

As proved in Appendix A, when K is sufficiently large we have

E (wx ) = - p2 u, E ( wy ) = - py,u

\\wa\2 = px

llwy| I2 = py

Similar to the facts shown in Appendix A, we can get

E (wc ) = -( p.2 + p^u

Similarly, we also have

uyi = -T

N(u, ,ct2)

'N (, ctC )

where avl = ctwvi. / p„ and CTci =( CTwx, + ctV2V, )/ (px +

py2 )2. As shown in Eq. (36), the optimal variance of the weighted inner product estimator is given by

ct2 =-

wxi wyi

----I--—

4 2 4 2

p,CTwvr + pvCTwx

2 2 CTci - CToi =

2 2 2 2 2 (pyCTwxi - p,CTwy, )

(p,2 + p2)2( pictii + pilCTwW xi)

>0 (B10)

in which the equality holds if and only if P2 < , = 0, i.e.,

22 px CTwy

= p,2 + py

Rewrite wx as wx = ^ wxiet and let the variance of wx