Candecomp/Parafac based Array Processing David Brie CRAN UMR 7039 - - PowerPoint PPT Presentation

Candecomp/Parafac based Array Processing

David Brie CRAN UMR 7039 - Universit´ e de Lorraine - CNRS Winter school : Search for Latent Variables : ICA, Tensors and NMF 2 – 4 February 2015, Villard de Lans

1

Candecomp/Parafac Decompositions CP based Array Processing CP-Based Vector Sensor Array Processing

2

Candecomp/Parafac Decompositions

3

Useful matrix operations

A = [a1 · · · aR], (I × R) B = [b1 · · · bR], (J × R) – Kronnecker product A ⊗ B =    A(1, 1)B A(1, 2)B · · · A(2, 1)B A(2, 2)B · · · . . .    , (IJ × R2) – Khatri-Rao product A ⊙ B = [a1 ⊗ b1 · · · aR ⊗ bR], (IJ × R) – Outer product a ◦ b, (I × J), a ◦ b(i, j) = a(i)b(j) i.e. a ◦ b = abT , rank 1 matrix. a ◦ b ◦ c, (I × J × K), a ◦ b ◦ c(i, j, k) = a(i)b(j)c(k) i.e. rank 1 tensor.

4

CP Decomposition of an order 3 tensor

CP : Candecomp/Parafac [Harshman, 1970-1972], [Carroll - Chang, 1970] X ≈

R

• r=1

ar ◦ br ◦ cr = [ [A, B, C] ] A(I × R), B(J × R), C(K × R)

= σ1 a1 b1 σR aR bR + + · · · = A B cR c1 C X

Tensor rank : minimal number of rank 1 tensors to represent X . Alternative writing Slice : Xk = ADk(C)BT , k = 1, · · · , K Unfolding : X(KJ×I) = (B ⊙ C)AT

5

CP based Array Processing

Antenna arrays

6

Direction Of Arrival (DOA)

Far field EM wave ⇒ plane wave Direction cosine k k =   sin(θ) cos(φ) sin(θ) sin(φ) cos(θ)  

b

x y z S k θ φ m

× × × × ×

Narrowband assumption Received signal on sensor 1 : s(t)ej 2π

λ t

Received signal on sensor m : s(t − tm)ej 2π

λ (t−tm) ≈ s(t)ej 2π λ tm

• complex envelope

ej 2π

λ t

carrier

Phase factor : ej 2π

λ tm = ej 2π λ kT dm where dm is the position of the mth sensor 7

Multiple Sources

P narrowband sources Steering vector of the pth source : a(kp) = [ej 2π

λ kT d1 · · · ej 2π λ kT dM ]T

X = AST + N dim(X) = (M × K) A =

• a(k1), . . . , a(kP)
• (M × P)

S =

• s1, s2, . . . , sP
• (K × P)

DOA estimation ⇒ estimation of A ⇒ estimation of kp, p = 1, · · · , P Unambiguous estimation ⇒ inter-element spacing ≤ λ/2 Eigen decomposition approaches (Music, Esprit)

8

Rotational Invariance

Basic idea of Esprit : the array can be decomposed into 2 translated but

• therwise identical subarrays ⇒ GEVD of stacked data corresponding to each

subarray

bc bc bc bc bc bc bc bc bc

Two non-overlapping subarrays

bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc

Two overlapping subarrays

bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc

The idea of Esprit is difficult to extend to multiple invariances !

9

CP-based Array Processing

Seminal work of Sidiropoulos, Bro and Giannakis (2000) ⇒ Handling multiple invariances ⇒ Joint estimation of both DOA and sources

bc bc bc bc bc bc bc bc bc

An array with multiple invariances

bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc

X =

• A1 ⊙ A2
• ST + N

10

Multi-scale Array

bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc Level 1 Level 2 Level 3

dl1,l2,...,lN =

N

• n=1

dn

ln

11

Data Model (1)

Phase factor for one narrowband source : al1,l2,...,lN (k) = exp

• j 2π

λ

N

• n=1

kT dn

ln

• =

N

• n=1

exp

• j 2π

λ kT dn

ln

• .

Array manifold for the entire sensor array a(k) = a1(k) ⊗ · · · ⊗ aN(k), with an(k) =     ej(2π/λ)kT dn

1

. . . ej(2π/λ)kT dn

Ln

    ⊗ Kronecker product

12

Data Model (2)

P narrowband sources A1 =

• a1(k1), . . . , a1(kP )
• (L1 × P)

. . . AN =

• aN(k1), . . . , aN(kP )
• (LN × P)

S =

• s1, s2, . . . , sP
• (K × P)

X =

• A1 ⊙ · · · ⊙ AN
• ST + N,

⇒ CP Model of order N + 1

Single snapshot ⇒ CP Model of order N

13

Model identifiability

– [Sidiropoulos and Bro 2000] : sufficient condition for the uniqueness of CP decomposition for N-way arrays – P sources with distinct DOAs and not fully correlated – Number of snapshots greater than number of sources (K > P)

N

• n=1

min(Ln, P) ≥ P + N. – Case of a single snapshot

N

• n=1

min(Ln, P) ≥ 2P + N − 1.

14

Parameter estimation (1)

Estimation of the steering vectors by CP decomposition of the data ˆ ap

n : nth level estimated steering vector for the pth source

Estimating the DOA parameters for the pth source ⇒ minimization of the following criterion : IN(kp) =

N

• n=1

Jn(kp). Jn(kp) = ˆ ap

n − an(kp)2, n = 1, . . . , N,

15

Parameter estimation (2)

IN(kp) : non-convex criterion highly non-linear ⇒ Sequential strategy similar to a GNC approach

b b b

I1(kp), ∆ ≤ λ/2 I2(kp) I3(kp)

⇒ The ordering of the sub-arrays is important !

16

Parameter estimation (3)

First Stage : Estimate A1, . . . , AN by CP decomposition of the data. Second Stage :

◮ For p = 1, . . . , P and

for n = 1, . . . , N compute k∗

p,n = argmin kp

In(kp).

◮ Output : The estimated parameters for the P sources :

ˆ kp = (ˆ up, ˆ vp, ˆ wp) = k∗

p,N with p = 1, . . . , P .

17

Simulation results (1)

bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc bc

Niveau 3 X Y O 10 λ λ 2

– Two scales L1 = 5, L2 = 4 – Three scales L1 = 5, L2 = 2, L3 = 2 ⇒ single snapshot case – Two sources with distinct DOAs – Comparison with ESPRIT [Wong, Zoltowski 1998]

18

Simulation results (2)

SNR dB 5 10 15 20 25 30 35 40 R.M.S.E for (u, v) estimation (in radians) 10-5 10-4 10-3 10-2 10-1 100

CP-based method for Source #1 CP-based method for Source #2 MUSIC-ESPRIT for Source #1 MUSIC-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

SNR dB 5 10 15 20 25 30 35 40 R.M.S.E for (u, v) estimation (in radians) 10-5 10-4 10-3 10-2

CP-based method for Source #1 CP-based method for Source #2 MUSIC-ESPRIT for Source #1 MUSIC-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

CRMSE , K= 5 CRMSE , K=20

# of snapshots 100 101 102 R.M.S.E for (u, v) estimation (in radians) 10-4 10-3 10-2 10-1

CP-based method for Source #1 CP-based method for Source #2 MUSIC-ESPRIT for Source #1 MUSIC-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

SNR dB 10 20 30 40 50 60 70 R.M.S.E for (u, v) estimation (in radians) 10-6 10-5 10-4 10-3 10-2 10-1 100

CP-based method for Source #1 CP-based method for Source #2 MUSIC-ESPRIT for Source #1 MUSIC-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

CRMSE , SNR = 15 dB CRMSE , K=1

19

Simulation results (3)

– Four scales L1 = L2 = L3 = L4 = 3 – Two sources with distinct DOAs – Comparison with ESPRIT [Wong, Zoltowski 1998] and Tensor-ESPRIT [Haardt et al. 2008]

20

Simulation results (4)

SNR dB 5 10 15 20 25 R.M.S.E for (u, v) estimation (in radians) 10-4 10-3 10-2

CP-based method for Source #1 CP-based method for Source #2 Standard Tensor-ESPRIT for Source #1 Standard Tensor-ESPRIT for Source #2 Unitary Tensor-ESPRIT for Source #1 Unitary Tensor-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

SNR dB 5 10 15 20 25 R.M.S.E for (u, v) estimation (in radians) 10-4 10-3

CP-based method for Source #1 CP-based method for Source #2 Standard Tensor-ESPRIT for Source #1 Standard Tensor-ESPRIT for Source #2 Unitary Tensor-ESPRIT for Source #1 Unitary Tensor-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

CRMSE , K= 10 CRMSE , K=50

# of snapshots 100 101 102 R.M.S.E for (u, v) estimation (in radians) 10-4 10-3 10-2 10-1

CP-based method for Source #1 CP-based method for Source #2 Standard Tensor-ESPRIT for Source #1 Standard Tensor-ESPRIT for Source #2 Unitary Tensor-ESPRIT for Source #1 Unitary Tensor-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2 Inter-grid spacing, in units of half-wavelengths 5 10 15 R.M.S.E for (u, v) estimation (in radians) 10-4 10-3 10-2 10-1 CP-based method for Source #1 CP-based method for Source #2 Standard Tensor-ESPRIT for Source #1 Standard Tensor-ESPRIT for Source #2 Unitary Tensor-ESPRIT for Source #1 Unitary Tensor-ESPRIT for Source #2 Dual-Size ESPRIT for Source #1 Dual-Size ESPRIT for Source #2 CRB for Source #1 CRB for Source #2

CRMSE , SNR = 6 dB CRMSE , K=5, SNR = 6 dB

21

CP-Based Vector Sensor Array Processing

22

Polarized EM wave

α β

α : orientation β : ellipciticity

23

Response of the vector sensor to a polarized wave

x y z ex ey ez hx hy hz rm u(φ, ψ) φ ψ

gk g(αk, βk) = gφ(αk, βk) gψ(αk, βk)

• =
• cos αk

sin αk − sin αk cos αk cos βk j sin βk

• bk

e(φk, ψk, αk, βk) h(φk, ψk, αk, βk)

• =

        − sin φk − cos φk sin ψk cos φk − sin φk sin ψk cos ψk − cos φk sin ψk sin φk − sin φk sin ψk − cos φk cos ψk        

• F(φk,ψk)

gk.

24

CP model of a Vector Sensor Array

✻ ✲ ✠x

y z sk φk ψk

✰

∆x

✠ ✒

X = (A ⊙ B)ST + N A = [a1(φ1, ψ1), . . . , aK(φK, ψK)] (M × K) steering matrix of the K sources B = [b1(φ1, ψ1, α1, β1), . . . , bK(φK, ψK, αK, βK)] (6 × K) polarization matrix of the K sources

25

Does polarization really matter ?

Two key factors directly affect the performance of polarized source estimation : – the polarization separation – the angular separation

−80 −60 −40 −20 20 40 60 80 10

−2

10

−1

10 Orientation Angle of the Variable Source RMSE ξ=0.03° ξ=0.3° ξ=0.6° ξ= 1° ξ= 2° ξ= 4° ξ=8° ξ=20°

(a) β1 = β2 = 0◦, α1 = 0◦ varying α2

−40 −30 −20 −10 10 20 30 40 10

−2

10

−1

10 Ellipticity Angle of the Variable Source RMSE ξ=0.03° ξ=0.3° ξ=0.6° ξ= 1° ξ= 2° ξ= 4° ξ=8° ξ=20°

(b) α1 = α2 = 0◦, β1 = −45◦ varying β2

Figure: The RMSE of source DOA estimation versus their polarization separation

26