DIRECTIONOFARRIVAL
ESTIMATION
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
Introduction
Finding the direction of arrival of a wavefront is identical to the problem of
finding power concentration in the frequency domain. In consequence.
DOA estimation is the root problem of spectral density estimation dBm/Hz
or just dBm/degree of solid angle.
Differences: 3D search, non uniform sampling, wave propagation effects,
no rational models, robutsness to missmatched is the must of DOA
estimation methods.
The problem: Giving a set of N snapshots from an aperture of Q sensors,
to estimate the number of sources and their angles of arrival. When N>
10.Q the covariance is almost stabiliced.
We will start with the narrowband problem for uncoherent point sources
in the far field. The effect of coherent sources will be mentioned for each
procedures and, finally, the wideband case will be presented.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
Thenarrowbandsnapshotand
covariance
NS
X n ai (n).S s w n S s .a n w n
i 1
NS
R S s .P s .S s R o P (r ).S r .S R o
H
r 1
Our goal: To estimate the
steering vector
contained in the matrix
of DOAs
H
r
Waveform and power
of the sopurces is not
our priority concern
Noise motivates the
problem and converts
a deterministic well
defined problem in an
estimation problem
To families of procedures: SCANNING METHODS
NULLING PROCEDURES (SUPERRESOLUTION)
Array Processing Miguel Angel
September 13
3
ML BASED
METHODS
(OPTIMUM)
Lagunas DOA
Estimation
SCANNINGFramework
Beamformer Response (to be designed A)
Scanning
90
120
1
0.8
60
Steering Vector S
(Direcyion to scan)
0.6
Polar response
150
30
0.4
0.2
180
210
330
240
300
In order to steer the beam to the desired
direction a constraint is mandatory in the
design
H
A .S 1
270
Elevation in degrees
Just in case others directions should be nulled or
attenuated (those where known stations or RF
equipment is located. In Array
this Processing
case te constrains
Miguel Angel
September 13
would be NC<Q
Lagunas DOA Estimation
A H .C H f
Ince a beamformer A have been selected, we measure the power at the output
of the beamformer as:
P ( S ) A .R. A mW .
H
The estimate is formed by the quotient of
the power measured divided by the
beamwidth of the beamformer.
TWO alternatives to derive a close expression of the effective beamwidth:
- From the equation that relates the output power with the beamformer
directivity and the incoming power density.
- From a calllibration problem.
P( S ) A .R. A
Field of View
Beamformer
directivty, steered on
the desired direction
September 13
A , d s d
2
Power density, i.e. the value to be
Miguel
Angel
estimated
by any DOA procedure
Array Processing
Lagunas DOA Estimation
Assuming that the beam is narrow with beamwidth Bw around the steered
angle
P( S ) A .R. A
Field of View
Estimate of s(), due to the
approximation, of the power density
A , d d
2
Field of View
A , d d A A
2
In summary, given the beamformer A the power level estimate and the power
density (local maxima will provide DOA estimates will be:
P ( S ) A .R. A
H
(S )
September 13
A .R. A
H
A A
Array Processing Miguel Angel
Lagunas DOA Estimation
The calliobration alternative consists on choosing the beamwidth Bw in such a
way that with undirectional noise only present in the scenario , i.e. Covariance
equal to 2.I, the estimat e will be equal to 2
H
(S )
A .R. A
Bw
R I
Bw A . A
H
Clearly the callibration requires that
This beamwidth can be defined as the bandwidth of an ideal
brick shape beamformer such that both, the original and
the ideal, produce the same output power when only nondirectional noise is present in the scenario.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
SCANNINGprocedures differ on
the way the scanning beamfomer
is designed
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
ThePhasedArray(PA)Method
The easiest procedure for beamforming is the so-called phased array
method. Imagine that we desire to steer the beam toward the broadside,
in this case all the entries of the beam are one resulting in an almost nooperation processing. In fact, the mathematical ground for this design is
just to select the beamfomer which produces the maximum dot product
with the steering vector of the direction to scan.
In summary, the
beamformer is selected from
a non-bias constraint in the
scanning direction and
minimum norm or response
to the leakage produced by
the non-directional noise
September 13
A .S 1
H
H
A .A
min
The solution is an
only-phase
beamformer known as
Array Processing Miguel Angel
the DOA
phased
array.
Lagunas
Estimation
1
A H .S
S .S Q9
Using this beamformer on the power and
the spectral density estimators we have:
P ( , ) P ( S ) A .R. A
H
1
Q2
.S H .R.S
A .R. A
1 H
( , ) ( S )
.S .R.S
H
Q
A .A
Note that the power level and the density level
differ only in a constant. The reason for that is that
this procedure is a constant-bandwidth scan since
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
1
A .A
Q
H
10
PA estimates for an ULA array of 15 elements and using
3000 snapshots to compute the array covariance matrix R
Resolution is limited by the
aperture size (beamwidth of the
scanning beam.
The major drawback of PA is
the large leakage suffered on
DOAs close where a source is
present
Estimacion de DOA: Phased array
15
10
dB
-5
-10
-100
September 13
-80
-60
-40
-20
0
20
40
Elevacion en
grados
Array Processing
Miguel
Angel
Lagunas DOA Estimation
60
80
100
11
For a planar aperture the performance is
even worse since in somes directions the
aperture is working down to 5 effective
sensors instead of 15.
Metodo phased array
1
0.8
0.6
Slowness/azimuth
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1
September 13
-0.5
Array Processing Miguel Angel
Lagunas DOA Estimation
0
Sur
0.5
12
Datadependentbeamforming
Since we desire to measure power impinging from the desired direction,
power entering to the array from other directions will introduce a positive bias
(leakage) on the power measured. In order to minimize leakage, the
objective should be:
P ( S ) A .R. A
H
MIN
In addition the constraint of 0
dB gain on the steering
direction have to be added.
September 13
Knowing that the bias will be positive we
have to minimize the output power.
A .S 1
H
Note that additional constrains can be added, i.e. null to
engine noise source in towed arrays, known RF stations
around, clutter, Array
etc. Processing Miguel Angel
Lagunas DOA Estimation
13
TheMLMbeamformerfromCapon
1
The resulting beamformer is:
R .S
H
S .R .S
Note that this beamformer suffers from desired degradation when a
coherent source is present on the scenario. The reason for coherent
sources degrading the performance of MLM is that they promote negative
bias on the power estimation, i.e. the beam steers the coherent source to
minimize the output power when steered to the desired.
On mathematical terms, coherent sources degrade the
rank of the covariance matrix of the sources degrading
severely the performance of the method.
Regardless the method was denominated by Capon as maximum likelihood,
this is not the case in general. ONLY when the beam is steered to one, and
single, existing source in a non directional noise the power measured is MLM.
Nevertheless, the optimum beamformer in this case coincides with a phased
Array Processing Miguel Angel
September
13
14
array
Lagunas DOA Estimation
In order to do not produce any miss-understanding related with the original name
of MLM, authors use to refer to this beamformer as minimum variance
beamformer in the sense that it minimizes the variance at the beamforer output.
After using this beamformer for the
power level and power density
estimates we obtain:
S S
P S
H
A .A
P S
1
H
S .R .S
A .R. A
H
watts
A A
S .R .S
S .R .S
Scanning, as in all the procedures described herein is done by changing
the steering angles within the steering vector S
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
15
Power level estimate MLM for an ULA array
Better resolution than PA
Almost no leakage from
existing sources
DOA Estimation: MLM Capon
15
10
Accurate power
level estimation
dB
0
-5
-10
-15
-100
September 13
-80
-60
-40
-20
0
20
40
Elevacion en grados
Array Processing Miguel Angel
Lagunas DOA Estimation
60
80
100
16
The improvement with respect PA is very evident for a planar aperture
Metodo MLM
1
0.8
0.6
Slowness/azimuth
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1
September 13
-0.5
0
Sur
Array Processing Miguel Angel
Lagunas DOA Estimation
0.5
1
17
For the power density estimate the resolution further improves The leakage
level is also improved. Of course power level is lost on density plots.
Metodo Normalizado-MLM
1
0.8
0.6
Slowness/azimuth
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1
September 13
-0.5
0
0.5
Sur Miguel Angel
Array Processing
Lagunas DOA Estimation
1
18
In order to show the performance of NMLM versus traditional MLM we need
just to go to scenarios with closely spaced sources
DOA estimacion: NML Lag
15
DOA Estimation: MLM Capon
15
10
10
dB
dB
5
0
-5
-10
-15
-100
-80
September 13
-60
-40
-20
0 -8 20
Elevacion en grados
-640
60 -4 80
100-2
0
Elevacion
en grados
Miguel
Angel
Array Processing
Lagunas DOA Estimation
6
19
Some authors, in order to use
traditional MLM for density proposed
(wrong) to use a constant bandwith
1/Q to produce a density estimate
from the power estimate.
Q
H
1
S .R .S
watts/degree
The formal proof of the superiority of NMLM versus MLM and PA, in terms of
resolution we just need to use the fact that for density estimates we have:
u S .R
vR
1/ 2
1 / 2
u u . v v u
H
S .R.S
Q
S .R.S S .R .S Q H 1
Q
S .R .S
H
and
uS
vR S
S .S S
H
September 13
u u . v v u
H
.R .S S .R .S
H
S .R .S
H
S .R .S
Thus, assuming all
get the same level at
the source location,
NMLM falls down
faster than the other
two methods.
Q
H
1
S .R .S
Array Processing Miguel Angel
Lagunas DOA Estimation
20
Note the degradation of NMLM (and MLM) for coherent sources. It is
worthwhile to remark that this effect is identical to the desired cancellation on
SLC and GSLC beamformers.
DOA estimation: Phased (blue actual location)
15
10
DOA estimation: NML Lag
COHERENT
SOURCES
15
10
DOA Estimation: MLM Capon
UNCOHERENT
SOURCES
10
dB
5
5
dB
dB
-5
0
-10
-15
-5
-15
-10
-5
September 13
0
5
10
Elevation in degrees
-5
-10
-5
15
0
5
10
Elevation in degrees
20
25
-10
-15
15
-10
Array Processing Miguel Angel
Lagunas DOA Estimation
20
-5
25
0
5
10
Elevacion en grados
15
20
25
21
Nullingprocedures(Super
resolutionmethods)
Remarks:
- The size of the aperture bound the minimum beam-width of a beamformer
response.
- Size do not preclude close located zeros of the beamforer response.
- The number of elements bounds the leakage suffered by a beamformer
- The number of elements (minus one) limits the number of zeros of the
beamformer response.
IDEA:
Instead of looking for maximum power associate to the presence of
sources, change the procedure such that minima (or zeros) are
associated to the presence sources SUPER-RESOLUTION
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
22
HOW TO DO IT:
Design a beamformer A such that minimices the output power. To do so the
objective is the same that in scanning methods.
Set some constraint to prevent the trivial solution, i.e. the zero beamformer.
Compute the response of the beamformer, since to do properly its job, it
will present minima or zeros response to those DOA where a source is
present in the scenario.
A .R. A
GENERAL FORMULATION:
Objective
f A c
Constraint
Source detector as the maxima
of the inverse of
Array Processing Miguel Angel
September 13
Lagunas DOA Estimation
the beamfomer spatial response.
MIN
s S
1
H
S A
2
23
Differenceofbeamformersfor
scanningorsuperresolution
MLM
Optimum as rx-1*sd
-10
PA
Quiescent
-15
10
SUPERRESOLUTION
eigenvector
# 1 eigenvalue...1.8376
-20
Array factor dB.
0
-5
-25
-30
-5
Array factor dB.
Array factor dB.
-10
-10
-35
-15
-80
-15
-80
-60
-40
-20
-20
0
Elevation
20
40
60
-60
-40
-20
0
Elevation
20
20
40
60
80
40
60
80
80
-25
-80
September 13
-60
-40
-20
Array ProcessingElevation
Miguel Angel
Lagunas DOA Estimation
24
Thespatiallinearpredictor
The constrain will be setting to one the
beamformer weight of one element of the
aperture.
A .R. A
MIN
A 1 1 with 1 0 .. 0 1 0 .. 0
H
The solution for the beamformer is:
R 1 .1
1H .R 1 .1
And, the corresponding estimate is:
Conformador
de
anulacion o
predictor
lineal
Error
prediccion
de
(S )
H
1 .R .S
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
25
Pisarenko
The constrain will be setting to one the
norm of the beamformer which is the
same that asking for minimum response
to the un-directional noise
H
A .R. A
MIN
A A 1
The solution for the beamformer is the
Minimum eigenvector of this problem:
A min eigenvector of
R.e .e
And, the corresponding estimate is:
( S )
e
September 13
H
min
.S
Array Processing Miguel Angel
Lagunas DOA Estimation
26
COMMENTS:
Theestimatesproducefalsepeaks.
Thenumberofelementsisgreaterthanthe
numberofexistingsources.
ForQelementstheresponsemayhaveupto
Q1zeros,beingNSthenumberofsources,it
willbeQ1NSfalsesources.
Pisarenko andProny solvedtheproblemof
labelingactualsources.
S,R,Ksolvedinanelegantmannerthis
problemofresolutionmethodswiththeso
calledMUSIC
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
27
MUSIC(Goniometre)
Super-resolution is valid for point sources and, in general for un-directional
front-end noise. For this reason, people use to refer these procedures as
point source detectors (far-field scenarios).
Under the above hypothesis, the
array covariance matrix has a
solid structure that can be used to
improve the point source location
performance.
NS
R Ps .S s .S sH 2 .I S .P.S H 2 .I
s 1
SIGNAL SUB-SPACE
- Dimenssion NS, or rank
NS, is formed from NS
independent steering
vectors
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
NOISE SUB-SPACE
- Full rank Q and formed
by Q vectors orthogonal
and equal to the
columns of the identity
matrix.
28
- The dimension of the array covariance is Q, thus the covariance matrix can
be described from an orthonormal base of dimension Q.
- Within this space of dimension Q there is a subspace of lower dimension NS
where the steering vectors we are interested for, form a non orthogonal base of it.
- Let us assume that vector eq (q=1,Q) is one of the Q vectors describing the full
space of the aperture, the power associated with this vector will be:
e R.e q q
H
q
- It is obvious that those vectors that take part of the signal subspace will have
more power associated than those which only, and just only, describe the noisesubspace.
- We assume that vectors form an orthogonal base, the noise vectors will be
orthogonal to the vectors describing the signal subspace and in consequence THEY
WILL BE ORTHOGONAL TO THE
VECTORS
OF THE SOURCES
ArraySTEERING
Processing Miguel
Angel
September 13
Lagunas DOA Estimation
29
The eigenvectors of the covariance matrix is the proper base to choose for our
analysis.
.e .e
q
H
q
E.D.E H
q 1
being
E e1 . e NS
e NS 1 . e Q
the eigenvectors satisfy that
R.e q q .e q
Note that the eigenvalue
is just the power
associated with the
corresponding
eigenvector when we
look at it as a
beamformer
As the largest eigenvalues are associated with the point sources we may say
that the corresponding NS eigenvectors fully describe the signal sub-space.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
30
R s e s , q s , q e s , q
q 1,..., N s
or
Ns
Ns
R s q S q S s , q e s , q e s , q
q 1
H
q
Note that the eigenvectors of the signal
subspace are also the eigenvector of the
array covariance and the eigenvalue is
increased in the noise level
q 1
The signal subspace is described
by the two largest
eigenvectors which
are orthogonal
Signal sub-space
(Since Ns=2 it is a
plane)
The steering vectors of
the sources, not
orthogonal, have to be
31
in the plane.
THE NOISE EIGENVECTORS, ALL OF THEM WITH
EIGENALUE EQUAL TO 2, WILL BE
ORTHOGONAL TO THE SOURCE STEERING
Array Processing Miguel Angel
September 13
VECTORS
Lagunas DOA Estimation
MUSIC
- Compute eigenvectors and eigenvalues
- Decide the dimension of the signal sub-space
- Form and find the maxima of the MUSIC estimate
(S )
Q
q NS 1
S .e q
Note that every beamformer or eigenvector have Q-1 minima but only at
the NS source location they will coincide. In this manner MUSIC
removes the problem of false peaks. Note also that Pisarenko reduces to
MUSIC when the dimension
of the noise subspace is one, i.e. is identical
Array Processing Miguel Angel
September
13
32
when
number of actual sources
is Q-1.
Lagunas DOA
Estimation
The major problem of MUSIC, general for any subspace-based method, is
to decide the length of the signal/noise subspace.
- When Sn (dimenssion)<NS (actual # of sources) the estimate tends to
show less peaks than actual sources.
- When Sn=NS the result is excellent even for very
small source separation
Music.-Sinusoides en ruido
35
30
- When Sn>NS the estimate
shows false peaks.
25
S(w) in dB
20
15
10
5
0
-5
-10
-15
September 13
200
Array Processing Miguel Angel
Lagunas DOA Estimation
400
600
frequency
800
1000
1200
33
SELECTING THE DIMENSSION OF THE SUBSPACES
-
A single source with low signal to noise ratio
- When sources are closely located the first eigenvalue increases and
the second decreases
Example: NS=2 (Two actual sources)
Note that eigenvectors are beams (Eigenbeams).
The maximum takes maximum energy from all sources present.
The minimum does the same but with the constraint of
being orthogonal to the previous one.
Eigenvalue 1
Max eigenbeam
Min. Eigenbeam
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
34
Difficult to decide the signal subspace dimenssion !!!
Ordered eigenvalues for MUSIC
10
Eigenvalue
September 13
6
7
Order
Array Processing Miguel Angel
Lagunas DOA Estimation
10
11
12
35
MusicPerformance
Autovalores
300
250
Metodo de estimacion de DOA: Music
50
45
200
40
35
150
dB
30
25
20
100
15
10
50
0
5
0
-8
September 13
-6
-4
10
Numero
Array
Processing Miguel Angel
Lagunas DOA Estimation
-2
0
Elevacion en grados
15
36
Metodo Music
1
0.8
0.6
Slowness/azimuth
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1
September 13
-0.5
0
Sur
Array Processing Miguel Angel
Lagunas DOA Estimation
0.5
37
DOAEnhancedMethodsfrom
noisesubspaceversions.
Pisarenko predicted MUSIC by further generalizing MUSIC (Almost 50 years
before MUSIC was reported) by the so-called potential density estimators.
First, note that any positive or negative
power of a covariance matrix is reflected in
its eigenvalues
R mq .e q .e qH
m
q 1
As a consequence taking the Capons power density estimate,
MLM
(S )
1
H
MLM
(S )
Can be written as: P
1
Q
1
q . S .e q
S .R .S
q 1
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
38
Knowing that noise eigenvalues are orthogonal to the steering vectors of the
sources, artificial resolution can be added by just using only the noise
subspace. This is refereed as the noise subspace version of Capons method
and was reported by Jhonson.
JH ( S )
1
Q
q NS 1
q1
S .e q
The noise subspace version of NMLM, also improving resolution is:
Q
NMLM ( S )
2
H
1
q . S .e q
q 1
Q
2
H
2
q . S .e q
q 1
In fact, Pisarenko proved the following statement:
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
39
For any function f(x) such that the inverse exist and is continuous g(f(x))=x the
family of estimates
lim Q ESTIMATE ( S ) g S f R S ACTUAL ( S )
H
Converges to the actual density when the size of the aperture goes to infinity
yet preserving average inter-element spacing. (Steering vectors become
eigenvectors of signal and the eigenvalues the actual density)
In consequence, the following estimates converge asymptotically to the actual
power distribution:
Q
POT 1 ( S )
POT 2 ( S )
S H .R m1 .S
S H .R m .S
1
m
2
H
m 1
q . S .e q
q 1
Q
2
H
m
q . S .e q
q 1
NOTE that MUSIC does
not have this property
since function f(x) is not
invertible. So Music is a
source detector but no an
estimate of the power
distribution versus angle
2
Q m H
q1 q . S .e q
2
Q
1
H
H
LOG
( S ) Ln exp S .R .S Ln exp( q ). S .e q
q 1 Miguel Angel
Array Processing
H
.R
September 13
.S
Lagunas DOA Estimation
40
MUSICforcolorednoisebackground
When the noise back-ground covariance is still full range and known a priori,
still MUSIC can be used with the following modification:
Assuming that the model for the measured
covariance is:
R S .P.S R o
H
NS
P( s) S s S R 0
Then--1/ 2
1/ 2
R 0 .R.R 0
being
1/ 2
bs R0 S s
September 13
NS
s 1
P( s).b s b s I
H
s
s 1
In consequence, the noise eigenvectors of matrix or
generaliced noise eigenvectors of the matrix pencil
R and R0 are orthogonal to the new vector bs
1/ 2
.R.R
1/ 2
Array Processing
Miguel Angel
0
0
Lagunas DOA Estimation
Re q (q) R 0 e q41
In summary the new estimate, will be constructed over the noise generalized
eigenvectors of the pencil matrix R, R0 and the sources will be located at the
local maxima, among S, of ..
(S )
1
2
q NS 1
b( S ) .u q
1
2
q NS 1
1 / 2
S .R 0 .u q
Note that the selection of the dimension stays on the generalized
problem as well as the difficulties on its estimation from data
covariance matrixes.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
42
ScanningcomplexityonDOA
estimation
For resolution close to 0.1, the scanning, trough vector S, in all the estimates
reported,
implies
the
computation
of
a
QxQ
quadratic
form
(360*10).(180*10)=6.480.000 times (!!!!!!).
This represents a hard limitation for implementing these procedures
on limited resource array stations.
- Limited field of view
- Selective scanning
- ESPRIT using moving platform or twin apertures.
The basic idea is to perform spatial linear prediction and solving it as an
exact problem using singular value decomposition SVD.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
43
ESPRIT
d
Xn
Yn
After target illumination, two set of snapshots, X
and Y are collected from two virtual apertures
separated a vector d, motivated by the constant
velocity movement of the platform.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
44
X n
Zn
Y n
Is the 2Q size global snapshot
Our target is to design a spatial linear predictor T such that the second
snapshot is predicted from the first minimizing the prediction error
T I Z
n T
X n
I n Yn T X n n
Y n
Regardless the procedure is coined as ESPRIT, the idea was original from J.
Munier (Grenoble) and it was named as Propagateur. The Propagateur was
used, not just for DOA estimation but also for callibration of towed sonar arrays.
Before solving the LP problem we will reduce
noise using the svd of the covariance of the
global snapshot.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
45
Using the svd of covariance Rz
1
Rz
N
q 1
H
n
E s s E s E n n E n
It is clear that the signal subspace
eigenvector can be decomposed
on the two apertures contributions
In addition, we know that the signal
eigenvectors of each aperture are
related to the corresponding DOAs by a
rotation matrix.
E x
Es
E y
E x SG
These two contributions are
QxNS size
September 13
Signal sub-space of
dimension NS (equal to the
actual number of sources
where
diag s 1, NS
Angel
Array Processing Miguel
Lagunas DOA Estimation
E y S G
2f c H
exp j
k s ( unitary ) d
c
46
THE PROBLEM NOW IS TO DESIGN THE LINEAR PREDICTION
BETWEEN THE TWO SUBSPACES. This can be formulated as finding F1
and F2 in this problem:
F1
E y .
F 2
To solve the problem let us take a look of svd of the composition of the two subspaces
E y U .D.V
H
V 11 V 12
H
V 21 V 22
September 13
2NS
Only the diagonal of the red
square contain the NS
Array Processing Miguel Angel
47
NS
Lagunas DOA Estimationeigenvalues different from zero
V 11 V 12
H
V 21 V 22
V 21
V
22
Multiplying by this matrix.
0
I
This product produces zero or minimum residual. In
consequence:
V 12
E y . 0
V 22
September 13
E x V 12 E y V 22
or
E x V 12 V 22 E y
E T Ey
x Angel
Array Processing Miguel
Lagunas DOA Estimation
48
E xT E y
Now, having T after 2 svd, we resort the relationship
between eigenvectors and steering vectors
E x SG
E y S G
. And we obtain
S GT S G
GT G
And, finally.
1
T G G
After 3 svd we get the NS DOAs from the eigenvalues
of matrix T. This completes the ESPRIT procedure
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
49
MAXIMUMLIKELIHOOD
PointSourceDOAestimation
First at all, let us solve the problem of a single source in white noise.
Pr X
an , S , R 0
1
exp X
det R 0
The ML estimate of the complex
envelope is:
a n
an S R 0
R0 X
R0 S
an S
X
Q
Note that to recover the envelope estimate we need the location or steering of the
source. Now using the previous estimate in the likelihood we have:
Pr X
S ,
September 13
2 Q
exp X
X
Q
Array Processing Miguel Angel
Lagunas DOA Estimation
X
Q
50
After some manipulations on
the exponent.
P is a projection operator.
P S
H
H
S Xt 1
S Xt
Xt S
Xt S
2
Q
Q
H
H
H
S S
SS
SS
1
H
H
X t 2 X t I
I
X t I
Q
Q
Q
P IS S
1
X t 2
This operator provides the projection on
the orthogonal subspace to the space
defined by matrix S
This P generates a projection of
the received snapshot on the
orthogonal subspace of S, i.e.
removes the content from
direction S from the data
snapshot.
Pv
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
SS
I
t
51
Since the trace of a scalar is the same scalar and the trace is circular then.
tr X
S
S
H
t I
Q
S
S
1
tr I
t
2
tr A B tr B A
tr a b scalar a b tr b a
H
X t X
H
t
To continue the estimation we need additional snapshots N. These snapshots are
independent ( since the noise is white on the time domain), and, in consequence
the probability will be the product of the N single probabilities.
Pr X
M
S, R0
N
Ln Pr X
1
2 QM
H
SS
exp I
Q
S , R 0 QNLn
1
2
M
X t X
1
H
t
tr N P R
being
September 13
1 N
SS
P I
y R
X
Array Processing Miguel Angel
Q Lagunas DOA Estimation
M 1
H
n
52
Now we can go for the estimate of the noise level just taking derivative of loglikelihood function
QM
2 2
tr M P R 0
The resulting estimate is
2
ML
tr P R / Q
Using the noise estimate at the log-likelihood
M
Ln Pr X
1
2
Q
S QMLn ML
Since the maximum among S, this implies that the proper S has to minimice
the noise power estimate.
2
MIN S tr P R / Q MIN
MIN S ML
S
S
R
1
tr R
Q
Q
H
SSH R
S RS
MAX S
MAX S tr
MAX
Processing Miguel Angel
Q ArrayLagunas
September 13
DOAQ
Estimation
1
QM
2
n
53
Whenthenumberofsourcesisgreaterthanone
complexitygrowsmakingalmostunpracticaltheML
estimate
X n S .a n w n
And the new
likelihood.
Pr( X / S , Q)
where
And the log-likelihood
L( S , ) k0 N .Q.Ln 2
2Q N
P I .S . S .S
n 0
.S . X n
.S
N 1
N 1 H
2
. exp X n .P .P . X n /
n 0
.P . X n
2
1 N 1
The ML estimate of
. P .X n
the noise
level is.
N .Q n 0
Array Processing Miguel Angel
54
2
September 13
a n S .S
The ML estimate
of the sources
waveforms is:
Lagunas DOA Estimation
At the end the function to be minimized is
quite similar to the case of a single source.
Nevertheless it implies a double search over
the to angles simultaneously. For NS sources
this implies a simultaneous search over NS
angles or 2NS angles for planar apertures
n 0
LS Traza P .R
N 1
L S P .X n
In addition, L(S) function is smooth which almost precludes any
accelerated search of local maxima. An example can be viewed
for two sources in an ULA array located at -21 and 10 degrees.
80
60
40
20
0
-20
-40
-60
-80
-80
September 13
-60
-40
-20
20
40
Array Processing Miguel Angel
Lagunas DOA Estimation
60
80
55
THEEMAlgorithm
PRELIMINAR: The MAP Estimate, the problem is given the statistics
(mean and covariance) of vectors Y and X, we like to have a MAP estimate
of Y given X.
The estimate is:
Y m y C .C 1 X m x
yx
xx
being : m y E Y and C
with covariance C C
yy
September 13
yx
yy
E Y m y X m x H
C C 1 C
yx
Array Processing Miguel Angel
Lagunas DOA Estimation
xx
xy
56
Using the model, we
have..
m x H .m y
Y m y C
H H C 1 X H m y
yy
yx
xx
When
Y my wy C
yy
yy
HH
2
y I Q. NS
xx
H .C
xx
x2
NS
Y m y H H C 1 X H m y
xx
September 13
yy
Array Processing Miguel Angel
Lagunas DOA Estimation
H H x2 I
H H H 2 I
57
TheEstimateandMaximice
Observed Data
Snapshot
Xn
UN-COMPLETE DATA
Desired Data
Single source snapshot
Yn
COMPLETE DATA
Modeling un-complete and complete data .
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
58
NS
X n a s n S n w n
s 1
NS
m x as n S n
s 1
C xx 2 I
Xn H
September 13
Y 1n a1 ( n) S 1 w1n
.. ..
Y n Y sn as ( n) S s w sn
.. ..
Y a ( n) S w
NS
NSn
NSn s
a1 ( n) S 1
..
m y as (n) S s
..
a (n) S
Y n I Q s .... NS I Q Y n
2
I Q. NS
C yy
NS
Array Processing Miguel Angel
Lagunas DOA Estimation
59
GoingbacktotheMLestimate
Let us group on vector all the parameters we like to estimate, i.e. the source
waveforms and the corresponding steering vectors, the likelihood we would like to
maximize is:
Pr n
In order to connect this with the complete data
we use Bayes theorem.
Pr n
Xn
Pr
Pr
Y
,
Pr Y n
X ,
n
Since, given the complete data the un-complete
data can be found directly, without the need of the
sources parameters, in consequence this term do
not impacts on the maximization
the likelihood
Arrayof
Processing
Miguel Angel
September 13
Lagunas DOA Estimation
60
In terms of the log-likelihood the function to be
maximized is:
X n Y n Y n
X
,
Let us assume we have a prior of the parameters (m) at step m this allows us
to have the mean of the un-complete data. Is to The objective is to obtain a
better estimate (m+1)
Y n
Y n
(m)
( m 1)
(m)
( m 1)
,
X n , ,
U
September 13
(m)
( m 1)
(m)
( m 1)
Array Processing Miguel Angel
Lagunas DOA Estimation
61
(m)
( m 1)
(m)
( m 1)
Since U entails the NS single source likelihood problems, solving these NS
problems-- MAXIMIZE we will find an improvement on U. In summary:
(m)
( m 1)
(m)
(m)
This is set in order to maximize U
Since V entails the MAP estimate of Y given X and the prior of the parameters,
whenever we use a different parameter vector from the prior the function will
decrease ESTIMATE. In summary:
(m)
( m 1)
(m)
(m)
ITERATING OVER ESTIMATE and MAXIMIZE STEPS WE WILL IMPROVE THE
MULTIPLE SOURCES LIKELIHOOD
solving the problem.
Array Processing Miguel Angel
September 13
Lagunas DOA Estimation
62
TheMAXIMIZEstep
This step reduces to solve the ML problem for a single source in white
noise.
- NOTE THAT the sources waveforms estimation implies to perform
EM steps at the single snapshot level. At the same time, noise in the
complex envelope may require several iterations in order to reduce
the effects of the estimation noise on the following estimate step.
- The need to estimate the source waveform is the major drawback of
the EM procedure.
- Convergence is fast for reduces number of sources
- The procedure is robust to coherent sources
- The weight 1/NS comes from assuming the noise on ion of the
complete data a fraction of the global noise. When it is assumed that
the noise on each complete data is equal to the global noise, the
coefficient reduces toArray
1. Processing Miguel Angel
September 13
Lagunas DOA Estimation
63
TheESTIMATEstep
Given the source parameters from the prior, we compose the mean of the
sources vector..
Then, we obtain the complete data as:
a1( m ) ( n) S 1( m )
( m 1)
..
1
(m)
H
Y n
my
H X n H m y
(m)
(m)
(m)
NS
m y as ( n) S s
..
Q
(m)
(m)
1
(
)
m
a NS ( n) S NS
my
.. X n H m y
NS
IQ
2
I Q. NS
C yy
NS
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
64
( m 1)
Yn
m (ym )
( m 1)
Y sn
I
Q
.. X n H m y
I Q
NS
Xn
a q( m ) S (qm )
q 1
qs
The EM becomes very intuitive: To generate the data for a single source the
rest of the sources contribution is subtracted (is removed) from the received
snapshot
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
65
EMArchitecture
Without loss of generality we can draw the processing architecture of the EM
algorithm for the case of two sources.
(m)
1n
a1(n), S1
Estimate s1
Xn
(m)
1
a
-
(m)
2
( n) S
(m)
2
Estimate s2
Y
September 13
(n) S
(m)
1
(m)
2 nArray Processing
Miguel Angel
Lagunas DOA Estimation
a2(n), S2
66
TheAlternateProjectionalgorithm
Some, less formal, versions of the EM has been used in practice in field like
radar (CLEAN) as well as sub-optimal procedures which overpass the source
waveform estimation like the AP. AP is a concept used on linear programing
methods an it is well known in mathematics.
The basic idea of the EM remains in the sense that passing from the original
snapshot to the single source problem we need to remove the other source effects.
Nevertheless, in order to remove the other sources we use blocking the spatial
directions instead of subtracting waveforms.
As an example, this matrix may be used to source s from a
snapshot
0
.
0
1 exp( j (u s 2 u s1 )
exp( j (u s 3 u s 2 )
1
.
0
.
.
.
.
0
.
exp( j (u sQ 1 u sQ 2 )
0
0
0
0
.
1
exp( j (u sQ u sQ 1 )
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
67
The AP architecture.
(m)
B2
(m)
1n
A Kalman tracker of the
source trajectory can be
added to further improve
the resulting performance
Estimate
DOA1
S1
Xn
(m)
B1
(m)
Y 2n
September 13
Estimate
DOA2
Array Processing Miguel Angel
Lagunas DOA Estimation
S2
68
The processing in every branch reduces to a blocking of the rest of the
sources followed by a phased array scanning as corresponds to the ML
estimate for a single source.
These two steps can be done simultaneously by a single
beamforming that steering the desired, nulls out the rest of the
sources.
Thus, the beamformer scanning for source s has to hold the
following constrain:
H
A . S, S
(m)
1
,.., S
(m)
s 1
,S
(m)
s 1
,.., S
(m)
NS
1,0,..,0
In addition, it has to shown minimum response to the white
noise, i.e. has to be minimum norm.
The new estimate of the DOA for
source S, results from the absolute
MIN
maxima of
( m 1)
H
Array Processing Miguel Angel s
x 69
S
Lagunas DOA Estimation
A A
September 13
max A .R . A
In general: For s=1:NS
C s S , S 1 ,.., S s 1 , S s 1 ,.., S NS
(m)
(m)
(m)
(m)
1 [1 0 .. 0]
T
1
H
A C s . C s .C s .1
( m 1)
s
A H .R . A 1H . C H .C 1.C H .R .C . C H .C
x
x
s
s
s
s
s
s
max
1
H
H
S
A .A
1
.
C
.C s .1
s
.1
end
Note that we divide by the noise bandwidth of the scanning beamformer to
form the DOA estimate
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
70
ExampleofAPperformance
Two sources. ULA array half wavelength.
The number of snapshots is ...4000
WAR source 1 is coherent with source 2
The number of sources is...2
Source #1--> 10dB 0 elevation 0
azimuth
The velocity in elevation is -0.005 /snap
Source #2--> 10dB -20 elevation 0
azimuth
The velocity in elevation is 0.005 /snap
The number of sensors is......8
Elevation Field of view from -25 up to 5
Scanning precission in elevation of 0.2
degrees
20
18
16
14
dB
12
10
8
6
4
2
0
-25
-20
September 13
-15
-10
-5
Elevation in degrees
Array Processing Miguel Angel
Lagunas DOA Estimation
71
Snapshot # 4000
proc.#1 doa -20 radial vel. -0.0049529
proc.#2 doa 0 radial vel. 0.0050821
Tracking plot (blue=actual)
5
DOA detected (red)
-5
-10
-15
-20
-25
September 13
50
100
150
200
scans
250
Array Processing Miguel Angel
Lagunas DOA Estimation
300
350
400
72
The number of snapshots is...3000
WAR source 1 is coherent with source 2
The number of sources is...2
Source #1--> 20dB 60 elevation 0 azimuth
Elevation velocity -0.015 /snapshot
Azimuth velocity -0.02 /snapshot
Source #2--> 20dB 70 elevation -130 azimuth
Elevation velocity -0.02 /snapshot
Azimuth velocity 0.03 /snapshot
The number of sensors is...13
Elevation Field of view from 5 up to 75
Scanning precission in elevation of 0.2 degrees
slowness/azimuth plot, actual blue, estimated green
90
120
1
0.8
60
0.6
150
30
0.4
0.2
180
330
210
240
300
270
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
73
Azimuth tracking plot
0
-20
-60
-80
-100
Elevation tracking plot
-120
70
-140
60
250
50
100
150
snapshots
200
300
50
elevation in degrees
azimuth in degrees
-40
40
30
20
10
September 13
0
Array Processing
Miguel
0
50 Angel 100
Lagunas DOA Estimation
150
snapshots
200
250
74300
APforBeamforming
Our goal: A beamforming steeresd to some angle (0)
With perfect nulling of unknown coherent or uncoherent interferences
Initial beamforming: Phased Array steered at the desired
The number of snapshots is...5000
WAR source 1 is coherent with source 4
The number of sources is...4
Source #1--> 10dB 0 elevation
Source #2--> 20dB 70 elevation
Source #3--> 20dB -20 elevation
Source #4--> 30dB 5 elevation
The number of sensors is...16
Elevation scanning data
Elevation of the desired 0
quiescent response
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
75
Optimum as rx-1*sd
Array factor dB.
-5
-10
-15
-20
-80
September 13
-60
-40
-20
0
Elevation
20
40
Array Processing Miguel Angel
Lagunas DOA Estimation
60
80
76
Beamformer steered to the desired blocking s1
Beam response. Blocking..1 directions
-5
Array factor dB.
-10
-15
-20
-25
-30
September 13
-80
-60
-40
-20
0
Elevation
20
Array Processing Miguel Angel
Lagunas DOA Estimation
40
60
80
77
Periodogram response. Blocking..4 directions
30
-5
0
0
10
Array
Arrayfactor
factordB.
dB.
Array factor dB.
-10
-5
-5
5
-10
-15
-10
0
-15
-15
-20
-5
Scan blocking desired and directions before->
Find absolute maxima
25
The
number
of snapshots
is...5000
The design beamformer steered to the desired, blocking
interferers
with
20
The number
of sources
is...4
Periodogram
after desired
Blocking..1
blocked
Periodogram
response.
Blocking..2
directions
minimum norm.
Periodogram
response.
Blocking..3
directions
50
35
35
30
Source
#1--> 10dB 0 elevation
15
CONTINUE until flat scanning
Source
#2--> 20dB 70 elevation
30
30
25
10
Beam response. Blocking..4 directions40
Source #3--> 20dB -20 elevation
25
25
20
5
Source
#4--> 30dB 5 elevation
30
20
20
15
The
number of sensors is...16
0
Elevation scanning data
15
15
10
-5
20
Beam
Quiescent
Blocking..1
Elevation
of the desired 0
Beam response.
response.
Blocking..2 directions
Blocking..3
directions
10
10
5
-10
Remain 13 degrees of freedom
10
55
0
-15
Remain
12
degrees
of 0freedom
-100
-80
-60
-40
-20
20
40
60
Elevation
degrees
00
Remain
11 degrees of freedom
-5
0
Remain
10 degrees of freedom
-5
-5
-10
Desired level
-10
-10
-15
Estimate...10.0591
-100 -80
-80
-60
-40
-20 dB.
20
40
60
-100
-60
-40
-20
00
20
40
60
Elevation
Elevation degrees
degrees
Actual.....10 dB.
80
100
80
80
100
100
-20
-20
-10
-25
-25
-25
-15
-80
September
13-40
-80
-60 -60
-40
-80 -60
-20
-20
-40
0
0-20
Elevation
Elevation
20
20
Array Processing
Miguel Angel
60
80
40
60
80 40
0 40
20
60
Lagunas
DOA Estimation
Elevation
80
78
WIDEBANDDOAESTIMATION
The estimate depending on DOAs and frequency is obtained by defining a
steering vector for the wideband beamformer (Q sensor and NT tap delays).
S S , S exp j 2fT ,..., S exp j 2 ( N 1) fT
T
where
sin
S exp j 2 . f
d cos . exp j 2mfT
c
The beamformer output is:
y ( n) a . X X
H
T
n
T
n 1
..... X
Thus, replacing the steering and the new sanpshot, all the procedures
previously can be extended
to the wideband case
Array Processing Miguel Angel
September 13
Lagunas DOA Estimation
T
T
n N 1
79
FocusingonWidebandestimation
All the wideband estimates needs of focusing the source. The reason is that,
given a source, on position (,), that impinges the aperture in a bandwidth B,
in order to properly detecting the source we need to obtain the power that the
source produces at the array output. This power is computed as the sum of
the powers measured on the bandwidt B.
( , , f )df
B
This averaging of all the map is also known as focusing (resume in a single
plot) the power per angle solid
September 13
frequency
Array Processing Miguel Angel
Lagunas DOA Estimation
Once power is focalized,
super-resolution methods
can be applied.
80
FocusingforULA
For an ULA array of Q sensors and M delay lags, the output of a wideband
phased array will be:
M
X (u , f , n) X n . exp( j 2fmT ). exp( j 2q.u )
m 1 q 1
Thus, in this case, the output is like a 2D-DFT one in the spatial domain with
frequency u and the other in the time domain with frequency fT
Let us imagine that we have in the scenario a single source located at elevation
s , located at f0 with bandwidth almost zero.
The 2D-DFT will show a delta function locates at spatial frequency equal to us
and frequency f0T. The next figure depicts the situation in the spatial time
frequency plot.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
81
fT=0.5
us
f0T
u=-0.5
u=0.5
fT=0
But, THERE IS A BASIC DIFFERENCE between wideband DOA detection and
traditional 2D spectral estimation namely THE TWO FREQUENCIES ARE
COUPLED
f .d
. sen( )
c
September 13
and
c
c
f
.u 90 .u
d .sen( )
d
f .d
u
Array Processing Miguel Angel
82
c
Lagunas DOA Estimation
fT=0.5
us
DARK
ZONE
DARK
f0T
ZONE
u=-0.5
u=0.5
fT=0
When the source extents its bandwidth the u
frequency moves along a line
u dfT / cT sin( )
September 13
In consequence, the
power have to be
focused along lines
with slop equal to:
Array Processing Miguel Angel
Lagunas DOA Estimation
f .d
u
c
d
sin
c
83
fT=0.5
=-90
=90
f0T
u=-0.5
u=0.5
fT=0
For a central frequency fc, the sampling frequency as T=1/2fmax. Selecting
the sensor separation d as d=c/2, we have:
f max
2 f max c
d
sin( )
sin( )
sin
cT
f c Angel
2c Array Processing Miguel
September 13
Lagunas DOA Estimation
84
Note that for a given response, the resolution is always lower at low
frequencies than at high frequencies
fT=0.5
=-90
=90
f0T
u=-0.5
u=0.5
fT=0
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
85
f
1 /2 T
f= c /(2 d )
- 0 .5
f
1 /2 T
- 0 .5
Top: Corect choice of separation d in order to have maximum field of view. Also
for sources at the maximum vertical focussing is a good approsimation. Botton
Wrong choice of the interelement separation (large d)
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
86
Focusingforsuperresolution
Let us imagine that we have the array covariance matrix at a given frequency f as:
R ( f ) S f P f .S ( f ) ( f ) I
H
Focusing this matrix at frequency fc implies to find out a transformation such that:
T R ( f )T
S f c P f .S ( f c ) ( f ) I
H
With this, summing up all the matrix the power (also the noise) is focused in a single
matrix. Using Music NMLM or any other super-resolution method will provide
accurate source location
T R( f )T
f
September 13
2
S f c P f .S ( f c ) ( f ) I
f
Array Processing Miguel Angel
Lagunas DOA Estimation
87
Back to the design of the proper transformation, we can see that there are two
design conditions:
T R ( f )T
S f c P f .S ( f c ) 2 ( f ) I
T S ( f ) S fc
or T S f S f
TT
Clearly using the pseudo-inverse is not a valid solution
T Sf Sf Sf
c
S
1
H
f
since T T S f S f S f
H
S
1
We need to relax the first constrain which is better an more easy
than to relax the white Array
spatial
noise one.
Processing Miguel Angel
September 13
Lagunas DOA Estimation
H
fc
88
The design will be:
TS f S f
T S
Taking derivatives
Sf
s.t.
c
TT
L S
fc
Sf
S f S f U DfU
H
With the svd of the DOAs matrixes. Note
That for matrix L this merely instrumental
Since L is unknown.
S f S f V D ff U
H
L U DlU
TU D f Dl U
September 13
V D ff U
Clearly, the
solution is:
Array Processing Miguel Angel
Lagunas DOA Estimation
then
T VU
D l D f D ff
89
EJERCICIO
Una manera diferente para derivar metodos de
estimacion de DOA consiste en metodos
denominados de covariance matching, es decir, de
ajuste de la matriz de covarianza observada.
Sea la matriz obtenida a partir de un conjunto de N
snapshots, siendo N suficiente como para estabilizar
el estimador.
Es claro que el contenido de potencia
que proviene de una direccion dada
viene dado por
X m X mH
m 1
R S .S H
c
Asi pues para saber que contenido de potencia
tiene nuestra observacion R en una direccion
dada, basta encontrar el valor de que mas
asemeja la observacion R con Rc. Es decir el
que minimiza la diferencia entre las dos
matrices
September 13
1
R
N
Array Processing Miguel Angel
Lagunas DOA Estimation
P S min f R, R
90
PRIMER CASO: f(,) es la norma de Frobenius de la matriz diferencia.
RR
Nota:
c F
R S .S H
Encontrar el estimador de la
potencia, para cada direccion
que minimiza esta distancia
traza A H . A
traza A H . A A H A
SEGUNDO CASO: Si a la observada le restamos la potencia que le llega de
una direccion, la maxima potencia que podemos restar seria la que es aquella
que no hace la matriz diferencia definida negativa. Esta matriz diferencia
tendria al menos un autovalor cero, por lo que la nueva definicion de f(,) para
encontrar seria:
P S max min R S S H 0
Encontrar el estimador de la potencia, para cada direccion como el valor
maximo que provica en la diferencia un autovalor cero.
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
91
September 13
Array Processing Miguel Angel
Lagunas DOA Estimation
92
