FA-';’
OA T t.
High Resolution Algorithms for Spectral Analysis and
A C C E P T S ' 0
Array Processing
JLTY - /! ■ E : ’r > v
-by
e / Z , k y - “ L ,N W e i x i u D u
/
j
B.S~., Slianclong Polytechnical University, China, 1983
A Dissertation Submitted in Partial Fulfillment of the
Requirements for the Degree of
DOCTOR OF PHILOSOPHY
in the Department of Electrical and Computer Engineering
We accept th is d isse rta tio n as conform ing to th e req u ired s ta n d a rd
D r. R. L v n n / J ^ l i n . S upervisor (D e p a rtm e n t of E C E )
Dr. A. Z i e f i n s D e p a r t m e n t a l M em ber (D e p a rtm e n t of EC E)
Dr. P. Driesisen, D e p a rtm e n ta l M em ber (D e p a rtm e n t of EC E)
D r. It. R. D avidson. O u tsid e M em ber fD e o a rtm e n t of M a th e m atic s)
Dr. R. V&ifldiecl£ A d ^ itio n a^ M ertib er (D e p a rtm e n t of E C E )
Dr. A. R . M o g h a d d ag ^ o o d E x tern aH E x am in er (U. of W isconsin)
© WEIXIU DU, 1992
University of Victoria
All rig h ts reserved. D issertatio n m ay n o t be rep ro d u ced in w hole or in pari by photocopying or o th e r m eans, w ith o u t th e perm ission of th e a u th o r.
A B S T R A C T
In th is d isse rta tio n a novel covariance m a trix e s tim a to r has been proposed for th e covarian ce-m atrix b ased high resolution s p e c tra l e stim a to rs. T h e proposed covariance m a trix e s tim a to r can fully exploit cross correlations e x iste n t am ong d istin ct sets of ran d o m vectors draw n from different processes to o b ta in a m ore stable e stim a te of th e covariance m a trix from a s h o rt d a ta record c o rru p te d by ad d itiv e noise. T h is e s tim a to r is derived from a th eo re m on th e L east Squares Linear P re d ic tio n of one ran d o m v ector from a n o th e r ran d o m vector. T h e th eo re m can also b e in te rp re te d as e s tim a tin g th e auto-covariance m a trix of th e first ran d o m vector from th e th e cross-correlation m a trix b etw een th e tw o ran d o m vectors a n d th e auto-covariance m a trix of th e second ra n d o m v ecto r such t h a t a given o p tim al criterion is satisfied. A pplying th is m e th o d in c o n ju n ctio n w ith a high resolution alg o rith m results in perfo rm an ce im p ro v em en t of th e sp e c tra l e stim a to r.
T h e new covariance m a trix e s tim a to r has been ap p lie d in th e following th re e areas:
1. S p a tia l sm oothing for th e d irectio n of arriv al e s tim a tio n in th e presence of coherent signals.
2. Covariance en h a n ce m en t tiy u tilizin g te m p o ra l co rrelations b etw een a rra y sn ap sh o t vectors.
3. S p ectral e stim a tio n for tim e sequences.
S im ulations show t h a t th e e x p e cte d perfo rm an ce im p ro v em en t can b e achieved in term s of resolution, e s tim a tio n errors a n d S N R th reshold.
In a d d itio n to th e covariance m a trix e stim a to r, we also p resen t som e o th e r research results in array processing a n d seism ic signal processing. A general tra n sfo rm a tio n
m a trix based on th e vector p -n o rm has b een proposed. T his new t ranslorm at i o n
m a trix provides options to satisfy different design specifications in array processing. F in ally th e velocity e stim a tio n problem in seismic signal processing is discussed. T h e conventional sem blance m e th o d is found to be th e conventional beam lorm iiig m e th o d for a fixed two-way tim e. A n o p tim a l velocity e stim a to r is proposed based on th e L inearly C o n strain ed M inim um V ariance beam form er. T h e o p tim al v e l o c i t y
e s tim a to r d e m o n stra tes th e high discrim in atio n pow er (resolution) both f o r noise free d a ta and noisy d a ta . W hen th e new covariance e stim a to r is used in conjunct i o n
w ith th e o p tim a l velocity e stim a to r, we can achieve a resolution w ith deeper notch. T his fac t once m ore d e m o n stra te s th e advantages of th e proposed covariance m a t ri \
e stim a to r.
D r. R. L ynp jJyiriin, S upervisor (D e p a rtm e n t of FCK)
D r. A. ^SdKnski, D e p a rtm e n ta l M em ber (D e p a rtm e n t ol I'd K]
D r. P. Driess4n, D e p a rtm e n ta l M em ber (D e p a rtm e n t of KCK)
D r. R. R. D avidson, O u tsid e M em ber (D e p a rtm e n t of Mat ]lemat
D r. R. VahlHTeck, A dditional M em ber (D e p a rtm e n t ol hiCK)
T it le P a g e i
A b stract ii
Table o f C ontents iv
L ist o f P rin cip le S ym bols and U ses xi
List o f Tables x iv
List o f Figures x v
A cknow ledgm ents x v iii
1 In trod u ction 1
1.1 S p ectral E stim a tio n a n d A rray P r o c e s s i n g ... 1
1.1.1 A pplicatio n s a n d P erfo rm an ce M e a s u r e m e n t ... 2
1.2 Role of th e C ovariance M a t r i x ... 3
C O N T E N T S
1.3.1 S am ple C ovariance M a t r i x ... 1.3.2 M axim um Likelihood E s tim a te of th e S tru c tu re d C ovariance
M a t r i x ... 1.3.3 P erio d icity C on strain ed E s t im a t o r ... 1.3.4 C o h erent Signal-subspace M e t h o d ... 1.3.5 S teered C ovariance M a t r i x ... 1.4 C o n trib u tio n s of th e D issertatio n ... 1.5 O rg an izatio n of th e D i s s e r t a t i o n ...
2 H isto ric R ev iew
2.1 W hy High R esolution S p ectral E s t i m a t o r s ... 2.1.1 C lassical S pectral E stim a to rs ... 2.1.2 L im ita tio n of Classical S p e c tra l E s t im a t o r s ... 2.2 D a ta M odeling and A ssum ptions ... 2.2.1 A rray P r o c e s s i n g ... 2.2.2 S p e c tra l E s t im a t i o n ... 2.2.3 S u m m a ry of A s s u m p t i o n s ... 2.2.4 Sam ple C ovariance M a t r i x ... 2.3 High R esolution S p e c tra l E s t i m a t o r s ... 2.3.1 C onventional B e a m fo rm in g ... 2.3.2 M axim um E n tro p y M eth o d ( M E M ) ...
2.3.3 M inim um V ariance M ethod (M V ) ... 18 2.3.4 P isarenko M e th o d ... 18 2.3.5 M U SIC ... 19 2.3.6 R o o t M U S I C ... 20 2.3.7 M inim um N orm M e t h o d ... 21 2.3.8 E S P R IT a n d M atrix P encil ... 22 2.3.9 M axim um Likelihood (M L) M e t h o d ... 23 2.4 C lassifications of S p ectral E s t i m a t o r s ... 25 2.4.1 C lassification by C h ro n o lo g y ... 26 2.4.2 C lassification by M o d e l ... 26 2.4.3 C lassification by V ector S p a c e ... 27 2.4.4 C lassification by N u m erical P ro c ed u re ... 28 2.5 D e te rm in a tio n of th e N u m b er of Signals ... 28
3 Covariance E stim a tio n 30 3.1 In tro d u c tio n and S ta tistic a l B a c k g r o u n d ... 30
3.1.1 N uisance P a r a m e t e r s ... 32
3.1.2 E stim a tio n Im provem ent by A veraging M u ltip le E stim a te s . 33 3.2 LSLP of R andom V ectors ... 34
3.3 Im p ro v in g the E stim a te s of C ovariance M atrices by In co rp o ratin g C ro s s -c o r re la tio n s ... 36
C O N T E N j vii
4 A p p lication s 42
4.1 S p atial Sm oothing ... 12
4.1.1 I n tr o d u c ti o n ... 42
4.1.2 S p a tia l Sm oothing for th e Subspace-based M ethod ... 44
4.1.3 Im proved S p a tia l Sm oothing M e th o d s ... 15
4.1.4 S im u latio n R e s u l t s ... 53
4.2 C ovariance Im provem ent by T aking A dvantage of T em poral C orre lations ... 55
4.2.1 I n t r o d u c ti o n ... 5f 4.2.2 E s tim a tin g C ovariance M atrices by E x p lo itin g T em poral ( Cor rela tio n s ... 57
4.2.3 M a trix E n h a n c em e n t U sing T em p o ral C o rrelations ... 5!)
4.2.4 E n h a n c em e n t for C oherent Signals ... (JO 4.2.5 E x tension to W ideband S i g n a l s ... Cl 4.2.6 S im u latio n R e s u l t s ... 62
4.3 S p ectral E stim a tio n ... 64
4.3.1 I n t r o d u c ti o n ... 61
4.3.2 M a trix Telescope S e rie s ... 67
4.3.3 P ro p o sed M e t h o d ... 67
4.3.4 O rd e r S e l e c t i o n ... 70
5 D esig n o f Transform ation M atrices 76
5.1 I n tr o d u c ti o n ... 76
5.2 B a c k g r o u n d ... 78
5.2.1 V ector P -n o rm and Frobenius N o r m ... 78
5.2.2 Sector Processing a n d In te rp o la te d A r r a y s ... 79
5.3 P -norm M e t h o d ... 81
5.3.1 I m p le m e n ta tio n ... 82
5.4 Sim ulation R e s u l t s ... 84
6 S eism ic Signal P rocessin g 86 6.1 I n tr o d u c ti o n ... 86
6.2 M odeling Seism ic Signals ... 87
6.2.1 W avefront M o d e ls ... 87
6.2.2 W ave S hape M o d e l ... 89
6.3 V elocity E s t i m a t i o n ... 90
6.3.1 C u rre n t T e c h n iq u e s ... 90
6.3.2 M odifications and Im provem ents O ver th e P revious M eth o d s 95 6.4 T h e R elationship B etw een Sem blance an d E ig e n stru c tu re V elocity E s t im a t o r s ... 97
6.4.1 Sem blance and C onventional B e a m fo r m in g ... 97
6.4.2 O p tim a l B e a m f o r m e r s ... 98
C O N T E N T S ix
6.5.1 Som e Special C o n s id e ra tio n s ... 99
6.5.2 E n h a n c em e n t of th e E s tim a te s of C ovariance M atrices . . . 100
6.5.3 L C M V V elocity E s t i m a t o r ...101
6.6 S i m u l a t i o n s ...103
6.6.1 C om parison of C oherency M easure T h resh o ld D iscrim ination 103 6.6.2 H andling N oisy Seism ic D a ta ...105
7 Su m m ary and C onclusions 122 7.1 I n tr o d u c ti o n ... 122 7.2 S u m m a ry of th e D i s s e r ta tio n ...122 7.3 C onclusions and F u tu re W o r k ... 124 A M a tr ix D erivative R u les 126 A .l D erivation of V e c t o r s ... 126 A. 1.1 D e f in itio n s ... 126 A .1.2 C h ain R u l e ... 127 A .1.3 Som e Useful F o r m u l a s ...127 A .2 D erivation of M a t r i c e s ...127 A .2.1 D e f in itio n s ... 128 A .2.2 P ro d u c t a n d C hain R ules ...128
A.2.3 Som e Useful F o r m u l a s ...129
B .l A n o th er Form for th e F isher Inform ation M a t r i x ... 130 B.2 Inform ation M a trix for G aussian R an d o m V e c t o r s ...130 B.3 In fo rm atio n M a trix for A rray P r o c e s s i n g ... 132
L I S T O F S Y M B O L S
List of Principle Symbols and Uses
a
90; — p a rtia l derivative w .r.t. 6-II' I If — m a trix Frobenius norm
II • Up — v ecto r p-norm
( .)# — pseudo-inverse of a m a trix
(•r1
— inverse of a m a trixor
— com plex conjugate(•)* — H e rm itia n tra n sp o se of a m a trix
Sij — K ronecker d e lta
A — sensor spacing
A R i — p e rtu rb a tio n of Ri
X — w avelength
A; — ith eigenvalue of covariance m a trix R
A — eigenvalue m a trix
— m ea n vector — K ronecker p ro d u ct
u> — frequency in ra d
a 2 — noise variance
Oi —■ b e a rin g angle of th e «th signal source <p,6 — p a ra m e te r to b e e s tim a te d
a(0) — ste e rin g vector
m — ste e rin g m a trix
A R — auto-regressive
A R M A — a u to regressive-m oving average
B — ste e rin g m a trix for th e v irtu a l array
c, V — w ave velocity of th e m ed iu m C D P — com m on d e p th p o in t
CSM — co h erent signal-subspace m eth o d
d — n u m b e r of signals
d ( x , t ) — m ea su red seism cg ram det[-] — d e te rm in a te of a m a trix D F T — d isc re te F ourier tra n sfo rm DOA — d irectio n of arrivals
E S P R IT — E s tim a tio n of signal p a ra m e te r via ro ta tio n invariance tec h n iq u e
fv — p eak frequency of th e R icker wave sp e ctru m
F F T — fast F ourier tra n sfo rm F IR — finite im pulse response
/ id e n tity m a trix
...
0 M x k , I m, 0m x(A T -M -A :)1
Im,n( k) ... I j Oj x M O jx (A T -M -j)
0 ( N - M - j ) x j 0 ( N - M - j ) x M i N - M - j
J (t d) ^isher inform ation m a trix as a fu n ctio n 0 ’ <PD Jr — reflection or exchange m a trix
Kc — coherency m easu re by K ey’s m eth o d
LCM V — L inearly con strain ed m in im u m variance b eam fo rm er
LP — linear p rediction
LS — least squares
M —- n u m b er of sensors
MA m oving-average
MEM — M axim um e n tro p y m eth o d
ML — M axim um likelihood
MV — M inim um variance
M SE — m ean squared erro r
n ( x , t ) — m ea su rem e n t noise in seism ogram P — root m ea n squared slowness P — signal covariance m a trix
Pa — p ro je c tio n o p e ra to r of m a trix A
P(u>) — pow er sp e c tru m density rank[-] — ran k of a m a trix
Vi — ith ro o t of a polynom ial
r ( k ) —- au to c o rre latio n function
II — covariance m a trix
R { t o , v a ) — in -g ate covariance m a trix
I P { t o , V 3) — im proved in -g ate covariance m a trix (squared) Ra — covariance m a trix of x a { t )
R j — forw ard sp a tia lly sm o othed covariance m a trix
Rfb — - forw ard-backw ard sp a tia lly sm o o th e d covariance m a trix k ) —- forw ard sp a tia lly sm o othed covariance m a trix
(squared) w ith im provem ent
R } b forw ard-backw ard sp atially sm o o th e d covariance
m a trix (squared) w ith im provem ent
R i3, R i j — forw ard cross-covariance m a trix of X { ( t ) a n d X j { t ) — backw ard cross-covariance m a trix of X j { t ) a n d X j ( t ) Ru — im proved e stim a te of Ru
Rh — es' 'm a te of R u from cross co rrelatio n s
IP-ti ( s tim a te of R u (squared) fro m cross co rrelations
L I S T O F S Y M B O L S xiii
s — signal vector
S i { t ) — f th in cid en t signal
span{-} — v ecto r space span ?d by colum ns of a m at rix
Sc — sem blance coefficients
Senhanced — visually enhanced sem blance coefficients
S — signal m a trix
SN R — signal to noise ra tio
to — two-way tra v e l tim e for zero-offset tra c e tx — two-way tra v e l tim e for a tra c e of offset x t r ( ' ) — tra c e of a m a trix
T — lin e a r tra n sfo rm a tio n m a trix
U ( p , t 0) — hy p erb o lic tra n sfo rm coefficients
U,v
— singular vector m a trixVi — i t h eigenvector of th e covariance m a trix II. Vrms — ro o t m ean squared velocity
vs
— signal eigenvector m a trixVn — noise eigenvector m a trix
w
— noise m a trixX — tra c e offset
x { t ) , y ( t ) — observed d a ta (vector) 3 'a i.t) — a u g m e r te d sam ple vector
X
-— d a ta m a trix6.1 E vents p a ra m e te rs for noise-free d a t a ... 105 6.2 E vents p a ra m e te rs for noisy d a t a ... 115
L ist o f F ig u res
4.1 T ran slatio n al equivalent s u b a r r a y s ... 43 4.2 P a rtitio n in g a uniform lin ear a rra y in to s u b a r r a y s ... 49 4.3 M U SIC sp e ctra using th e conventional and th e im proved forw
ard-b ack w ard sp a tia l sm oothing m e t h o d s ... 54 4.4 R oot M S E vs. SN R for th e conventional a n d th e im proved sp a tia l
sm o o th in g m e t h o d s ... 56 4.5 F o rm a tio n of au g m e n ted sam ple v e c to r s ... 58 4.6 M U SIC sp e c tra using th e conventional and th e enhanced covariance
m a trix e stim a tio n m eth o d s (using te m p o ra l co .- r e la tio n s ) ... 63 4.7 R oot M S E vs. SN R for th e conventional a n d th e en hanced covari
ance m a tr ix e stim a tio n m eth o d s (using te m p o ra l correlations) . . . 65 4.8 R oot M S E vs. SN R for a u g m e n ted sam ple vectors of different sizes 66 4.9 T elescoping series of a m a t r i x ... 68 4.10 P isare n k o s p e c tra using th e conventional a n d th e new covariance
m a tr ix e stim a to rs (frequency e stim a tio n ) ... 73 4.11 R oot M S E vs. SN R for th e conventional a n d th e new covariance
e s tim a to rs (frequency e s t i m a t i o n ) ... 74
4.12 T h e average R oot M SE vs. order of th e covariance m a tr ix ... 75
5.1 C onvergence of th e N ew ton’s m e t h o d ... 85
6.1 A m odel for th e reflection seism ology s y s t e m ... 88
6.2 An h y p o th e tic a l C D P g a t h e r ... 91
6.3 Form ing a tim e g a te for sem blance a n a l y s i s ... 93
6.4 N oise-free s y n th e tic C D P g a th e r for d iscrim in atio n th re sh o ld te s t . 106 6.5 V elocity contour using sem blance clean d a t a ) ...107
6.6 V elocity con to u r for LCM V w ith sp a tia l sm o o th in g (clean d a ta ) . . 107
6.7 M arginal p lo t for region 1 (clean d a t a ) ... 108
6.8 M arginal p lo t for region 2 (clean d a t a ) ... 109
6.9 M arginal p lo t for region 3 (clean d a t a ) ... 110
6.10 V elocity contour using th e disp lay -en h an ced sem blance (clean d a ta ) 111 6.11 M arginal p lo t of th e display-enhanced sem blance velocity sp e c tru m for region 1 (clean d a t a ) ... I l l 6.12 M arginal p lo t of th e d isp lay-enhanced sem blance velocity sp e c tru m for region 2 (clean d a t a ) ... 112
6.13 M arginal p lo t of th e display-enhanced sem blance velocity sp e c tru m for region 3 (clean d a t a ) ... 112
6.14 Noisy sy n th e tic C D P g a t h e r ...115
L I S T O F F I G U R E S xvii
6.16 V elocity contour for K ey’s m eth o d (noisy d a ta ) ... 116 6.17 V elocity co n tour for LCM V w ith sp a tia l sm oothing (noisy d a ta ) . . 117 6.18 V elocity con to u r for LC M V w ith th e im proved sp atial sm oothing
(noisy d a t a ) ... 117 6.19 M arginal p lot of th e sem blance velocity sp e ctru m in region 3 (noisy
d a t a ) ... 118 6.20 M arginal p lo t of th e velocity sp e c tru m by th e LCM V w ith th e con
ven tio n al sp a tia l sm o o th in g in region 3 (noisy d a t a ) ... 118 6.21 M arginal p lo t of th e velocity sp e c tru m by th e LC M V w ith th e im
proved sp a tia l sm o o th in g in region 3 (noisy a t a ) ...119 6.22 V elocity con to u r for th e d isp lay-enhanced nem blance (noisy d a ta ) . 119 6.23 M arginal p lo t of th e display-enhanced sem blance velocity sp e ctru m
for region 3 (noisy d a t a ) ... 120 6.24 O v erlaid m arg in al plot of th e sem blance velocity s p e c tra in region
3 (noisy d a t a ) ... 120 0.25 O verlaid m arg in al plots of th e velocity s p e c tra by th e L C M V w ith
th e conventional sp a tia l sm o o th in g in region 3 (noisy d a t a ) ... 121 6.26 O verlaid m arg in al plots of th e velocity s p e c tra by th e L C M V w ith
A c k n o w le d g m e n ts
I wish to express m y sincere g ra titu d e to m y supervisor, Dr. R . Lynn K irlin for his c o n sta n t en co u rag em en t, guidance, su p p o rt a n d assistan ce th ro u g h o u t th e d isse rta tio n research. His deep insig h t in a rra y signal processing a n d invaluable suggestions m ad e th is research all possible.
I would also like to th a n k Drs. A. Zielinski, P. D riessen, R. R . D avidson and R. V ahldieck for th e ir services as m y supervisory c o m m itte e m em b ers, a n d Dr. A. R. M oghaddam joo for his service as th e e x te rn a l ex a m in er in m y P h. D. oral ex am ination.
M y acknow ledgm ents a re also d u e to all th e facu lty m em b e rs a n d staff in th e E C E d e p a rtm e n t. T h e help th e y offered to m e d u rin g m y P h . D. p ro g ra m is deeply ap p reciated .
T h is d isse rta tio n is d e d ic a te d to m y p a re n ts, M r. X ingw u Du a n d M rs. Y ishu Z hang, m y wife Jessica, a n d m y son M elvin. T h e ir love a n d care m a d e it w o rth while.
T h is research is p a rtly su p p o rte d by th e U niversity of V ic to ria Fellow ship, an d th e P e tch R esearch Fellowship.
C h a p te r 1
I n tr o d u c tio n
1.1
S p e c tr a l E s t im a t io n a n d A r r a y P r o c e s s in g
S p e c tra l e stim a tio n or sp e c tra l analysis has for several decades been a tra d itio n a l research area for sta tistic ia n s. T h e in tro d u c tio n of th e fast Fourier tra n sfo rm ( F F T ) a lg o rith m two decades ago, however, has ex p an d ed th e role of spectral e s tim a tio n from one of research novelty to one of p ra c tic a l u tility. T h e digital signal processing c o m m u n ity has, as a resu lt, ta k e n an ever-increasing in te rest in spectral e s tim a tio n research an d ap p licatio n s. In th e signal processing co m m u n ity a sp ec tra l e s tim a to r is referred to as any signal processing m eth o d which characterizes th e frequency c o n te n t of a m easu red signal. T h e p rac tic a l problem s in sp ectral analysis are o fte n to o b ta in an e stim a te of th e tru e sp e c tru m from a finite interval of o b servation w hich is usu ally c o rru p ted by a d d itiv e noise.
A rra y processing deals w ith th e processing of signals c arried by p ro p ag a tin g wave p h e n o m e n a . T h e received signal is o b ta in e d by a sensor array located in th e field of in te re st. T h e aim of a rra y processing is from th e sensed d a ta to e x tra c t useful in fo rm a tio n a b o u t th e received signal field (e.g. its sig n a tu re , directio n of a rriv al (D O A ), speed of pro p ag atio n ). A lthough sp e ctral e s tim a tio n and a rra y processing a p p e a r to b e d istin c t, m an y of th e alg o rith m s proposed are equally a p
plicable to b o th problem s [1], For th is reason, two problem s can often b e discussed in a unified approach.
Several recen t books by Kay [2] M arple [3], H aykin [4, 5] provide an excellent review on these topics.
1.1.1
A p p lic a tio n s and P er fo r m a n ce M ea su rem en t
Spectral analysis a n d a rra y processing have a w ide range of ap p licatio n s such as system id en tificatio n , ex p lo ratio n seism ology, sonar, ra d a r, rad io a stro n o m y a n d tom ography. T h e m ea su rem e n t of th e p erfo rm an ce for a sp e c tra l e s tim a to r depends on ap p licatio n s. Som e useful sta tis tic a l m ea su rem e n ts of s p e c tra l e stim a to rs a re bias, e rro r covariance, resolution a n d signal to noise (SN R ) th resh o ld . L et p 6 C dxl be th e tr u e p a ra m e te r v ecto r a n d p be th e e s tim a te o b ta in e d using a sp e c tra l e stim a to r from a finite observation c o rru p te d w ith noise. T h e n we define bias, e rro r covariance , reso lu tio n a n d SN R th re sh o ld as follows, w ith E d e n o tin g ensem ble or e x p ected value a n d th e su p e rsc rip t H d e n o tin g m a trix H e rm itia n tran sp o se.
D efin ition 1.1 (B ias): The bias o f an estim ate by a spectral estim a to r is defined
by the sta tistica l average o f the difference between the estim ated pa ra m eter vector p and the true p aram eter vector ip or
A p = E [ p - tp) = E[p] - <p. (1.1)
D efinition 1.2 (Error Covariance): The covariance o f the estim a tio n error
associated with a spectral estim a to r is defined by the covariance o f the error vector — ip or
C H A P T E R 1. I N T R O D U C T I O N
For the unbiased estim ator, w t have
Cov{(p) = E { { f > - < p ) ( $ ~ v ) U} (1.3)
D e f i n i t i o n 1 .3 ( R e s o l u t i o n ) : The capability o f the spectral estim a to r to resolve
very close frequency com ponents is referred to as resolution, which is com m only quantified with its probability.
D e f i n i t i o n 1 .4 ( S N R T h r e s h o l d ) : The S N R above which the spectral estim ator
can achieve expected resolution is denoted the S N R threshold.
1.2
R o le o f t h e C o v a r ia n c e M a tr ix
T h e u n iv a ria te G au ssian or norm al d istrib u tio n occupies a cen tral position in th e s ta tis tic a l th eo ry of analy zin g ran d o m variables since th is s ta tis tic a l m odel is s u it able fo r such a large n u m b er of cases w hen m u ltip le m easu rem en ts a re tre a te d . T his s itu a tio n is even m ore pronounced in m u ltiv a ria te analysis due to th e p aucity of a n a ly tic a lly tra c ta b le m u ltiv a ria te d istrib u tio n s — one n o tab le exception being th e m u ltiv a ria te G aussian d istrib u tio n . C onsider a com plex G aussian random vec to r x € C Mxl w ith m ean v ector p = E[x] a n d nonsingular H e rm itia n covariance m a trix R - E [ ( x — p ) ( x — p ) H]. T h e n th e density of th e x is
/(X ) = irM d e t(i?) 6XP H * ~ ~ /*)] (1-4)
For G au ssian ra n d o m vectors, m ean vectors an d covariance m atric e s are sufficient s ta tis tic s , which can be used to uniquely d e te rm in e th e d istrib u tio n d en sity func tions. F u rth e r, in m an y p ra c tic a l situ a tio n s, m ean vectors of m u ltiv a ria te random processes are e ith e r zero or can be rem oved from th e d a ta vectors, th u s covariance m a tric e s becom e th e only s ta tis tic s to b e e s tim a te d for G aussian ran d o m processes.
Sim ilarly, covariance m atric e s also play a c e n tral role in high resolution algo rith m s developed for sp e c tra l e s tim a tio n a n d a rra y processing since we usually use th e G aussian d istrib u tio n to m odel ran d o m vectors o b ta in e d in th e s e applications. As a m a tte r of fact, m an y of th ese high resolution algorithm s req u ire a tw o step p rocedure. F irs t an e s tim a te of th e covariance m a trix is o b ta in e d from th e given d a ta set; th e n th e sp e c tra l analysis algorithm s a re applied. T h e results of th e second ste p a re highly sensitive to th e q u a lity of th e resu lts of th e first ste p . For s ta tio n a ry ran d o m processes, th e e s tim a te of th e covariance m a trix w ith required sta b ility can be o b tain e d by th e sam p le covariance m a trix using a sufficient long o b servation tim e period. In th e o ry th ese e stim a te s (sam ple covariance m atrices) a p p ro a c h th e tru e covariance m atric e s w ith p ro b ab ility one w hen th e observation tim e is infinitely long, assum ing th a t th e ran d o m processes are sta tio n a ry a n d er- godic. How ever, in p rac tic e th e observation tim e is usu ally lim ite d since e ith e r th e ran d o m processes can only be tre a te d as sta tio n a ry in a lim ite d tim e in terval or th e available d a ta sam ples from th e ran d o m process are finite. T h u s a p rac tic a l pro b lem is to find th e b e st possible e s tim a to r for a covariance m a trix from a finite observation.
1 .3
S u m m a r y o f O th e r W o rk
1.3.1
S a m p le C ovarian ce M a trix
T h e sam p le covariance m a trix has been w idely used in th e e stim a tio n of th e co- variance m a trix of m u ltiv a ria te ran d o m processes. It is show n th a t th e sam ple covariance m a trix is th e m ax im u m likelihood (M L) e s tim a te of th e tru e covariance for i.i.d. G aussian m u ltiv a ria te ran d o m vectors assum ing t h a t th e tru e covariance m a trix does n o t have any specific s tru c tu re o th e r th a n s y m m e try [6]. T h is e s ti m a to r is no longer o p tim a l w hen ad d itio n a l a p rio ri in fo rm a tio n regarding th e s tr u c tu r e of th e covariance m a trix or signal sources a re available.
C H A P T E R 1. I N T R O D U C T I O N
1 .3 .2
M a x im u m L ik elih o o d E s tim a te o f th e S tr u c tu re d
C ovarian ce M a trix
In th e case of a uniform lin e a r array a n d un co rrelated sources, we have additional
a p rio ri s tru c tu ra l in fo rm atio n of th e covariance m a trix . O ne such piece of infor
m a tio n is th e T oeplitz s tr u c tu r e of th e covariance m a trix . To exploit this in form a tio n , B u rg et al [7] derived an ite ra tiv e alg o rith m to co m p u te th e ML e stim a te of a T o ep litz co n stra in e d covariance m a trix .
1.3.3
P e r io d ic ity C o n str a in e d E stim a to r
A n o th e r piece of s tru c tu ra l in fo rm atio n is th e periodicity of th e sp a tia l covariance fu n ctio n for a uniform lin e a r array a n d u n c o rre la ted sources. Ziskind and Wax [8] proposed a schem e to ex p lo it th e p e rio d icity p ro p e rty of th e sp a tia l covariance fu n ctio n e x iste n t in these scenarios. W h e n ap p ly in g th e perio d icity con strain ed es tim a tio n m e th o d in c o n ju n ctio n w ith th e high resolution a lg o rith m , th e resolution en h a n ce m en t can be achieved especially in th e th resh o ld region.
1 .3 .4
C o h eren t S ig n a l-su b sp a c e M e th o d
In th e scenarios of tem p o rally w ide-band sources W ang a n d K aveh [9] proposed th e C o h eren t Signal-subspace M eth o d (C SM ) to m ore effectively exploit larger source ba n d w id th s. In CSM , focusing m atric e s are used to reduce m u ltip le narrow -band covariance e stim a te s m ad e over th e receiver b a n d w id th io a single (hopefully) focussed covariance m a trix . CSM offers a significant im provem ent over th e inco h e re n t m e th o d s in th e w ide-band se ttin g .
1.3.5
S te e r e d C ovarian ce M a trix
Also in th e w ide-band settin g , K rolik and Sw ingler [10] proposed to use steered covariance m atrices. It is shown th a t th e steered covariance m a trix has th e a d vantage th a t it can be e stim a te d w ith m uch g re a te r s ta tis tic a l s ta b ility for a given sam ple size an d ex h ib its lower SN R th resh o ld .
1.4
C o n tr ib u tio n s o f t h e D is s e r t a t io n
In th is d isse rta tio n a novel covariance m a trix e s tim a to r [ l i , 12] will b e p resented. T his e s tim a to r takes advantages of cross correlations e x isten t am ong d istin c t sets of ran d o m vectors draw n from different processes. F or exam ple, co rrelatio n s can be th o se betw een o u tp u ts from different su b array s in sp atial sm o o th in g , or b e tween sn ap sh o ts of th e a rra y o u tp u ts or even betw een tim e sliced vectors from a d a ta sequence in sp e c tra l e stim a tio n . A pplying this m eth o d In co n ju n ctio n w ith a given high reso lu tio n alg o rith m resu lts in en h an ced perform ance im p rovem ent for th e sp e c tra l e stim a to r. T h e im proved covar'ance e s tim a to r is derived from th e theorem on th e least squares lin ear p red ic tio n of m u ltiv a ria te ran d o m processes; it is fu rth e r sim plified for covariance-based (an d especially su b space-based) high resolution sp e ctral e stim a to rs. It has ap p licatio n s in th e following sp e c tra l analysis and a rra y processing areas, and som e of th ese ap p licatio n s will b e stu d ie d in d e ta il later:
1. S p a tia l sm o o th in g in DOA e stim a tio n .
2. S p e c tru m e s tim a tio n in tim e sequence analysis. 3. E x p lo itin g te m p o ra l correlations in a rra y processing •1. B eam form ing
C H A P T E R 1. I N T R O D U C T I O N 7
6. Focusing tra n sfo rm a tio n design in CSM
In a d d itio n to th e new covariance m a trix e s tim a to r we will also p rese n t in this d isse rta tio n som e o th e r w ork we have done in th e areas of sp ectral analysis and array processing.
1.5
O r g a n iz a tio n o f t h e D is s e r t a t io n
T his d isse rta tio n is organized as follows:
• C hapter 2: provides a h isto ric review of m a jo r classical and m o d ern sp ectral e s tim a tio n m eth o d s. A unified signal an d noise m odel is e stab lish ed for th e cases of b o th array processing an d sp e ctral analysis.
• C hapter 3: presen ts a covariance m a trix e stim a tio n m eth o d w hich can fully e x ploit th e cross co rrelations of m u ltiv a ria te ran d o m vectors. T h e e s tim a tc r is fu rth e r sim plified for covariance-based high resolution sp ectral e stim a to rs, w here signal subspace is of d o m in a n t im p o rta n c e . A sy m p to tically th e re is no difference betw een th e new covariance m a trix e stim a to r an d o th e r es tim a to rs, b u t th e proposed m e th o d shows its ad v antage over th e sam ple covariance m eth o d w h en th e observation tim e is sh o rt. As a consequence of app ly in g th e new m e th o d , th e perfo rm an ce of th e following high resolution alg o rith m s can be enh an ced . Even th o u g h for a specific sp e ctral e s tim a to r a b e tte r covariance m a trix e stim a te will yield sim ultaneously sm aller e rro r vari ance, lower S N R th re sh o ld , a n d higher resolution, any one of th ese m easure's justifies th e m eth o d .
• C hapter 4'• applies th e th eo ry developed in C h a p te r 3 to different, areas of a rra y processing a n d sp e c tra l analysis. For each ap p lic atio n th e th e o ry of C h a p te r 3 is a d a p te d to th e specific problem ; sim u latio n resu lts verify th e achievem ent of th e e x p e cte d goal.
• C hapter 5: presents a tra n sfo rm a tio n m a trix m eth o d based on th e vector p-norm w hich generalizes th e com m only used F robenius norm m eth o d . • C hapter 6: is focussed on th e applicatio n of th e m o d ern a rra y processing
m ethodology to reflection seism ic signal processing. E specially th e velocity estim at ion problem has been studied.
C h a p te r 2
H isto r ic R e v ie w o f S p e c tr a l
E s tim a to r s
2 .1
W h y H ig h R e s o lu t io n S p e c tr a l E s tim a to r s
2.1 .1
C la ssica l S p e c tr a l E stim a to rs
C lassical dig ital sp e c tra l e s tim a tio n m e th o d s are usually based on th e D iscrete F ourier T ransform (D F T ) or th e Wiener-Khintchine th eorem . Suppose we o b serve a s ta tio n a ry ra n d o m process for a v e ry long tim e, so t h a t we o b tain a tim e series y ( n ) for n = 1, 2 , . . . , N , w here N is very large. T h e n th e associated power
s p e c tru m can be o b ta in e d by th e p eriodogram
P (u ,) = (2 .1 )
w here Y (uj) is th e D F T of y ( n )
n = 1
Today Y(u>) is c o m p u ted very rap id ly by m eans of th e C ooley-T ukey fast F ourier tran sfo rm (F F T ) [13].
W iener in his 1930 p a p e r [14] gave th e following m e th o d , w hich was sta n d a rd until th e work of T ukey in 1949 [15]. W ie n e r’s m e th o d was in te n d e d for very long tim e series. It consisted of co m p u tin g th e autc correlation fu n ctio n as th e tim e average
r(*)= 77 Y, y*(n)y(n +
’ n
w ith th e asterisk s ta r * d e n o tin g com plex co n ju g ate, for — p < k < p, w here p less th a n th e d a ta len g th N , and th e n c o m p u tin g th e pow er s p e c tru m P(tu) as th e Fourier tra n sfo rm
OO
P ( u , ) =
£ r(t)e-''"*,
(2.2)
k ——oo
w here r ( k ) is tru n c a te d a fte r p lags. T h is F ourier tra n sfo rm rela tio n sh ip betw een th e a u to c o rre la tio n and th e power sp e c tru m is now called th e Wiener-Khint chi ne th eorem .
2.1 .2
L im ita tio n o f C la ssica l S p e c tr a l E stim a to r s
In sp ite of th e c o m p u ta tio n a l efficiency of th e F F T , reso lu tio n s of th ese spec tra l e stim a to rs a re lim ited to th e reciprocal of te m p o ra l or s p a tia l a p e rtu re s of th e recorded d a ta . W hen high resolution sp e c tra l analysis is req u ired , th e d a ta record len g th sh o u ld be long in th e case of tim e series analysis, or th e sensor a rra y should cover a larg e sp a tia l range in th e case of a rra y processing. As we m en tio n ed earlier, in p ra c tic a l rea lity large te m p o ra l or sp a tia i a p e rtu re m ay b e difficult to o b tain . For e x am p le th e sizes of sensor array s m ay n o t be easily ex te n d e d ; th u s
C H A P T E R 2. H I S T O R I C R E V I E W 11
th e s p a tia l a p e rtu re is often m uch sm aller th a n req u ired to o b ta in desired high resolution using th e conventional m eth o d s. F u rth e r, th e resolution of a spectral e s tim a to r is in som e sense closely rela te d to th e cap ab ility of d etection of signals by th e e s tim a to r. Increasing signal d ete c tio n cap ab ility is a com m on need in rad a r, sonar a n d m a n y o th er a rra y processing applications. T h is need has been m o tiv a t ing th e developm ent o f so-called high resolution sp e ctral e stim a to rs. Discussions of th e d isse rta tio n will b e focused on recently developed high resolution spectral estim a to rs.
2.2
D a t a M o d e lin g a n d A s s u m p tio n s
Since s p e c tra l e s tim a tio n an d a rra y processing are equivalent in m any asp ects, the signal m odels u sed in b o th ap p licatio n s can b e ch aracterized in a unified approach. In th is section th e signal m odel for b earin g (D O A ) e stim a tio n problem in array processing is considered first; la te r it is show n th a t th is signal m odel will be equally ap p licab le to s p e c tra l e stim a tio n . To sim plify th e co m p lex ity of fo rm ulation, we m o stly consider c ne-dim ensional e s tim a tio n problem s from e ith e r tem p o rally or sp a tia lly un ifo rm ly sa m p le d d a ta . T h e gen eralizatio n of th e one-dim ensional and uniform ly sa m p le d case to m u lti-dim ensional and non-unir’o rm ly sam pled cases is a s tra ig h tfo rw a rd ex ten sio n of th e sim plified case.
2.2 .1 A rray P r o c e ssin g
C onsider a u niform lin e a r sensor a rra y com posed of M om n i-d irectio n al sensors im m ersed in a hom ogeneous dispersive m ed iu m w ith wave p ro p ag a tio n velocity c. A ssum e t h a t d n a rro w -b an d sources (d < M ) cen tered a t frequency u>0 are
im p inging on th e a rra y from directio n s 6 1,6 2, ■■■ ,&d a n d are asso ciated w ith signals s i ( t ) , ^2(1), • • ■, sd(t). A ssum e fu rth e r th a t th e signals e m itte d by th e sources are
vector is specified by th e lin ear com bination of steering vectors , w hich are m u tu a lly linearly indep en d en t
d
*(*) = + w (t )> (2-3)
It=x
where w ( t ) is th e a d d itiv e noise v ecto r an d a{0k) is th e steering v ecto r for th e k t h
signal, which is given by
a(0k) = [1, eJ'u'oTfc, ej2w°T k (2.4)
where r k = A sin 0k/ c w ith c den o tin g th e p ro p a g a tio n speed of th e signal in th e m edia, w ith A rep re sen tin g sensor spacing satisfying co n d itio n A < A/2, a n d w ith A d en o tin g th e w avelength. We use su p e rsc rip ts T a n d — 1 to rep re sen t respectively
m a trix tra n sp o se and inverse o p erations. a{0) is also called th e direction vector or array manifold. W e use th e sym bol 0 w ith o u t a su b scrip t to rep re sen t a possible
d irectio n of a rriv al a n d th e su b scrip ted sym bol 0k, k = 1,2, . . . , d, to rep re sen t th e
tru e directio n of arrival in th e noise-free d a ta . A lte rn a tiv e ly eq u a tio n (2.3) can b e w ritte n by a v e c to r-m a trix form
x ( t ) = A s ( t ) + w ( t ) , (2.5)
w here A = [a(#i), a i Of ) , . . . , a(0d)] is th e steering matrix, an d s( t) — [s i(t), ^ ( t ) , . . . , Sd(t)]T. T h e a d d itiv e noise is assum ed to b e a sta tio n a ry , zero-m ean ran d o m process th a t is te m p o ra lly an d sp a tia lly w h ite a n d u n c o rre la ted w ith th e signals. W ith th ese assu m p tio n s we get th e covariance m a tr ix of th e a rra y o u tp u t vectors
R = E [ x ( i ) x H{t)} = A P A h + a2I , (2.6)
w here P = E [ s ( t ) s H(t)] is th e signal covariance m at r i x a n d <r2I = E [ w ( t ) w H( t )] is th e noise covariance matrix. W h en P is a d x d diagonal m a trix , th e signals
C H A P T E R 2. H I S T O R I C R E V I E W 13
are uncorrelated; w hen P is nonsingular an d nondiagonal, th e signals are partially
correlated; a n d w hen P is singular, a t least two signals are coherent, ( th a t is, a t
least one source is a scaled and delayed version of th e o th er source(s)).
In th e above, t should b e replaced w ith t,, i = 1 ,2, . . . , N w hen th e array o u tp u t
is d ig itized, a n d N denotes th e n u m b er of a rra y o u tp u t vectors available. We also call N th e n u m b e r of snapshots .
2 .2 .2
S p e c tr a l E s tim a tio n
Let a scalar record of d a ta sequence y ( t , ) , i — 1 ,2 , be com posed of a uniform ly sam p led sum of d cisoid signals. N = N r — M + I d a ta vectors x ( t t) £ C M xl can b e form ed by
d
x { t i ) = [y(ti),y{ti + l ) , . . . , y ( t i + M - 1)]T = ^ 2 a ( u k) s k(ti) + w(ti), (2.7)
k= 1
w here a( uk) a n d w( ti) have th e sim ilar definition as in th e case of a rra y processing, i.e. a(o>fc) = [1, e^Wk, e ^2wk, . . . , e ^ M~ ^ Wk]T w ith wjt den o tin g th e frequency in r a d / s
of th e fcth cisoid, a n d w( ti) € C Mxl is th e a d d itiv e noise vector form ed in th e sam e way as x ( 't). T h e m a trix form of e q u a tio n (2.7) is given by
x(t i) = A s ( U ) + w(t{), (2.8)
w here A = [a(oq), a(u>2) , . . •, a(w<i)] is th e signal model matrix, an d s(t,) = [.si(fj), s2(ifi), • • • ,Sd(U)]T . W ith th e sam e assu m p tio n s on th e a d d itiv e noise, we have
R = E[ x { U ) x H(ti)] = A P A h + a2I . (2.9)
C o m p a rin g eq u atio n s in th e above tw o subsections, we conclude th e I)O A e s ti m a tio n a n d frequency e s tim a tio n problem s are equivalent in th e form of signal
m odeling and fo rm u latio n . D ue to th e equivalence, e ith e r sp ectral e s tim a tio n or the D O A e stim a tio n or b o th are used in pro b lem form ulation.
2 .2 .3
S u m m a ry o f A ssu m p tio n s
We su m m a riz e a ssu m p tio n s often im posed in a rra y processing and s p e c tra l a n a ly sis, a n d /o r in th e previous signal a n d noise m odeling:
• T h e n u m b er of signals is less th a n th e n u m b e r of sensors, nam ely, d < M . • T h e g eom etry of th e sensor a rra y is such th a t a{0k) for k = 1 , 2 , . . . , d are
linearly in d ep e n d e n t. For uniform linear arrays a n d u niform ly sam p led tim e sequences, this assu m p tio n is a u to m a tic a lly satisfied due to th e V anderm onde s tru c tu re of th e A m a trix .
• T h e noises w(t{) a re in d ep e n d e n t of signals Sk(U) a n d of each o th e r, an d are id en tically d istrib u te d com plex, zero-m ean, G aussian vectors w ith covariance m a trix a2 1, w here a2 is an unknow n scalar.
• T h e sensor spacing A is less th a n or eq ual to half of th e sm allest w avelength to avoid sp a tia l aliasing effects; nam ely, A < Am,n /2 .
2 .2 .4
S a m p le C ovarian ce M a tr ix
In a p rac tic a l d ig ital sy ste m a tru e covariance m a trix of a rra y o u tp u t vectors is not available; on ly an e s tim a te is o b ta in e d using N sam p led o u tp u t vectors x(t{):
£ = 7F (2. i0)
JV i = l
R is te rm e d th e sample covariance ma t ri x of x(t{), w here we use th e ca ret sign
C H A P T E R 2. H I S T O R I C R E V I E W 15
2 .3
H ig h R e s o lu t io n S p e c tr a l E s tim a to r s
W ith th e e s tim a te of th e a rra y o u tp u t covariance m a trix given, we a re able to p resen t so-called high resolution s p e c tra l e stim a to rs in te rm s of th e e stim a te d covariance m a trix . In his recent p a p e r, K irlin [16] fo rm u la ted a unified p re se n ta tio n for several recen tly developed high resolution a lg o rith m s based on th e eigen-decom position of R . It is well know n th a t th e covariance m a trix R is pos itive (sem i-)definite, th u s it can b e w ritte n in term s of its ordered eigenvalues Ai > A2 > • • • > Am and asso ciated eigenvectors u,-, i = 1,2, . . . , M:
M
R = v k V H = £ X i V iv[r ,
« '= l
( 2 . 1 1 )
w here V is an M x M u n ita ry m a trix com posed of M eigenvectors
V = [ui , u2, . . . , v m],
and A is a diagonal m a trix form ed w ith elem ents of A;
A = diag(A i,A2, . . . , A M).
S im ilarly th e inverse of R can be w ritte n as:
M
R - ' = Y i A r W t'=l
(2.12)
B ased on th e eigen-decom position of R , we can p a rtitio n signal and noise subspace in te rm s of eigenvectors, or
R = [ V sVn]
' k
o ' \ V aH. 0 An . X H .
w here V, com posed of th e first d eigenvectors form s th e signal subspace and Vn com posed of th e rem ain in g eigenvectors form s th e noise subspace, respectively. S im ilarly A., is th e eigenvalues associated w ith signal su bspace a n d An associated w ith th e noise subspace.
T h e c a re t sign m ay b e neglected for n o ta tio n b re v ity in th e rem aining p a rt of th e d isse rta tio n if it can b e inferred from co n tex t or if th e e stim a te does not need to be d istin g u ish ed from th e tru e value.
2.3.1
C o n v e n tio n a l B ea m fo rm in g
T h e conventional m e th o d of m ap p in g th e m o n o ch ro m atic field pow er as a function of signal a rriv in g angle is th e b eam form ing o p e r a t e . v h ic h em ploys a proced u re known as d elay-and-sum processing to steer a b e a m in a p a rtic u la r directio n . G iven m easu rem en ts from a sensor a rra y and th e field of view 0 G 0 , th e b eam fo rm er
scans th e region of in te re sts in w hich sources m ay be p re se n t a n d calcu lates a rra y o u tp u t pow er
a » ( O ) R a (0)
“ ' a " (0) a(0) '
*
*
T h e peaks of th e s p e c tru m rep resen t th e e stim a te s of D O A . T h e conventional beam fo rm in g is c o m p u ta tio n a lly sim ple, b u t it has poor p erfo rm a n c e a t low SN R or in presence of m u ltip le sources.
2 .3 .2
M a x im u m E n trop y M e th o d (M E M )
T h e Wiener-Khintchine th eo re m is n ot s u ita b le for s h o rt d a ta records since th e a u to c o rre la tio n fu n ctio n is tru n c a te d after p lags. T his re su lts in a sm o o th in g effect in th e frequency d om ain. W hen th e p a u to c o rre latio n sequence r ( l ) , r ( 2 ) , . . . , r(p)
C H A P T E R 2. H I S T O R I C R E V I E W 17
is assu m ed know n, th e n a question m ay be posed as to how th e rem ain in g unknown lags r ( p + l ) , r ( p + 2) , . . . should b e specified in ste a d of a rb itra rily se ttin g th em
to zero. M otivated by th e results of in fo rm atio n theory, B urg [17] argued th a t th e e x tra p o la tio n should b e m ad e in such a way as to m axim ize th e e n tro p y of th e tim e series ch ara c te riz e d by th e e x tra p o la te d a u to c o rre latio n sequence. T his tim e series w ould th e n b e m o st random , in an e n tro p y sense, of all series t h a t have tin ’ know n first p lags. H e proposed th e m ax im u m e n tro p y m eth o d which begins with m axim izing th e e n tro p y r a te of a G aussian ran d o m process
f l n P M E M { f ) d f (‘2.15)
J - v r/ 2
su b je c t to th e c o n stra in ts
F P M E M ( f ) e j2*Jkd f = r ( k ) (2.16)
J —tt/2
over 0 < k < p. T h e so lu tio n found by th e L agrange m u ltip lie r techniques can be w ritte n as follows:
i T d -1i
w / ) - i u ^ - „ , (; ) i r <*•»>
w here 1^ = ( 1 , 0 , . . . , 0 ,0 ) a n d R is th e T oeplitz au to c o rre latio n m a trix of know n lags.
I t h a s been show n t h a t for G aussian ran d o m processes a n d a know n a u to c o r re la tio n sequence of un ifo rm spacing th e M EM s p e c tru m is equivalent to t h a t of o rd er p auto-regressive (A R ) m odel based on linear p red ic tio n th eo ry [3]. A good review on th e rela tio n sh ip of M EM an d A R m odel can be found in P apoulis [18].
2.3.3
M in im u m V arian ce M e th o d (M V )
T h e m in im u m variance m eth o d was developed by C apon [19]; it is also know n as C a p o n ’s M axim um Likelihood m eth o d . In a filtering in te rp re ta tio n , th e M V e s tim a to r is an F in ite Im pulse R esponse (F IR ) d igital filter w hich m inim izes th e variance of th e a rra y o u tp u t (an d th u s interferences) u n d e r th e c o n stra in t t h a t its response to th e specified d irection is unity:
m in E [ ||tu Hx(< )||2] su b je c t to w Ha (0) = 1. (2.18)
T h e o p tim a l solution, easily found using Lagrange m u ltip lie r techniques, is
2 .3 .4
P isa r en k o M e th o d
T h e P isaren k o m e th o d [20] is th e first a lg o rith m to em ploy a subspace of th e e stim a te d covariance m a trix . In th e P isaren k o m e th o d , we first c alcu late th e
(d -f 1) x (d + 1) covariance m a trix a n d vj..\ i, th e eigenvector associated w ith
th e sm allest eigenvalue. It L easy to show th a t first d eigenvectors form th e sam e vector su bspace as th e steering m a trix A . Since tv+ 1 is o rth o g o n al to all th e first d eigenvectors a n d th e signal subspace, it is also o rth o g o n al to all directio n v ecto r a(0) a t th e tru e DO A. T h e o rth o g o n a lity re la tio n sh ip is given by
a H(0k)vd+i = 0 for k = 1,2, . . . , d. (2.2 0)
T hen th e e x trem a-search in g a lg o rith m of th e P isaren k o m e th o d is given by
C H A P T E R 2. H I S T O R I C R E V I E W 19
Ideally, a t th e tru e signal incident d irectio n , Ppisarenko{0) has a value of infinity. H ow ever, since we only have e s tim a te s of th e eigenvector iv+ ii th e orthogonality rela tio n sh ip does not hold any m ore. In ste a d th e cosine of th e angle betw een (i(0*),
k = 1,2, . . . , d , and v d+i will b e p ro b ab ly close to zero, or Ppiaarmko(0) will have
peak s a t th e e stim a te s of DO A.
For a uniform lin e a r a rra y th e P isarenko m eth o d can be represented in a poly n o m ial-ro o tin g fo rm a t
d
D ( z ) = a H( z ) vd + 1 = r0 ~ n ~ - 1 ), (2.22)
i = l
w h ere z — eJ w°Asm0l c, Xhe d irectio n s of arrival are found from th e angles of th e d ro o ts of th e polynom ial.
2 .3 .5
M U S IC
T h e M U SIC or M u ltip le S ignal C lassification algorithm was first developed by S c h m id t [21] a n d B ienvenu et al [2 2]. T h e M U SIC alg o rith m can b e view ed as an
ex ten sio n of th e P isaren k o m e th o d since it uses a covariance m a trix o f size M x A/, w here M c an be any n u m b er as long as M > d is satisfied. T h e M U SIC alg o rith m , to o , em ploys th e o rth o g o n a lity rela tio n sh ip of th e noise subspace to steerin g vector a t t h e tru e signal directio n s 0^ for k = 1,2, . . . , d or
a H(0k)Vn = 0. (2.23)
T h u s th e M U SIC pow er sp e c tru m is w ritte n as
w ith th e peaks of th e sp e ctru m rep resen tin g th e DOA e stim a te s. Like o th e r subspace-based m eth o d s, M U SIC has th e ad v antage th a t it can b e ap p lied to a sensor array of an a rb itra ry geom etry, w hich is know n or c a lib ra te d . However it will to ta lly fail w hen coherent signals are present.
2.3.6
R o o t M U S IC
R oot M U SIC, first suggested in [23], is a v ariatio n of M U SIC , a n d is ap p licab le to uniform linear sensor arrays. T h e polynom ial of R o o t M U SIC is given by
Pr m(z) = aH( z ) V nVnHa( z )
M -l
= r0 n (1 - r « '2 - 1 )( 1 ~ r *z)- (2.25)
i=i
T h ere a re 2( M — 1) ro o ts associated w ith th e polynom ial Pr m' th e d roots w ith i
th e larg est am p litu d e in sid e u n it circle are chosen as signal zeros, a n d th e ir angles are th e DOA e stim a tes. T h e ad v an tag e of R o o t M U SIC over S p e c tra l M U SIC can be u n derstood by considering th e effects of a n e rro r Ar,-. A ssum ing th e erro r A r; is a ran d o m variable w ith m ean r,-, it can b e decom posed in to ra d ia l a n d ta n g e n t com ponents. Since th e D O A e stim a te s ta k e angels of th e d ro o ts w ith th e larg est a m p litu d e inside u n it circle, th e ra d ia l co m ponent of A?’,' will cause no erro r in th e DOA e stim a te by R ^ot M U SIC. How ever, such ra d ia l errors do afreet th e sp e c tru m by S p ectral M USIC. In fac t, th e sh arp n ess of th e sp e c tra l peak s is d e te rm in e d by th e rad ia l com ponent of Ar,- or th e ra d ii! d ista n c e betw een r,- a n d th e u n it circle. T his is p a rtic u la rly critical for closely spaced roots as it m ay resu lt in only one peak causing an a p p a re n t loss in resolution. So S p e c tra l m e th o d s alw ays have less resolution com pared to R o o t form s.
C H A P T E R 2. H I S T O R I C R E V I E W 21
2 .3 .7
M in im u m N o rm M e th o d
T h t vlinim um N orm m e th o d [24] is also a subspace-based alg o rith m . It uiuls th e
M X 1 v ector dj^N w ith a u n it first elem en t, which is en tirely n th e noise subspace and h as th e m in im u m 2-norm . W e p a rtitio n th e signal and noise subspace vectors as follows
V„ =
gH
v :
c*
V'y n
w here g H and c11 are th e first rows of Vs and Vn , respectively. Since lies in th e ran g e of Vn , dj^N will b e orthogonal to th e colum ns of Vs, i.e.,
VsHdMN = 0. (2.26)
Solving e q u a tio n (2.26) v ia th e le a st squares technique, we have
dMN = 1
JOB-. 1 - g H g
(2.27)
A lte rn a tiv e ly d/vfjv can b e represen ted u sing noise eigenvectors as follows
dMN =
1_ CHC (2.28)
r or a u niform lin e a r a rra y th e polynom ial - ooting schem e for th e DO A e stim a tio n using th e M in im u m N orm m eth o d is given by
M - l
D ( z ) = aH ( z ) d MN = I I ( i - r i z ~ l ).
»=i
Ideally D ( z) has d zeros on th e u n it circle and th e rem ain in g zeros lie inside th e unit circle. In th e case of noisy d a ta th e d ro o ts closest to th e u n it circle are chosen as th e signal roots; th e ir angles are th e D O A e stim a tes.
T h e extrem a-search in g form of th e M inim um N o rm m e th o d is given by
(2.30)
whose peaks rep re sen t th e D O A e stim a tes. Pm n( 0 ) =
2.3.8
E S P R IT and M a trix P e n c il
E S P R IT (E stim a tio n of Signal P a ra m e te r v ia R o ta tio n Invariance T echnique) [25] is a technique th a t m ay be used to e stim a te frequency or d ire c tio n of arrival. As com pared to M U SIC, E S P R IT reta in s m ost of th e essen tial fea tu re s of th e a rb itra ry array of sensors, b u t it achieves a significant re d u c tio n in c o m p u ta tio n com plexity of M U SIC associated w ith ex trem a-search in g proced u res by im p o sin g a c o n stra in t on th e s tru c tu re of th e sensor array. In o th e r w ords, th e re should b e tw o id en tical su b arrays available; one is a know n-displacem ent 6 (no ro ta tio n ) version of th e o th er. T he o u tp u t vectors from th e two su b array s a re
x (t ) = A xv(t) + n x (t) = A s ( t ) + n x(f); (2.31)
y{t) = A ys{t) + n„(f) = AQ s{ t ) + n y{t), (2.32)
w here th e m a trix $ is a diagonal d x d m a trix of th e p h ase delays betw een th e doublet sensors for th e d wavefronts:
C H A P T E R 2. H I S T O R I C R E V I E W 23
L et E sx and E sy be th e signal subspace o b tain e d from su b a rra y o u tp u t vector .r(t) a n d y ( t ) , respectively. W e can form a m a trix pencil as follows:
Vay$ = Vsx. (2.34)
W h e n e stim a te s of a n d Vsy are available, $ can be d e te rm in e d using a least, sq u ares or a to ta l le a st squares approach. T h e m eth o d using a least squares a p p ro ac h to e s tim a te $ is referred to as L east Squares (LS) E S P R IT . A sym ptotically, L east Squares and T o ta l L east Squares (T L S) E S P R IT will have th e sam e perfor m a n c e [26]. F or LS E S P R IT
* = ( V ? y Vs yY l VsyVsx. (2.35)
E S P R IT can also b e c h a ra c te riz e d by th e m a trix pencil from linear algebra theory. It is shown [27] th a t th e D O A e stim a tio n problem can be view ed as th e generalized eigenvalue problem .
2 .3 .9
M a x im u m L ik elih ood (M L ) M e th o d
M L e stim a to rs possess salient p ro p erties of a sy m p to tic unbiasness and a sy m p to tic efficiency. A ssum e t h a t th e noise vector w ( t ) is a sta tio n a ry and ergodic com plex v alu ed G aussian process of zero m ean an d covariance m a trix a2I , w here a2 is an unknow n scalar. T h e noise sam ples w ( U ) , i = 1 , 2 , . . . , N a re s ta tistic a lly in d ep e n d e n t. If we p ack N in d e p e n d e n t sn ap sh o t vectors in to an M x N m a trix colum n by colum n
X = A S + W, (2.36)
w h e re X a n d W are th e M x N m atrices
W = [to(*i),tu(<2),...,u>(<Ar)], and S is th e d x N m a trix
S = [a(fi),jj(f2) , . . . , s ( f j v ) ] ,
th e n th e jo in t d istrib u tio n of th e sam p led d a ta is given by
f ( X \ 0 ) = n ^ - ^ :2- j j e x p { - [ x ( t ,) - A s ( Z i)]/ / [tr2/ ] "1^ ( ^ - A s ( i ; ) ] } . ( 2 . 3 7 )
T h e log likelihood, ignoring c o n sta n t te rm s, is given by
L = —N M l o g c r2 - - I £ ll*(*f) ~ M U ) \ \ 2. (2.38)
0-2 ,=i
To m axim ize L , we have to m axim ize th e log likelihood w ith resp e c t to all u n known p a ra m e te rs, th o u g h only th e D O A e s tim a te s a re of in te re st. To o b ta in th e com pressed likelihood [28], we m axim ize L w ith resp ect to a2 w ith A an d s fixed and get th e following
** = a In ^ ~ (2 -39)
S u b s titu tin g this resu lt back in to L a n d ignoring c o n sta n t te rm s, we get th e fol lowing m in im izatio n problem
C H A P T E R 2. H I S T O R I C R E V I E W 25
A pplying th e com pressed likelihood m eth o d once m ore by fixing A and m inim izing th e a b o v e fu n ctio n w ith respect to s yields
s (t i ) = ( A H A ) ~1A Hx ( t i). (2.41)
S u b s titu tin g s (t i ) into th e m in im izatio n fu n ctio n , we can o b ta in th e final log- likelihood fu n ctio n
L ( A ) = L ( 0 ) = t r ( P AR ), (2,12)
w here P A m— A ( A H A ^ A 11 is th e p ro jectio n o p e ra to r of A , R is th e sam ple covari ance m a trix e s tim a te d fro m th e given d a ta a n d tr[ ] is th e tra c e of th e bracketed m a trix .
In s p ite of t h e o p tim a lity of M L m e th o d , m axim izing th e log-likelihood func tion (2.42) involves m u ltiv a ria te no n lin ear p ro g ram m in g problem , which is com p u ta tio n a lly d e m a n d in g . F u rth e r, th e log-likelihood function is non-concave; th is req u ires a careful selection of th e in itia l value for th e m ax im izatio n . In c o n tra st, th e su b sp ace-b ased e x trem a-search in g m eth o d s in tro d u ced previously ju s t need th e one-dim ensional line search p ro ced u re, w hich is c o m p u ta tio n a lly m uch less com plex th a n ML. A lth o u g h several researchers [29] [30] [31] have proposed different fast n u m erica l a lg o rith m s, th e M L e s tim a to r a t this stag e is still m ore a theo retical s ta n d a rd to w hich th e su b o p tim a l e stim a to rs are com pared th a n a practical tool.
2 .4
C la s s ific a tio n s o f S p e c tr a l E s t im a t o r s
T h e re a re m an y classification m e th o d s for sp e ctral e stim a to rs dep en d in g on th e classification c rite ria used. T h e following m eth o d s are suggested.
2.4.1
C la ssifica tio n b y C h ro n o lo g y
Chronologically, we can roughly classify sp e ctral e s tim a to rs in to tw o categories:
1. Classical M ethods :
P eriodogram a n d th e Wiener-Khint chi ne th e o re m ty p ify th e classical m e th ods.
2. M odern M ethods:
All th e high resolution m eth o d s review ed p reviously in th is c h a p te r can be included in th is category.
2 .4.2
C la ssifica tio n b y M o d e l
S p ectral e stim a to rs m ay or m ay not b e m odel-based; th a t is, th e y m ay b e p a ra m etric or non -p aram etric:
1. P a ra m e tric M ethods:
In p a ra m e tric m eth o d s of sp e ctral analysis a m o d el is assu m e d in th e form u latio n of th e problem , a n d th e p a ra m e te rs of th e m odel a re e s tim a te d from th e lim ite d observation interval. A lth o u g h th e s e m odels can ta k e on a vari e ty of different form s, a ra tio n a l m odel w ith fin ite p a ra m e te rs is am ong th e m ost p o p u lar in th e co n te m p o ra ry sp e ctral e s tim a tio n . W h e n em ploying a ratio n al m odel, th e pow er sp e ctru m of th e signal processes has th e following form
p , jw) = b0 + bxe - i “ + • • • + 2
^ ’ 1 + a ie -J'w + . . . + ane - t t n- x)“
H erein, only m + n + 1 m odel p a ra m e te rs a re to b e d e te rm in e d . As we discussed previously, th e o re tic ally th e Wiener -Khint chi ne th eo re m requires infinitely long d a ta record. T h e n if th e re is on ly a sh o rt d a ta record available,